EPA-660/3-75-005

MARCH 1975

Ecological Research Series

MATHEMATICAL MODELING 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
MATHEMATICAL MODELING OF PHYTOPLANKTON

IN LAKE ONTARIO

1. MODEL DEVELOPMENT AND VERIFICATION

By

Robert V. Thomann

Dominic M. Di Toro

Richard P. Winfield

Donald J. O’Connor

Manhattan College, Riverdale, New York

Grant No. R800610

Program Element 1BA026

ROAP/Task No. 21AKP/15

Project Officer

William Richardson

Grosse Ile Laboratory

National Environmental Research Center

Grosse Ile, Michigan 48138

NATIONAL ENVIRONMENTAL RESEARCH CENTER

OFFICE OF RESEARCH AND DEVELOPMENT

U.S. ENVIRONMENTAL PROTECTION AGENCY

CORVALLIS, OREGON 97330

EPA-660/3-75-005

March 1975
ii

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 reviewed

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

interactions and kinetic behavior of the various components of each subsystem. A more

detailed, three-dimensional model is then used to describe open lake and nearshore

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

nonlinear differential equations must be solved while for the three-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

temperature. 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 reduction

programs on Lake Ontario.

This report was submitted in partial fulfillment of Grant Number R800610 by the

Environmental Engineering and Science Program, Manhattan College, Riverdale, New

York, under the sponsorship of the U.S. Environmental Protection Agency. Work was

completed as of October 1974.
iii

CONTENTS

Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii 

List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v 

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vi 

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x 

1. Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

2. Recommendations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

3. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

3.1 Purpose of Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

3.2 Scope of Research . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

4. Lake Ontario - Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

4.1 Morphometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

4.2 Hydrology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

4.3 Nutrients Inputs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

4.4 Available Lake Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

5. Lake Ontario - Model Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

5.1 Theoretical Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

5.2 Phytoplankton Chlorophyll . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

5.3 Nitrogen System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

5.4 Phosphorus System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5.5 Chemical Models - Dissolved Carbon Dioxide . . . . . . . . . . . . . . . . . . . . 20

6. Lake 1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

6.1 Basic Structure and Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

6.2 Model Calibration and Verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

6.3 Results of “Final” Verification Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . 50

7. Lake 2 Analytical Investigation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

7.2 Theoretical Considerations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

7.3 Asymptotic Population Growth Rates . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

7.4 Solutions for Specific Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

7.5 The Single Layer Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

7.6 The Effect of Increasing Growth Rate . . . . . . . . . . . . . . . . . . . . . . . . . . 73

7.7 The Effect of Increasing Settling Velocity . . . . . . . . . . . . . . . . . . . . . . . . 73

7.8 The Effect of Increasing Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

7.9 The Effect of Buoyancy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
iv

7.10 Preliminary Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

7.11 Summary of Mathematical Facts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

8. Preliminary Lake 3 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

8.1 Circulation and Dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

8.2 Validity of Assumed Circulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

8.3 Preliminary Phytoplankton Computations . . . . . . . . . . . . . . . . . . . . . . . . 92

Appendix A. Lake 1 - Input and Program Listing . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

A.1 Input Information . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

A.2 Morphometry and Hydrodynamic Regime . . . . . . . . . . . . . . . . . . . . . . . 99

A.3 Physical Exogenous Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

A.4 Biological and Chemical Exogenous Variables . . . . . . . . . . . . . . . . . . . 101

A.5 Chemical and Biological Endogenous Variables . . . . . . . . . . . . . . . . . . 106

References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
v

LIST OF TABLES

 Page

1. Flow Budget . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 

2. Nutrient Loadings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 

3. Nutrient Load Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 

4. CCIW Limnological Data Report Summary (Number of Cruises With

 Adequate Spatial Definition) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 

5. Thermodynamic Properties at 25°C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 

6. Thermodynamic Properties at 25°C . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 

7. Temperature Dependence Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 

8. Basic Physical Data of the Lake 1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 

9. Principal Parameter Values Used in Lake 1 Model Output (Shown in

 Figures 21-23) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 

10. Log10 (G/d + μ) for Exponentially Decreasing Growth Rate . . . . . . . . . . . . . . 71 

11. Monthly Heat Flux Values . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 

12. Regional to Lake 3 Segmentation Equivalence . . . . . . . . . . . . . . . . . . . . . . . 92 

A-1. Lake 1 Physical Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 

A-2. Lake 1 Mass Loadings and Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . 106 

A-3. System Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 
vi

LIST OF FIGURES

 Page

1. General basin map of Lake Ontario . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 

2. Typical CCIW cruise tracks, April 1969 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 

3. Typical output from data reduction - main lake average at 1 meter . . . . . . . . . 11 

4. Eutrophication modeling strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 

5. System diagram - Lake 1 model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 

6. Relationship of pH and the ratio, Alk/CT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 

7. Percent error in assuming [CO2 - Acy] = [CO2]. . . . . . . . . . . . . . . . . . . . . . . . . 25 

8. Computations for CO2 - H2O system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 

9. Major physical features included in Lake 1 model . . . . . . . . . . . . . . . . . . . . . . 33 

10. Effect of vertical mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 

11. Effect of phytoplankton settling velocity, Kmp = 10 μg/l . . . . . . . . . . . . . . . . . . . 38 

12. Effect of phytoplankton settling velocity Kmp = 1 μg/l . . . . . . . . . . . . . . . . . . . . 39 

13. Effect of phytoplankton settling velocity on nutrients and nutrient

limitation Kmp = 10 μg/l . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 

14. Effect of Michaelis constant for phosphorus on primary variables . . . . . . . . . . 42 

15. Effect of Michaelis constant for phosphorus on nutrient limitation

of nitrogen and phosphorus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 

16. Responses due to two phytoplankton growth formulations . . . . . . . . . . . . . . . 44 
vii

17. Sixteen-year computed output of phytoplankton chlorophyll . . . . . . . . . . . . . . 46 

18. Long-term behavior of model components under two phosphorus

Michaelis constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 

19. Preliminary comparison of model output and observed data . . . . . . . . . . . . . . 48 

20. Preliminary comparison of model output and observed data, continued . . . . . 49 

21. Lake 1 model verification, 0 -17 meters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 

22. Lake 1 model verification, continued, 0-17 meters . . . . . . . . . . . . . . . . . . . . . . 52 

23. Lake 1 model verification, 50-150 meters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 

24. Dynamics of the phytoplankton growth rate . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 

25. Dynamics of the phytoplankton death rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 

26. Dynamics of the phytoplankton net production . . . . . . . . . . . . . . . . . . . . . . . . 58 

27. Dynamics of nitrogen and phosphorus limitations . . . . . . . . . . . . . . . . . . . . . . 59 

28. Biological and hypolimnetic recycling of phosphorus . . . . . . . . . . . . . . . . . . . 60 

29. Lake Ontario phytoplankton biomass, horizontal averages, IFGYL,

 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1973 65 

30. Central Lake Erie phytoplankton biomass, four station log averages

Project Hypo, 1970 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 

31. Onondaga Lake phytoplankton biomass, two stations, 1969 . . . . . . . . . . . . . . 67 

32. St. Margaret’s Bay phytoplankton biomass, Station A, 1969 . . . . . . . . . . . . . . 68 

33. Eigenvalue equation solutions for constant growth layer and 

exponentially decreasing growth rate cases . . . . . . . . . . . . . . . . . . . . . . . . . . 70 

34. Settling velocity-growth rate interaction - the effect of dispersion.

Exponentially decreasing growth rate case . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 

35. Population growth rate versus maximum growth rate - the effect of

euphotic zone. The exponentially decreasing growth rate case . . . . . . . . . . . 74 
viii

36. Population growth rate versus settling velocity - the effect of euphotic

zone depth. The exponentially decreasing growth rate case . . . . . . . . . . . . . 75 

37. Population growth rate versus dispersion coefficient - the effect of

euphotic zone depth. The exponentially decreasing growth rate case . . . . . . 76 

38. The effect of negative sinking velocity. The exponentially decreasing

growth rate case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 

39. Application to lakes and oceans - the range of biologically relevant

dimensionless numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79 

40. Lake 3 segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 

41. Assumed summer circulation for Lake 3 model . . . . . . . . . . . . . . . . . . . . . . . . 84 

42. Assumed winter circulation for Lake 3 model . . . . . . . . . . . . . . . . . . . . . . . . . . 85 

43. Assumed horizontal and vertical dispersion coefficients for Lake 3 

model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 

44. Regional division of Lake Ontario by Lee (1972) . . . . . . . . . . . . . . . . . . . . . . . 89 

45. Comparison of computed temperatures from Lake 3 model with data

from Lee (1972), winter and fall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 

46. Comparison of computed temperatures from Lake 3 model with data

from Lee (1972), summer and spring . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 

47. Principal IFYGL stations and the Lake 3 grid . . . . . . . . . . . . . . . . . . . . . . . . . . 94 

48. Preliminary results of Lake 3, 0-4 meters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 

49. Preliminary comparison of Lake 3 model output to observed data at two

segments for June . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 

50. Data output flow diagram used for verification and display purposes . . . . . . . 97 

51. Three-dimensional plot of phytoplankton chlorophyll calculated from

Lake 3 model, June, 0-4 meters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 

A-1. Variation of volume with depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 
ix

A-2. Variation of area with depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 

A-3. Temporal variation of vertical dispersion coefficient - interface 1-2 . . . . . . . . 102 

A-4. Temporal variation of water temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 

A-5. Temporal variation of solar radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 

A-6. Temporal variation of photoperiod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 

A-7. Determination of background water clarity - Ke . . . . . . . . . . . . . . . . . . . . . . . 105 
x

ACKNOWLEDGMENTS

The cooperation of many International Field Year for the Great Lakes (IFYGL) investigators

who generally supplied their data for use in this report is acknowledged. Special thanks

are due to U.S. Environmental Protection Agency (USEPA) personnel at the Grosse Ile

Laboratory, including Messrs. William Richardson, Nelson Thomas, and Dr. Tudor Davies.

Finally, the excellent cooperation and assistance of the Canada Centre for Inland Waters

(CCIW) is also gratefully acknowledged.
1

CHAPTER 1

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 describing

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 interaction 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-wide model has been

established to permit the use of the model for

preliminary simulations of the effects of various

levels of nutrient reduction.
2

CHAPTER 2

RECOMMENDATIONS

Preliminary emphasis in this report is on the

scientific foundation and validity of the model,

principally the Lake 3 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 three-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 nearshore 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 nearshore

temperature and waste discharge.

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 phytoplankton 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 expanded 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 International

Field Year for the Great Lakes (IFYGL) effort.
3

CHAPTER 3

INTRODUCTION

3.1 PURPOSE OF RESEARCH

The overall purpose of this research is to

structure a mathematical modeling framework of

the major features of eutrophication in large

lakes. Lake Ontario, the subject of intensive

field work as part of the IFYGL is used as the

problem setting (Figure 1). The overall

objectives of the research include:

1. Determination of important interactions in

lake eutrophication.

2. Analysis of lake water quality and biological

responses to natural and man-made inputs.

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

3.2 SCOPE OF RESEARCH

The scope of this research is, therefore, lakewide and attention is directed primarily to the

behavior of the phytoplankton for the lakes as a

whole. The smallest scale considered in the

more advanced model (Lake 3) is on the order of

10-40 km. Detailed nearshore 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 research 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.
4

Figure 1. General basin map of Lake Ontario.
5

The primary thrust of the research is to provide

a more rational and scientifically defensible

structure for environmental 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 development and strategy, data

reduction and compilation, sensitivity analyses

and verification analyses. Simulations of

conditions under varying levels of nutrient

reduction are not included in this report and will

be part of a separate document.
6

CHAPTER 4

LAKE ONTARIO - BACKGROUND

With the Niagara River, the outflow of Lake Erie

and, hence, the cumulation of all the upper Great

Lakes outflow as its principle 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

corrosion (Hutchinson, 1957).

4.1 MORPHOMETRY

Lake Ontario is the smallest of the Great Lakes

in terms of surface area, 19,477 km2

 (7,520 mi2),

with a drainage area of 90,132 km2

 (34,880 mi2)

(Beeton and Chandler, 1963). 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 features (Beeton and Chandler, 1963;

Canada Centre for Inland Waters, personal

communication). Lake Ontario’s volume is 1,669

km3 (401 mi3), with a maximum depth of 244

meters (801 feet) (Canada Centre for Inland

Waters, 1974; Canadian Hydrographic Service,

1970). The lake’s elevation above sea level is

74.01 meters (242.8 feet) (Canadian

Hydrographic Service, 1970) which makes its

depth of crypto-depression to be 170 meters

(558 feet). Lake Ontario’s length is 307 km (191

mi), with a maximum width of 87 km (54 mi)

(Canadian Hydrographic Service, 1970). The

lake’s shoreline length is 1,380 km (858 mi)

(Beeton and Chandler, 1963). The variation of

volume and area with depth is illustrated in

Appendix A (Figures A-1 and A-2).

4.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 (International Joint

Commission, 1969), 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.28 cm (32.79 in) and the average

annual evaporation is 71 cm (28 in) (Great Lakes

Basin Commission, 1971).

Table 1. Flow Budget

Source Flow (cfs) % of Discharge

Major Tributaries

 Niagara River

 Twelve Mile Creek

 Oswego River

 Trent River Region

 Black River

 Genesse River

 Welland Canal

Other Tributaries

Other Sources

 Precipitation -

 Evaporation and

 Other)

Discharge

 St. Lawrence

 River

 195,000

 6,400

 6,200

 4,200

 3,800

 2,700

 1,200

 5,900

 6,600

 232,000

 84.1

 2.8

 2.7

 1.8

 1.6

 1.2

 0.5

 2.5

 2.8
7

The thermal bar development which precedes

the development of the full thermocline begins in

late April or early May in Lake Ontario

(International Joint Commission, 1969). 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 (International Joint

Commission, 1969). 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

Chapter 6).

4.3 NUTRIENT INPUTS

Nitrogen and phosphorus are the nutrients used

in the present configuration 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 (IJC,

1969). Table 2 indicates 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 average inflow of 232,000 cfs. The

inorganic form of the nutrients were set at the

lake’s winter concentrations with the remaining

load inputted in the non-living organic form.

4.4 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) (CCIW, 1969).

2. STORET, U.S. Environmental Protection

Agency (USEPA).

3. Report to the International Joint

Commission on the Pollution of Lake

Ontario and the International Section of the

St. Lawrence River (IJC, 1969).

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

adequate dense spatial grids (see Figure 2 for a

typical grid) but also comprised good temporal

coverage for the years surveyed (Table 4).

STORET, the USEPA’s water quality storage

and retrieval system is the prime residence of all

United States collected water quality data. The

pre-IFYGL STORET database consisted mainly

of data from three 1965 synoptic cruises. The

temporal coverage by these cruises was not

sufficient to be of direct use for verification

analysis. The data were, therefore, considered

only as supplemental to the basic data set.

STORET’s utility is for storage and analysis of

IFYGL data and as a database for tributary

nutrient concentration data. It was, therefore,

decided to use the data contained in the LDR as

a verification database for the preliminary

models and the STORET data as a verification

database for the 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 expanded

on this task, with the principal product being
8

Table 2. Nutrient Loadings

Source Nitrogen Phosphorus

lbs/day kg/day % lbs/day kg/day %

Niagara

Tributaries

Municipal

Industrial

 522,200

 191,000

 72,600

 97,300

 236,900

 86,600

 32,900

 44,100

 59

 22

 8

 11

 42,200

 15,600

 16,200

 1,000

 19,100

 7,100

 7,300

 500

 56

 21

 22

 1

Total 883,100 400,500 75,000 34,000

Table 3. Nutrient Load Distribution

System Nutrient Load As Boundary Concentration

lbs/day kg/day mg/l

Nitrogen

 Non-Living Organic N

 Ammonia N

 Nitrate N

Phosphorus

 Non-Living Organic P

 Inorganic P

 551,400

 18,800

 312,900

 57,100

 17,900

 250,100

 8,500

 141,900

 25,900

 8,100

 0.443

 0.015

 0.250

 0.457

 0.0143
9

Figure 2. Typical CCIW cruise tracks, April 1969.
10

Table 4. CCIW Limnological Data Report Summary (Number of Cruises With Adequate Spatial

Definition)

1966

(June-Sept)

1967

(June-Nov)

1968

(April-Nov)

1969

(April-Dec)

Nitrogen:

 Organic or Kjeldahl-N

 Ammonia Nitrogen

 Nitrate Nitrogen

Phosphorus:

 Total Phosphate

 Reactive Phosphate

Temperature

Chlorophyll a

-

-

4

-

5

9

-

11

10

10

4

10

11

11

-

5

6

5

6

9

6

-

9

9

8

9

9

9

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

reduction.

The Lake 1 model was designed to aid in the

understanding of the principal mechanism which

underlie the process of eutrophication. A key

criterion in evaluating the validity of a model is to

compare the model’s 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 fo

the evaluation of model output (see Chapter 6).

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

Figure 3. Typical output from data reduction - main lake average at 1 meter.
12

CHAPTER 5

LAKE ONTARIO - MODEL STRUCTURE

5.1 THEORETICAL BACKGROUND

The basic theory used for the Lake Ontario

model makes use of previous model structures

(Di Toro et al., 1971; Hydroscience, Inc., 1973;

Thomann et al., 1974) 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, and

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:

 (1)

where sk is the kth dependent variable (biological

or chemical), u, v, w are the velocity components

in the x, y and z directions, respectively, Ex. Ey

and Ez are the dispersion coefficients in each

respective spatial direction, Sk represents the

kinetic interactions between the k variables and

W is the direct inputs of the substance, sk. It

should be noted that Equation (1) as written,

does not include any interaction between the

variables sk and 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 Equation (1) is:

 (2)

where [E] is a diagonal matrix of dispensing

coefficients, (sk) is a column vector, U is the

velocity vector and 
13

In Equation (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 Equation (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 Chapter 8. Primary interest, therefore,

centers about the reaction kinetics, their function

forms and numerical values.

A finite difference approximation to Equation (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 Vj

, then it can

be shown that 

 (3)

where [V] is an n x n diagonal matrix of volumes,

(s)k is an n x 1 vector of the variables sk, [A] is an

n x n matrix of advection and dispersion terms;

(S)k is an n x 1 vector of kinetic interactions and

(W)k is an n x 1 vector of inputs of variable sk.

For m interacting variables and n physical

locations, Equation (3) indicates that a total of m

x n equations must be solved. In general, the

equations are both linear and nonlinear.

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

preparation necessary to obtain solutions are of

a substantial magnitude. 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 Figure 4.

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

chloride 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

transport 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

ten 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 models 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 a

representation of the carbon cycle.
14

Figure 4. Eutrophication modeling strategy.
15

Figure 5. System diagram - Lake 1 model.
16

5.2 PHYTOPLANKTON CHLOROPHYLL

The basic kinetic interactions, for phytoplankton

chlorophyll (a measure of biomass), is given for

a single volume, j, as

 (4)

where P is the phytoplankton chlorophyll (μg/l),

Gp and Dp are the growth and death rate (1/day),

respectively, and wji and wij are the sinking

velocities of the phytoplankton between

segments j and i.

5.2.1 Phytoplankton Growth Rate

The growth rate expression is similar to that

developed previously (Di Toro et al., 1971) and

as used in this work is given by

 

 = (temperature effect)  (light effect) 

 (nutrient effect) (5)

where Gmax is the maximum growth rate, T is

water temperature, Is is light intensity at Gmax, Ia

is incoming solar radiation, ke is the extinction

coefficient, H is the depth of the segment, kmn is

the half-saturation constant for total inorganic

nitrogen, N’, N2 is ammonia nitrogen, N3 is nitrate

nitrogen, Kmp is the half-saturation constant for

phosphorus and p2 is the available phosphorus.

The temperature is given by Eppley (1972) as:

 (6)

where Eppley reports Gmax at 0.587 (1/day). This

form permits phytoplankton growth at low

temperatures and at 20°C, results in a growth

rate of 2.1/day.

The function form for the light effect is as given

in previous work (Di Toro et al., 1971) and

incorporates vertical extinction of solar radiation

and self-shading effects.

The form is

 (7)

where

and k’

e is the light extinction coefficient at zero

chlorophyll 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

 (8)

The range of values for kmn and kmp, the halfsaturation constants for nitrogen and

phosphorus, respectively, has been documented

previously (Di Toro et al., 1971). More recent

work (Caperon and Meyer, 1972; Hendrey and

Welch, 1973) tends to support the original range;

namely kmn in the range 5-50 μg N/l and kmp in

the range 1-10 μg/l. There is an obvious species

dependence of this constant as well as a

possible dependence on the eutrophic state of

the lake. In the Lake 1 model, kmn is used at 25

μg N/l and kmp of 1 and 10 μg P/l is used as a

range.
17

At the present time, silica limitation is not used in

this model. Goering et al. (1973) and Paasche

(1973) have estimated the half-saturation value

for silica at about .02 - .08 mg Si/l for marine

species. Under peak growing conditions in Lake

Ontario, surface silica values are about 0.1-0.2

mg Si/l 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, Gpj is then

given by the product of Equations (6), (7), and

(8). Although there is no a priori reason for using

a product formulation, various experiments (see

Equation (1)) have tended to support the

approach. Another expression has been

proposed on the grounds that Equation (5) is too

severe in reducing growth due to the product

term. The alternate construct is Equation (9):

 (9)

or in general

 (10)

where k is a normalization factor and μi

represent the general reduction factors. In

general, from Equation (5) to Equation (10):

the extent of the inequality depending on μi

.

Undoubtedly then, Equation (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. Equation

(10), therefore, appears to be less severe than

Equation (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 (1934) (cited in Di Toro et al.,

1971) and more recently by Di Toro (1974a) for

the flask experiments. The Lake 1 model was

modified to test the sensitivity of the solutions to

Equation (5) or Equation (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 cross growth rate of

the phytoplankton biomass, Gpj Pj

 is equivalent to

the daily average productivity in the jth segment.

Thus,

 (11)

where vj

 is the average rate of carbon fixation

(mg C m2

-day) and acc is the carbon to

chlorophyll ratio (range 20-100 mg C/mg Chlor).

5.2.2 Phytoplankton Death Rate

The expression for the death rate of the

phytoplankton includes two primary effects:

endogenous respiration and predation by

herbivorous zooplankton. The death rate is

given by

 (12)

where K2 is a constant (1/day), Cg is the grazing

rate of the herbivorous 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 mg

C-day-°C appears appropriate and at 20°C,

represents about 2-6 ml filtered per individual per

day.
18

5.2.3 Phytoplankton Sinking

In Equation (4), the sinking of phytoplankton

from one segment to another segment can be a

potentially significant sink or source of

phytoplankton 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

(Hutchinson, 1967) and apparently the settling

velocity of nutrient rich cells is somewhat less

than cells that are nutritionally deficient (Yentsch,

1962). 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 (1974) have investigated the

settling velocity profile of particulate organic

carbon in Lake Ontario with an ingenious in situ

measuring device. Their results show a general

decrease in settling velocity with depth and a

marked seasonal variation of settling velocity.

For example, on July 19, 1972, the settling

velocity in the thermocline (about 10 m) was

estimated at about -0.3 m/day, while at 50

meters depth, the settling velocity was about

+1.3 m/day. But in March 1973, the velocity

profile was all positive ranging from 0.9 m/day at

20 meters to 0.1 m/day at 140 meters. Overall,

the values given by Burns and Pashley (1974)

varied between -0.4 and +2.0 m/day. These

results and others tend to indicate that the

settling velocity should be related in some ways

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.05 m/day to 0.5 m/day. The

sensitivity of the results to the phytoplankton

settling velocity is discussed below.

5.2.4 Zooplankton Carbon

The measure of the herbivorous zooplankton

populations z(1) is taken as the carbon content of

the biomass and all sources of food for z(1) are

the phytoplankton chlorophyll. The general

expression for a single volume j is therefore:

 (13)

where G(1)

z and D(1)z are the growth and death

rates respectively of z(1)

z. The form of G(1)z is

 (14)

where azp is the zooplankton carbon to

phytoplankton chlorophyll efficiency, Cgzp is the

grazing rate of the herbivorous zooplankton and

Kg is the half-saturation concentration of

phytoplankton at which the growth of z(1) 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,

 (15)

where K(1)

z is the endogenous respiration rate as

a function of temperature (1/day-°C), Cg21 is the

grazing rate (liters/day mg C-°C) and z(2) 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 Equations (13) and (14).

Thus:
19

 (16)

where 

 (17)

and

These expressions are an obvious

oversimplification of a complex food web and

reflect only a “linear” type of food chain (see

Figure 5). Indeed the effect of the zooplankton

on the phytoplankton as formulated above is to

always decrease the population, yet it has been

shown (Porter, 1973) that zooplankton can also

increase phytoplankton population or have no

effect. The mechanism for an increase in the

population is apparently through the passage of

gelatinous green algae through the zooplankton

gut with subsequent breakup 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

characteristics of Lake Ontario.

The zooplankton of Lake Ontario are dominated

by the crustaceans specifically of the order

Copepoda and Cladocera. Most hauls are made

over 0-50 meter depth range. Principal copepod

species given by Patalas (1969) are Cyclops

bicuspidatus and Tropocyclops prasinus and of

the cladocerans: Daphnia retrocurva and

Bosmina longirostris. In Patalas’ work in 1969,

these four species accounted for 91% of the total

zooplankton and averaged about 60/liter from

June-October.

Others (Glooschenko et al., 1972; McNaught and

Buzzard, 1973; Watson and Carpenter, 1974)

generally showed values in similar range, and in

some cases, up to a maximum of almost

200/liter. For use in biomass computations,

conversions must be made to equivalent carbon.

Such conversion depends on age and

accompanying weight of the animal and can,

therefore, vary over a wide range. Watson and

Carpenter (1974) in their conversions used

weights of species to arrive at total biomass and

equivalently averaged about 2.8 μg dry

weight/ind. The results using an average of 40%

carbon by dry weight show peak values generally

in August at lakewide averages of between 0.07-

0.12 mg carbon/liter.

At the average dry weight of about 1-5 μg/ind

grazing rates of about 0.8-2 ml/ind-day are

appropriate (Di Toro et al., 1971; Kibby, 1971).

Indeed the rate for Bosmina longirostris, a

dominant species in Lake Ontario has been

given as 1-3 ml/ind-day (Hutchinson, 1967).

Grazing rates expressed in liters per mg

zooplankton carbon per day have been assumed

to be linearly related to temperature and have

been used in the range of .03-.06 (1/mg C-day-

°C) corresponding to the reported grazing ranges

at 20°C.

5.3 NITROGEN SYSTEM

The basic equations for the nitrogen sub-system

follow previous work (Hydroscience, Inc., 1973;

Thomann et al., 1974) and include non-living

organic nitrogen, N1, ammonia nitrogen, N2 and

nitrate nitrogen, N3. The equation for the sourcesink term for organic nitrogen is:

 (18)

where a1 is nitrogen-to-chlorophyll ratio, and ac is

the carbon-to-chlorophyll ratio and K1 is the

overall decay of N1.

The first and third terms represent organic

nitrogen released through endogenous

respiration by the zooplankton and
20

phytoplankton, respectively; and the second term

represents the organic nitrogen of the grazed but

unassimilated phytoplankton. The last term

represents the decomposition, settling and other

effects that contribute to the decay of organic

nitrogen.

For ammonia,

 (19)

where K12 is the rate of production of ammonia

from organic nitrogen, K22 is the rate of oxidation

of ammonia to nitrate, and  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,

 (21)

where K23 = K22, the production of nitrate by

ammonia oxidation.

5.4 PHOSPHORUS SYSTEM

The equations for the phosphorus cycle used in

the models are as previously developed

(Hydroscience, Inc., 1973; Thomann et al., 1974)

and include non-living organic phosphorus and

phosphorus available for phytoplankton growth.

For organic phosphorus, the equation for the

source and sink for p1 is:

 (22)

where ap is the phosphorus-to-chlorophyll ratio

and K1 represents a general decay of organic

forms.

For the inorganic phosphorus, assumed as

phosphorus available to the phytoplankton, the

source-sink equation for p2 is:

 (23)

where K12 represents the rate of decomposition

of organic phosphorus to forms available for

phytoplankton utilization and K22 represents an

overall loss of inorganic phosphorus.

5.5 CHEMICAL MODELS - DISSOLVED

CARBON DIOXIDE

The inclusion of chemical reactions into models

of Great Lakes water quality, by contrast to the

biochemical and biological reactions considered

above, poses a new set of problems by virtue of

the rapid rates and the highly nonlinear forms of

the governing 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 presentation of these

results is available (Stumm and Morgan, 1970).

Application of certain models to the Great Lakes

water chemistry has been undertaken by Kramer

(1964, 1967) 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
21

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 alumino-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 as a prototype for all other

such reactions. In fact, the procedure is the

same for the inclusion of all such reactions within

the context of conservation of mass equations.

5.5.1 Equilibrium Analysis

The carbonate equilibria is treated in the

conventional way (Stumm and Morgan, 1970;

Trussell and Thomas, 1971). The major species

considered are aqueous or dissolved [CO2] (the

concentration of hydrated CO2 being small,

bicarbonate [HCO3

-

], and carbonate [CO3

=

] ions,

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 complicated by the

fact that they undergo reversible ionization

reactions, e.g., CO2 + H2O = H+

 + HCO3

-

 which

occur very rapidly relative to the time scales of

mixing and gas transfer (Kern, 1960). In fact, it

is clear that the individual species are each

extremely reactive so that a direct mass balance

formation results in equations which are

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

relative to the ionization reactions. A formal

mathematical discussion is given subsequently.

However, it is easy to see the principle involved.

Consider the total inorganic carbon

concentrations.

 (24)

It is clear that any reversible ionization reaction

will not affect the CT of the system (we assume

no carbonate precipitation are forming or

dissolving). Thus, a mass balance equation for

CT can ignore the 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 ligands, e.g., Na

+

, K+, Ca

++, ..., and

Cl-, SO4

=

, ..., which are present and assumed not

to react appreciably with the carbonate or

hydrogen-hydroxide species. Then, it is clear

that a charge balance equation:

 

 (25)

is unaffected by ionization reactions so that the

quantity called alkalinity

 (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

quantities for this system, the other common

conservative quantities (Stumm and Morgan,

1970) are linear combinations of CT and A.

If CT and Alk are known, the concentration of all

the species being considered is easily obtained.

In fact, using the equilibrium equations (Stumm

and Morgan, 1970; Trussell and Thomas, 1971)
22

 (27)

it is true that (Trussell and Thomas, 1971)

 (28)

where

 (29)

and

 (30)

It does not seem to have been noticed that if the

alkalinity concentration is large relative to (H)

= [OH-] - [H+], which is the common situation

then, Equation (28) becomes

 (31)

so that only the ratio of alkalinity to total

inorganic carbon is important in determining the

resulting pH whenever Alk = Alk + [H+

] - [OH-].

Using tabulated values of K1 and K2 for pore

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 critical conservation 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.

The equations which govern the distribution of

Alk and CT are available as mass balances,

independent of the ionization reactions. For

example, in a one-dimensional vertical model the

equations would be

 (32)

 (33)

with a zero mass flux boundary condition at z =

0 for alkalinity and an air-water exchange

condition for CT:

 (34)

representing the exchange of carbon dioxide.

The difficulty in solving these equations is due to

this nonlinear boundary condition since CO2 (aq)

is a nonlinear function of the dependent

variables of the problem: [Alk] and [CT], 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 reaction (if bacterial

synthesis is ignored) can be represented as:

 (35)

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

Similarly as phytoplankton grow they utilize

inorganic carbon and this sink would be included

in SCT. These mass balance equations would

then be solved in conjunction with the equations

for the other species of interest.
23

Figure 6. Relationship of pH and the ratio, Alk/CT. 

 
24

5.5.2 Linear Approximation

It is clear that mass balance equations which

ignore the ionization reactions are only correct if

the dependent variable is a conservative quantity

relative to these reactions. Consider the

conservative quantity [CT] - [Alk] = [CO2 - Acidity].

Clearly it is conservative since it is the difference

of conservative quantities. From Equations (24),

(25) and (26), it is the case that

 (36)

But this equation implies that over a certain

range of pH and CT, the major component of

CO2-acidity is [CO2] itself. Figure 7 presents the

% error in the approximation [CO2-Acy] = [CO2].

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 CO2 does not ionize to HCO3

-

 and H+

but rather remains as CO2. Hence, ionization

reactions can be ignored. Thus, from Equations

(32) and (33) the mass balance equation for

CO2-acidity is

 (37)

and the reaeration source of [CO2-Acy] can be

approximated as KL ([CO2]sat - [CO2-Acy]) so that

within the errors in Figure 7, the equation for

[CO2-acidity] is linear. 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 contrast to the

nonlinear CT equation which must be solved

with all source 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 simplification 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.

5.5.3 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 implementation of this

dependency poses no serious problem since

standard curves fitting techniques can be applied

to the experimental data. However, such an

approach backs 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 reaction at standard

state, Gr

° via the equation:

 (38)

where

 P = universal gas constant

 T = temperature °K
25

Figure 7. Percent error in assuming [CO2 - Acy] [CO2].
26

and further that the free energy change is related

to the change in enthalpy, Gr

°

 and entrophy,

Sr

°

 as follows:

 (39)

Hence, the equilibrium constant is given by:

 (40)

This equation is strictly true at any temperature;

the difficulty is that Hr

°

 and Sr

°

 are also

functions of temperature. The conventional

temperature at which they are available is

298.15°K or 25°C. Table 5 presents the relevant

thermodynamic properties for the carbonate

system and Table 6 gives the properties for the

reactions of interest.

Table 5. Thermodynamic Properties at 25°C

Compound

HGf

°

(kcal/-

mole)

HGf

°

(kcal/-

mole)

HGf

°

(kcal/-

mole)

CO2(g)

H2CO3(aq)*

HCO3

CO3

=

H+

H2O

OH-

 -94.051

-167.22

-165.39

-161.84

 0

 -68.315

-54.970

 -94.254

-148.94

-140.26

-126.17

 0

 -56.687

 -37.594

 51.06

 44.8

 21.8

-13.6

 0

 16.71

 -2.57

Source: Wagman et al. (1968).

*H2CO3(aq) = H2CO3 + CO2(aq);

 [H2CO3(aq)] = [CO2/(aq)].

Table 6. Thermodynamic Properties at 25°C

Reaction

Hr

°

(kcal/-

mole)

Sr

°

(kcal/-

mole)

1. H2O = H+

 + OH2. CO2(g) + H2O = H2CO3

3. H2CO3 = H+ + HCO3

-

4. HCO3 = H+ + CO3

=

13.345

 -4.854

 1.830

 3.55 

-19.28

-22.97

-23.0

-35.4

Two approximations are investigated:

1. Assume that Hr

°

 and Sr

°

 are constants

over the temperature range of interest

(Stumm and Morgan, 1970).

2. Apply a correction which attempts to

evaluate the change in the specific heat as

a function of temperature (Helgeson, 1967).

The equation for the second approximation is:

 (41)

where

 (42)

and Tr = 298.15 °K

  = 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 nonlinear relationship giving

somewhat better results at the lower

temperatures. Thus, to within a few percent,

these approximations are quite adequate.
27

Table 7. Temperature Dependence Comparison

Reaction = HCO3

-

 = H+ + CO3

=

Enthalpy of Reaction = 3550.0000 cal/mole

Entropy of Reaction = 55.4000 cal/deg-mole

Temp °C Experimental Pk Method 1 Pk* Percent Error Method 2 Pk** Percent Error

0.0

5.0

10.0

15.0

20.0

25.0

30.0

10.6250

10.5570

10.4900

10.4300

10.3770

10.3290

10.2900

10.5770

10.5260

10.4767

10.4292

10.3832

10.3389

10.2959

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

Reaction = H2CO3

*

 = HCO3

-

 + H+

Enthalpy of Reaction = 1830.0000 cal/mole

Entropy of Reaction = -23.0000 cal/deg-mole

Temp °C Experimental Pk Method 1 Pk* Percent Error Method 2 Pk** Percent Error

0.0

5.0

10.0

15.0

20.0

25.0

30.0

6.5790

6.5170

6.4640

6.4190

6.3810

6.3520

6.3270

6.4908

6.4645

6.4391

6.4146

6.3909

6.3681

6.3459

-1.34

-.81

-.38

-.07

.16

.25

.30

6.5175

6.4816

6.4489

6.4191

6.3922

6.3681

6.3466

-.93

-.54

-.25

.00

.18

.25

.31

Reaction = CO2(g) + H2O = H2CO3

*

Enthalpy of Reaction = 4854.0000 cal/mole

Entropy of Reaction = 22.9700 cal/deg-mole

Temp °C Experimental Pk Method 1 Pk* Percent Error Method 2 Pk** Percent Error

0.0

5.0

10.0

15.0

20.0

25.0

30.0

1.1100

1.1900

1.2700

1.3200

1.4100

1.4700

1.5300

1.1364

1.2062

1.2735

1.3385

1.4013

1.4620

1.5207

2.37

1.36

.28

1.40

-.62

-.54

-.61

1.1630

1.2233

1.2832

1.3430

1.4026

1.4620

1.5214

4.78

2.80

1.04

1.74

-.53

-.54

-.56
28

Table 7. Temperature Dependence Comparison (Continued)

Reaction = H2O = H+

 + OHEnthalpy of Reaction = 13345.0000 cal/mole

Entropy of Reaction = -19.2800 cal/deg-mole

Temp °C Experimental Pk Method 1 Pk* Percent Error Method 2 Pk** Percent Error

0.0

0.5

10.0

15.0

20.0

25.0

30.0

14.9435

14.7338

14.5346

14.3463

14.1669

13.9965

13.8330

14.8911

14.6991

14.6140

14.3352

14.1626

13.9958

13.8344

-.35

-.24

-.14

-.08

-.03

-.01

.01

14.9134

14.7135

14.5221

14.3390

14.1636

13.9958

13.8350

-.20

-.14

-.09

-.05

-.02

-.01

.01

 *Stumm and Morgan, 1970

**Helgeson, 1967

5.5.4 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 (Shapley and Cutler,

1970). 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 CO2 - H2O

for which the partial pressure of CO2(g) is fixed

and varying amounts of acidity and alkalinity are

added. The results are computed for the

temperature range 0° to 30°C using the linear

approximation for the temperature dependence

of Gr

°

 for the three reactions of interest. Figure

8a gives the total inorganic carbon that is in

equilibrium with an atmosphere for which the

partial pressure of CO2(g) is 10-3.5

atmospheres. As alkalinity is added, the CT

increases and it is approximately true that [CT] =

[Alk] 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 Figure 8b.

For varying partial pressures of CO2(g) and Alk

= 0, the results are shown in Figure 8c and 8d.

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.

5.5.5 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:
29

Figure 8. Computations for CO2 - H2O system.
30

 (43)

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

sides of Equations (32) and (33) are denoted by

D[Alk] and D[CT]. The conservation of mass

equations for each species is obtained by

equating mass transport to sources and sinks:

 (44)

 

 (45)

 (46)

 (47)

 (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 kv’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 kv’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 variables which are

appropriate sums of the species in such a way

that the reaction terms cancel out. Thus if CT =

[CO2] + [HCO3

-

] + [CO3

=

] then by adding

Equations (44), (46), and (47) the result is:

 (49)

The other variable is, of course,

 

 (50)

and since D is a linear operator,

 (51)

and using Equations (45), (46), (47), and (48) the

result is Equation (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 CT are exact, regardless of

the magnitude of the kv’s. They are the invariant

of the reactions (Shapiro, 1962) 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,

Equation (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 reaction of

the form.
31

was not explicitly considered. It is a simple

matter to verify that the inclusion of this reaction

does not affect the validity of the resulting

equations for Alk and CT but no guarantee has

been presented that such a reaction does not

exist. This can be done (Shapiro, 1962)

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 CT, 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 invariant quantities are always

easy to find has been structured (Di Toro,

1974b). The results are essentially a

generalization of the previous derivations.
32

CHAPTER 6

LAKE 1 MODEL

6.1 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 Figure 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 nearshore versus 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

indicate substantial horizontal gradients,

although variations do exist in certain areas.

(Lake Ontario can be contrasted in this regard to

Lake Huron which 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. Figure 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 (c) vertical

settling of the phytoplankton. The lake shown in

Figure 9 has a defined epilimnion 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 the summer and then

goes through a fall overturn. 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 Figure

9 and Table 8 and using the general Equation

(3), the equations for Lake 1 are:

Segment Number 1 - Epilimnion:

 (52)

Segment Number 2 - Hypolimnion:

 (53)

Segment Number 3 - Sediment:

 (54)
33

Figure 9, Major physical features included in Lake 1 model.
34

Table 8. Basic Physical Data of the Lake 1 Model

Segment

Number

Segment

Interface

Volume

(m3 x 106) %

Depth

(m)

Surface

Area (m2) cfs

Flow

m3/sec %

1

2

3 (sediment)

1-2

2-3

 297,000

1,373,000

 -

19

81

17

73.3

0.15*

1.6 x 1010

0.89 x 1010

 43,500

188,500

1232

5323

19

81

Note: Vertical dispersion coefficient between segments number 1 and 2 varied from 0-90 cm2

/sec 

 (0-778 m2/day).

*Segment number 3 depth is arbitrary.

where sk,B1 and sk,B2 are the concentrations of the

kth system inputted into the epilimnion and

hypolimnion, respectively, QB1 and QB2 are the

flows into segments 1 and 2, V1 is the volume of

segment i, E12 is the vertical turbulent exchange,

A12 is the cross-sectional area between the

epilimnion and hypolimnion and z is the

average depth of the two segments. As shown

in Figure 5, the Lake 1 model consists of ten

systems (k = 1...10). The kinetics are as given

previously in Chapter 5. 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.

6.2 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 sensitivity 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 the IFYGL.

These comparisons are qualitative in the sense

that an assessment of the degree to which the

computed output 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 system diagram

(Figure 5) and secondary variables constructed

for the behavior of the model. The variables

examined, therefore, for calibration and

verification are:

1. Phytoplankton chlorophyll “a” - g/l

2. Primary productivity - mg carbon/m2

-day

3. Total zooplankton carbon - mg C/l

4. Bottom deposition rate - mg C/m2

-year

5. Total phosphorus - mg/l

6. Available phosphorus (orthophosphate) -

mg/l

7. Total Kjeldahl nitrogen - mg/l

8. Ammonia nitrogen - mg/l

9. Nitrate nitrogen - mg/l

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 rates) the results are

generally always consistent with the observed

data. The results of the comparison to observed
35

data are presented after the results of the

preliminary sensitivity and calibration runs.

6.2.1 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 seven 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: (1) test

program elements, (2) study the behavior of the

system and its sensitivity to various system

parameters and inputs, and (3) prepare a

preliminary verification of the Lake 1 model using

data collected prior to IFYGL. These early runs,

therefore, 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 example, variable levels of spring

and fall vertical mixing, and settling velocities for

phytoplankton in addition to zooplankton grazing

rates, and effect of half-saturation constant for

phosphorus. Some of the output from these

areas is discussed below.

6.2.2 Variable Vertical Mixing

Vertical mixing plays an important role in the

dynamics of phytoplankton 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 et

al. (1974) estimated values from 0.1 cm2/sec to

22 cm2

/sec depending on effects such as wind

conditions and density stratification. In the Lake

1 model, 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 90 cm2/sec. More

extensive analysis 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 this chapter.

Figure 10 shows the effect of vertical dispersion

on chlorophyll a. The runs include a sinking

velocity of 0.05 m/day. A bimodel distribution of

phytoplankton occurs in all cases due to the

interaction of the two higher zooplankton levels

used in these runs. As populations build up in

the spring, the herbivorous 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 Figure 10, the main effect of vertical

mixing is in the spring, where populations

increase more rapidly when no mixing is allowed

and reaches 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.

6.2.3 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

phytoplankton sinking in quiescent waters (see

Chapter 3) 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
36

Figure 10. Effect of vertical mixing.
37

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 grazing was set at 0.06

1/mg C day °C, no vertical mixing was used and

the Michaelis constants were 25 μg/l and 10 μg/l

for nitrogen and phosphorus, respectively. At a

velocity of 0.5 m/day phytoplankton populations

never exceed about 3 μg chlor/l and total

zooplankton carbon (Figure 8) never exceeds

about 0.08 mg C/l. 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.05 m/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, Kmp,

has an important effect, however. An order of

magnitude reduction in this coefficient to 1 μg/l

does permit growth of the phytoplankton at the

higher settling velocity of 0.5 m/day. Figure 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 Figure 13, which is a plot of nitrate and

available phosphorus concentrations calculated

under the conditions of Kmp = 10 μg/l. In

addition, Figure 13 shows the 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 Equation 8). It can be

seen that nutrient uptake and growth limitation is

minimal for the case of sinking velocity = 0.5

m/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 Figure 13,

nitrogen does not significantly limit growth for the

case of Kmp = 10 μg/l.

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 0.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 calculated 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 Figure 11). Two

effects are occurring: (1) herbivorous

zooplankton are growing rapidly and (2)

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

6.2.4 Effect of Michaelis Constant

As indicated in Equation (5), the growth of the

phytoplankton is considered to be limited by a

Michaelis-Menten product term of total inorganic

nitrogen and available phosphorus given by

Equation (8). The dynamics of growth, therefore,
38

Figure 11. Effect of phytoplankton settling velocity, Kmp = 10 μg/l.
39

Figure 12. Effect of phytoplankton settling velocity, Kmp = 1 μg/l.
40

Figure 13. Effect of phytoplankton settling velocity on nutrients and nutrient limitation, Kmp = 10 μg/l.
41

are affected by the levels of Kmn and Kmp, the

half-saturation constants for nitrogen and

phosphorus, respectively. Several sensitivity

runs have been made, therefore, specifically

examining two levels of Kmp, 1 and 10 μg/l both

of which encompassed the range of published

values. Phosphorus is emphasized because of

the waste removal programs presently underway

in the Great Lakes basin.

Figures 14 and 15 show the results of two

computations at an order of magnitude

difference in Kmp. For both runs, settling velocity

was 0.05 m/day, Kmn (nitrogen half-saturation

constant) is 25 μg/l and zooplankton are

included as before. As seen in Figure 14 at the

1 μg/l level, phytoplankton biomass grows earlier

more rapidly and reaches a peak value of 14 μg/l

chlorophyll/l at the end of April. At the 10 μg/l

value, growth is later, slower and reaches a

maximum value of about 9 μg/l almost a month

later at the end of May. Values in the late

summer are comparable. The nutrient plots

shown in Figure 14 indicate that in May,

phosphorus is limiting under the 1 μg/l case

while in August-September, nitrogen appears to

be limiting. This is shown more clearly in Figure

15 - which is a plot of the nitrogen and

phosphorus limitation terms under both

conditions. It can be noted that at 1 μg/l, 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 hand, at a

phosphorus Michaelis constant of 10 μg/l, 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 halfsaturation 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 program.

6.2.5 Effect of Alternate Growth Formulation

As discussed in Chapter 3, Equations (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 Equation (9). The results for three

output variables are shown in Figure 16. Both

runs displayed therein are for identical

conditions, except for the alternate growth

formulation as indicated in the figure.

Equation (9) results in immediate and rapid

growth of the phytoplankton with subsequent

total depletion of the phosphorus. The dynamic

behavior of the results using the formulation of

Equation (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 consistent with the observed data (see

Figure 3). From this one comparison, then,

Equation (5) the original product formulation is

still to be preferred.

6.2.6 Time to Dynamic Equilibrium

With the relatively long detention time of about

eight 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

 (55)

where so is the initial condition of s and to is the

detention time given by V/Q. For K = 0, the

conservative case, the time to reach 95% of a

new steady value is about three detention times,

which for the whole of Lake Ontario is about 24
42

Figure 14. Effect of Michaelis constant for phosphorus on primary variables.
43

Figure 15. Effect of Michaelis constant for phosphorus on nutrient limitation of nitrogen and

phosphorus.
44

Figure 16. Responses due to two phytoplankton growth formulations.
45

years. Furthermore, 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 inputs,

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

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

μg/l. As shown, peak phytoplankton levels in the

spring gradually decrease from greater than 9

μg/l in the first year to about 7 μg/l 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 Figure 18 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 3 μg/l, the difference

is about half of that at equilibrium.

The lower two plots of Figure 18 are quite

interesting and illustrate further the difference in

the dynamic behavior under the two phosphorus

Michaelis levels. Under the lower level of 1 μg/l,

equilibrium is reached almost immediately in

both phosphorus and total inorganic nitrogen.

Under an order of magnitude increase in the

level to 10 μg/l, dynamic steady-state is not

reached for the inorganic nitrogen until almost

two detention times or 16 years. With the lower

level of 1 μg P/l, 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 of

Figure 18, the system becomes nitrogen limited.

The times of the year of the nitrogen and

phosphorus limitation are different, however, as

shown in Figure 15.

In essence then, these runs illustrate the need to

carry out any simulations of a proposed nutrient

reduction program about 8-16 years to determine

the expected phytoplankton levels under a

different waste load input. If, however, the loads

are changing from year-to-year, as for example,

in a linear increasing or decreasing fashion, then

even longer times may be required. Indeed,

under changing loads, a dynamic equilibrium

may never be theoretically possible. For many

practiced purposes, the system will probably

reach steady-state in a 10-20 year time horizon

after significant changes have been made.

6.2.7 Preliminary Comparison

Figures 19 and 20 show one of the early

comparisons of observed data and model output.

The Michaelis constant for phosphorus was 10

μg/l and settling velocity was 0.05 m/day. In the

preliminary comparisons, a variety of similar

plots was generated. Figures 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 Figures 19 and 20, however,

remained to be exploited in greater depth before

a “final” verification could be obtained. For
46

Figure 17. Sixteen-year computed output of phytoplankton chlorophyll.
47

Figure 18. Long-term behavior of model components under two phosphorus Michaelis constants.
48

Figure 19. Preliminary comparison of model output and observed data.
49

Figure 20. Preliminary comparison of model output and observed data (continued).
50

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

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

good a comparison as possible to available data.

6.3 RESULTS OF “FINAL” VERIFICATION

ANALYSIS

Data were available for the years 1967-70

collected by the Canada Centre for Inland

Waters (CCIW) 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, a

consistent set of parameter values which

satisfactorily explains the observed data on a

variety of different variables has emerged. A

plausible explanation for the dynamic behavior of

the phytoplankton is, therefore, possible and is

discussed below.

Figures 21 and 22 show the verification of the

Lake 1 model upper layer (0-17 meters) while

Figure 23 shows the verification for the

hypolimnion (data in range, 50-150 meters).

Table 9 shows the principal coefficients used for

the runs. As shown, the verification for plankton

in the epilimnion is quite good and satisfactorily

duplicates the spring peak, subsequent midsummer die-off and fall bloom. Five other

variables in addition to the phytoplankton

chlorophyll are satisfactorily verified in the

comparison. It should be noted here that it is

generally not possible to obtain the verification

shown in Figures 21-23 by arbitrary specification

of the coefficients, so that the results shown do

not simply represent a “curve fitting” exercise.

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

computation of the components of the dynamic

behavior, considerable insight can be gained by

examining the influence of various phenomena

on the growth and death rates of the

phytoplankton.

Figure 24 shows the dynamics of the kinetic

growth rate of the phytoplankton (Gp in Equation

(4)). 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 substantial 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, Gp is lowest at the

end of May and peaks at about 0.25/day in

August. The peak value represents an overall

reduction of 90% from saturated growth

conditions. Primary productivity values

calculated from the growth rate and

phytoplankton biomass, Equation (11), give peak

values of about 600 mg C/m2-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 Glooschenko et al.

(1974) while the latter is somewhat lower than

the Glooschenko’s estimated spring values of

500-1100 mg C/m2

-day.

Figure 25 shows the dynamics of the kinetic

death rate of the phytoplankton (Dp in Equation

(4)). Three effects are to be noted: (1)

phytoplankton settling, (2) water temperature

effect on endogenous respiration, and (3) effect

of zooplankton grazing. Peak resultant death

rate is 0.25/day at the beginning of August and
51

Figure 21. Lake 1 model verification, 0-17 meters.
52

Figure 22. Lake 1 model verification, 0-17 meters (continued).
53

Figure 23. Lake 1 model verification, 50-150 meters.
54

Table 9. Principal Parameter Values Used in Lake 1 Model Output (Shown in Figures 21-23)

Notation Name Value Unit

Kmn

Kmp

Cg

Kg

K2

azp

Kz

K1

K22

Kop

a1

ap

ac

w

Half-saturation constant-nitrogen

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

25

2

.06

10

0.1

0.6

.001

.00175

.002

.007

10

1

50

0.1

μg N/l

μG P/l

1/mg C-day-°C

μG Chlor/1

Days-1 (at 20°C)

-

(Days-°C)-1

(Days-°C)-1

(Days-°C)-1

(Days-°C)-1

μG N/μg Chlor

μG P/μg Chlor

μG C/μg Chlor

m/day
55

Figure 24. Dynamics of the phytoplankton growth rate.
56

Figure 25. Dynamics of the phytoplankton death rate.
57

is primarily due to zooplankton predation and

quantitatively explains the mid-summer decline

in phytoplankton biomass.

The dynamics of phytoplankton net production,

i.e., chlorophyll/liter-day, dP/dt, are shown in

Figure 26. Net production here includes the

kinetic interactions of growth in Figure 21, death

and predation in Figure 25, and the effect of

vertical mixing and lake outflow. The results

shown in Figure 23 summarize the basic

hypotheses of phytoplankton dynamics in Lake

Ontario.

The spring growth phase to approximately midMay is due primarily to increasing light and

temperature. There are little zooplankton 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 interactions which lead to a midsummer decline and subsequent broad fall peak.

It is hypothesized, as computed in Figure 26,

that during mid-summer, and following the

nutrient limitation effect, zooplankton grazing

becomes important and accounts for the

minimum values of biomass in July. However, at

increased zooplankton grazing, biological

recycling of nutrients becomes more significant.

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.

Figures 27 and 28 explore these dynamic effects

at another level of detail. The model permits

calculations of nutrient limitation effects as well

as the flux of nutrients due to zooplankton

excretion and release of nutrients due to

endogenous respiration by the plankton. Figure

27 shows the nutrient effects. As seen, the

spring bloom is halted essentially by a reduction

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

a limiting effect on growth. Figure 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).

Figure 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 summerearly 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 Figure 21). The results

shown in Figure 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,
58

Figure 26. Dynamics of the phytoplankton net production.
59

Figure 27. Dynamics of nitrogen and phosphorus limitation.
60

Figure 28. Biological and hypolimnetic recycling of phosphorus.
61

therefore, may not necessarily result in a

concomitant reduction in biomass. Further, the

time to reach a new equilibrium 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 behavior 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 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),

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 late

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 program is instituted. The long

detention times, recycling effects and large store

of nutrients in the hypolimnion are the principal

reason for this long response. Therefore, the

Lake Ontario presently being observed

represents conditions of about 10-20 years ago

and since loads into the upper basin (e.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 tenyear lag. The results also indicate 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.
62

CHAPTER 7

LAKE 2 ANALYTICAL INVESTIGATIONS

7.1 INTRODUCTION

As illustrated in Figure 4, a strategy of essentially

parallel investigations has been pursued in the

development of Lake Ontario phytoplankton

models. The Lake 2 investigation has results 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 distribution of phytoplankton

biomass have been known for some time (Riley

et al., 1949) and their applications to various

Great Lakes problem settings is ongoing (Di

Toro et al., 1973; Canale et al., 1973; Thomann

et al., 1973). 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 (Kierstead

and Slobodkin, 1953) 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 decrease. What has not been

recognized is that if the population does

increase, it asymptotically approaches a

condition which is a substantial simplification of

the general solution of the biomass conservation

equation. In particular, it can be shown that

there exists an asymptotic population growth

rate, μ, which is independent of position, and an

asymptotic population distribution, which

completely 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 that, in fact,

this prediction is born out by observed

phytoplankton biomass distributions.

The specific spatial setting for the analysis

presented subsequently is in the vertical

dimension. The motivation for the investigation

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 manageable level.

A final product of this analysis is a demonstration

that mathematical theory and biological fact can

be brought into accord if the complexities 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 a the primary dependent

variable is a biological example; the use of a

dispersion coefficient is a physical example. The
63

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.

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

 (56)

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:

 (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. The magnitude

of G is quite well known (Di Toro et al., 1971;

Eppley, 1972) from numerous short-term growth

experiments. For the case of excess nutrients

and constant optimal light intensity, G is in the

range of 1.5-2.5 day-1 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-1 at 20°C (Riley et

al., 1949). 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-1 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 reduce 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 extinction.

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 distribution will be of

the form:

 (58)

that is, the entire population will be growing at an

asymptotic population log growth rate, μ,

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 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 substantial insights into the

behavior of phytoplankton populations during

their growth phase.
64

7.3 ASYMPTOTIC POPULATION GROWTH

RATES

To establish that, in fact, phytoplankton biomass

behaves as predicted by the mathematical

analysis, a series of plots of log P(z,t) versus

time at various depths are presented in Figures

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 Figure

29. The data collected during the 1973 IFGYL

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/μ, 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 (Burns and Ross, 1972), the

phytoplankton population increased

approximately five-fold in 35 days. As shown in

Figure 30, the asymptotic population growth rate

was μ 0.043 day-1 and again the population

increased at this exponential rate at 1, 12 and

17-20 meters, with the population in the deeper

layers (21-25 meters) 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.

Onondaga Lake is a small (12 km2

 surface area,

20 m maximum depth), highly eutrophic lake in

New York State (O’Brien and Gere Engineers,

Inc., 1971). The duration of the growth period

for the phytoplankton population is 150 days

during which the population increases 1,000-fold

to concentrations of 105 cells/ml corresponding

to a population growth rate of μ = 0.048 day-1.

As shown in Figure 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 (Platt and Irwin,

1973). Similar behavior is observed during the

period of exponential growth. The population

growth rate is 0.16 day-1, 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 μ, ten 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 transport

parameters are constant, an asymptotic

population growth rate occurs which characterize

the growth of the entire population, independent

of position.

7.4 SOLUTIONS FOR SPECIFIC CASES

The existence of an asymptotic population

growth rate, independent 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 analytical solutions:

A two-layer approximation:

 (59)
65

Figure 29. Lake Ontario phytoplankton biomass, horizontal averages, IFYGL, 1973.
66

Figure 30. Central Lake Erie phytoplankton biomass, four station log averages, Project Hypo, 1970.
67

Figure 31. Onondaga Lake phytoplankton biomass, two stations, 1969.
68

Figure 32. St. Margaret’s Bay phytoplankton biomass, Station A, 1969.
69

and an exponential approximation:

 (60)

The solution proceeds directly using the

asymptotic solution (Equation 58). Substituting

it into the conservation of mass equation (56)

yields:

 (61)

This equation is solved together with the

boundary condition at the water surface:

 (62)

which requires that no biomass cross the airwater interface, and the boundary condition at

infinite depths, for which we require that the

biomass remain finite. The result is an

eigenvalue problem to find μ and P(z) for the

above ordinary differential equation with the

homogenous boundary conditions.

For the exponential case, the differential

equation (56) becomes:

 (63)

which has a solution (Murphy, 1960):

 (64)

where

and J<(z) is a Bessel function of the first kind. At

large depths, the solution becomes:

 (65)

which approaches zero exponentially. The

eigenvalue equation is obtained from the

boundary condition at z = 0 (Equation 62):

 (66)

A multiplicity of μ’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

certain sets of  and  no positive μ 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:

 (67)

where: 

and 

The solutions of these eigenvalue equations is

shown in Figure 33 in terms of the dimensionless

ratio of growth rate of the sum of death rate and

population growth rate G/D + ; a penetration

depth which characterizes the size of

the euphotic zone in terms of the growth rate

and the dispersion coefficient; and a peclet

number wL/E, which is the ratio of advective
70

Figure 33. Eigenvalue equation solutions for constant growth layer and exponentially decreasing

growth rate cases.
71

transport to dispersive 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, consider the effect of

settling velocity for the range of G/D +  < 10m,

it is clear that a peclet number, wL/E, of greater

than 0.1 is required to affect the growth rate of

the population.

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.

7.5 THE SINGLE LAYER APPROXIMATION

It is common practice in one-layer models

(Steele, 1956) 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 approximation comes from the

resulting differential equation for the euphotic

zone layer:

 (68)

or E’ the bulk mixing coefficient and P1 the

hypolimnion concentration. 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 Figure 34, a

plot of w/L versus G, both normalized by D + .

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,

Table 10. Log10 (G/D + ) for Exponentially Decreasing Growth Rate

Log10()

 - 1.0 - 0.9 - 0.8 - 0.7 - 0.6 - 0.5 - 0.4 - 0.3 - 0.2 - 0.1

0.50

0.20

0.10

0.05

0.02

0.01

0.00 2.0170

2.2545

1.8226

2.3352

1.8509

1.6325

1.7515

1.5721

1.4464

1.9903

1.4364

1.3456

1.2724

1.4123

1.2021

1.1512

1.1068

1.4158

1.1192

1.0105

0.9816

0.9542

1.7007

1.0455

0.9129

0.8519

0.8337

0.8166

1.0358

0.8251

0.7536

0.7172

0.7059

0.6951

2.0281

0.7731

0.6681

0.6263

0.6038

0.5966

0.5897

Log10()

 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

5.00

3.00

2.00

1.00

0.50

0.20

0.10

0.00

0.9140

0.6083

0.5486

0.4994

1.4778

0.6411

0.4908

0.4545

0.4228

0.7376

0.4880

0.4020

0.3789

0.3580

1.4192

0.5114

0.3856

0.3326

0.3175

0.3034

0.6510

0.3842

0.3113

0.2771

0.2670

0.2575

0.7034

0.4352

0.3002

0.2548

0.2322

0.2253

0.2187

1.2401

0.4394

0.3186

0.2401

0.2106

0.1953

0.1906

0.1860

0.5484

0.3106

0.2440

0.1950

0.1754

0.1649

0.1615

0.1583

0.3512

0.2322

0.1920

0.1602

0.1469

0.1395

0.1372

0.1349

0.2479

0.0795

0.1540

0.1327

0.1235

0.1183

0.1167

0.1150
72

Figure 34. Settling velocity-growth rate interaction - the effect of dispersion. Exponentially

decreasing growth rate case.
73

assume that D +  and L are fixed, and the plot

is w/L versus G for varying dispersion coefficient.

For large vertical dispersion, 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

 constant. What is occurring is that the settling

biomass is being rapidly remixed into the

euphotic zone so that the loss 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 one,

indicating that for this situation it requires larger

than proportional increases in G to compensate

for increases in w/L. The meaning of this

unexpected result will become clear

subsequently.

7.6 THE EFFECT OF INCREASING GROWTH

RATE

For various euphotic zone depths (assuming the

transport parameter w and E are constant) the

effect on  of increasing G is shown in Figure

35. It is immediately clear that no solutions are

possible for which  > 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 the 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

+ . As G is further increased, however, the

population growth rate asymptotically

approaches G - D as expected, and since there

are nonzero values of G for which D +  is zero,

it is reasonable that the slopes are greater than

one to enable D +  to, in effect, catch up with

G. The condition G = D +  occurs at smaller G

for deeper euphotic zone, as expected.

7.7 THE EFFECT OF INCREASING SETTLING

VELOCITY

For various euphotic zone depths (assuming G

and E are constant) the effect on  of increasing

w is shown in Figure 36. As w is increased, the

population growth rate is not affected until a

critical range is reached at which time the effect

begins. Further increases cause a sharp drop

off in  until population growth is no longer

possible.

7.8 THE EFFECT OF INCREASING

DISPERSION

The effect on  of increasing the dispersion

coefficient, E, for various euphotic zone depths

(assuming G and w are constant) is shown in

Figure 37. A most surprising result is that for the

range of the parameter, EG/w2 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 population 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 Figure 34 for small

dispersion coefficients.

7.9 THE EFFECT OF BUOYANCY

It is well known that certain species of

phytoplankton exhibit negative sinking velocity

and an in situ measurement has confirmed an

example (Burns and Pashley, 1974) so that

instead of a liability their buoyancy is an asset

which should increase population growth. Figure

38 presents the results. For negative peclet

numbers with magnitude less than 0.2, not much

increase in  is calculated. However, as |w|L/E

increases further, the population growth rate

approaches the maximum possible, i.e.,  + G -

D, independent of as indicated by the

horizontal lines in the figure. The effect of the

upward advection is to pile the biomass at the
74

Figure 35. Population growth rate versus maximum growth rate - the effect of euphotic zone. The

exponentially decreasing growth rate case.
75

Figure 36. Population growth rate verus settling velocity - the effect of euphotic zone depth. The

exponentially decreasing growth rate case.
76

Figure 37. Population growth rate versus dispersion coefficient - effect of euphotic zone depth. The

exponentially decreasing growth rate case.
77

Figure 38. The effect of negative sinking velocity. The exponentially decreasing growth rate case.
78

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 exceeds 0.3.

7.10 PRELIMINARY APPLICATIONS

As a first step in exploring the utility of the

preceding analysis, a series of representative

values for the dimensionless parameters are

shown in Figure 39. The majority come from the

coastal ocean (Riley et al., 1949) with two lakes

included. It is interesting to note that the values

of G/D +  are comparable whereas the value of

 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;  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 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 pursue this point.

7.11 SUMMARY OF MATHEMATICAL FACTS

The properties of Equation (56) are best

illustrated by examining a spatially discretized

version.

Let Pi

(t) = P[(i + 1⁄2) z,t] for a finite set of points

starting at z/2 and spaced z apart thereafter.

The spatially discrete version of the differential

equation follows by using the approximations:

 (69)

 (70)

There results a set of ordinary differential

equation which can be written in matrix form as:

 (71)

Two properties of A are important:

1. A is essentially nonnegative, i.e.,

aij > 0 I  j (72)

2. A is irreducible, i.e., for any i,j there exists

a finite sequence of indices k(0) = k, k(1),

..., k(r) = j, such that ak(h-1),k(h)  0 for h=1,...,

r.

Property (1) follows directly from the fact that the

off diagonal elements of A are made up of sums

of E/z2 and w/z which are positive. Property

(2) essentially requires that mass from any

segment 1 can eventually be transported to any

segment j, which is clear from the geometry of

any reasonable problem.

The theorem (Birkoff and Varga, 1958) of

interest is:

The matrix A satisfying (1) and (2) has a

unique (up to a multiplicative constant)

strictly positive eigenvector 1, with real

simple eigenvalue 1, i.e., 

 (73)
79

Figure 39. Application to lakes and oceans - the range of biologically relevant dimensionless

numbers.
80

and 1 is the largest eigenvalue in the

sense that 

 (74)

for all the other eigenvalues, j, of A.

It follows directly from this fact that the

asymptotic solution of the differential equation

(71) has the form:

 (75)

where

 

and k is a constant.

This fact can be proven directly (Birkoff and

Varga, 1958). However, to see what is

occurring assume that all eigenvalue j

 are

distinct so that n eigenvalue of the solutions

exists (Bellman, 1960) of the form.

 (76)

where kj

 are constants related to the initial

condition P(z,0). If 1 is positive and Re [j

]

negative for all other j, then clearly the positive

eigenvalue part of the solution dominates after

a time on the order of 1/1. If both 1 and 2 are

positive (and 2 real for simplicity), the solution

in segment i approaches:

 (77)

but 2 < 1 so that the exponent is negative and

it becomes small relative to one after a time on

the order of 1/(1 - 2) and again the maximum

eigenvalue solution dominates.
81

CHAPTER 8

PRELIMINARY LAKE 3 MODEL

As shown in Figure 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 threedimensional, time-variable model which does

permit some analysis of phenomena in the

nearshore region as compared to the open lake

region. The structure of Lake 3 represents a

further step beyond the simple Lake 1 model, the

extent 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 details 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 Figure 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 nearshore

environment. All cross-sectional areas for

interfaces in three-dimensions between the

segments were planimetered and volumes were

computed and summed to within 10% of total

lake volume.

8.1 CIRCULATION AND DISPERSION

The emphasis in Lake 3 is on the biological and

chemical kinetics rather than the hydrodynamic

circulation, although the interaction between the

kinetics and transport is important. The

hydrodynamic circulation is, therefore,

considered to be externally supplied through

inferences from observations and/or

hydrodynamic model output. An extensive

literature exists on the circulation of Lake Ontario

and relatively sophisticated mathematical models

have been developed (Simons, 1971, 1972,

1973). A review of the complete development is

given by Simons (1973). 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 discussion will, therefore, outline the

procedure used for establishing circulation

patterns and flow exchanges in the 67-segment,

five-layer, Lake 3 model. Two flow balances are

presented, one corresponding to “summer

circulation” and the other to “winter circulation”.

It is known that Lake Ontario has two distinct net

surface circulation patterns. During summer and

fall the prevailing winds, the dominant factor

governing water circulation, are from the west
82

Figure 40. Lake 3 segmentation.
83

and southwest and one distinct circulation

pattern is set up (“summer circulation”). During

winter and spring, the prevailing 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 with

the stratified and isothermal temperature states

of the water, respectively.

The IJC (1969) and the USEPA (Casey, 1965)

have published summer circulation patterns for

Lake Ontario. Mean monthly resultant velocity

vectors that were observed and used by USEPA

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 (1969) as a guideline.

The water budget information contained in the

IJC report was used for establishing tributary

inflows to the various shoreline segments. It was

decided to split the major inflow to Lake Ontario

from the Niagara River (195,000 cfs) over the

first and second vertical layers (0-4 and 4-17

meters) while restricting all other tributary inflows

(Oswego River, Genesse River, Twelve Mile

Creek, etc.) to the first layer (0-4 meters). The

Niagara inflow was split volumetrically between

segments 5 and 32 each receiving 4/17 and

10/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 by IJC. The final

circulation pattern, as given by intersegment

velocity values can be seen in Figure 41. It

should be noted that since the Niagara River

inflow was split between 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 the top

two layers), the outflow segment, is only 20-25

meters while the third layer extends 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 another 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 USEPA general circulation

pattern corresponding to winds out of the west

was used for establishing a flow balance.

Bonham-Carter et al. (1973) wind-driven

numerical model outputs a fairly similar

circulation for a west wind. The model, however,

indicates easterly flow along the north shore,

while the USEPA 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 USEPA 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 indicate how Bonham-Carter et al.

(1973) modeled circulation was set up to this

segmentation scheme. The balanced velocity

vectors are 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 5 km) grid on to the

larger grid of the phytoplankton model. As

additional information becomes available,

adjustments to eh circulation structure can be

made.
84

Figure 41. Assumed summer circulation for Lake 3 model.
85

Figure 42. Assumed winter circulation for Lake 3 model.
86

Time and space varying horizontal and vertical

dispersion coefficients are incorporated to

simulate some of the more predominant features

such as the thermal bar and vertical stratification.

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

Figure 43 illustrates the time varying horizontal

and vertical exchange coefficients to be used in

this model.

8.2 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 reasonable in its broadest outlines.

Since the segments of the Lake 3 model are

large (10-40 km), 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 required 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 temperature

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 (1971) for IFYGL, who

have computed values of monthly heat flux

terms for 1965 through 1968. Though

discrepancies between their results and

observed data are reported, the computation of

a total heat storage term, Qt

, proves to be ideal

for coupling to the Lake 3 model and permits

inputting temperature as a lake surface forcing

function.

Using Lake Ontario meteorological data and

empirical relationships, then heat storage values

were computed using the formula:

 (78)

Qt

 is the heat storage, Qs the incident solar

radiation, Qb the net back radiation, Qh the heat

transfer to atmosphere and Qe is the evaporation

energy where all units are cal/m2

-day. With the

monthly heat flux terms, values of Qt

 were

calculated for the five-year period studied.

These values are listed in Table 11 and the fiveyear average values were used in the Lake 3

runs.

Once Qt, the total heat flux term or total heat

storage in the lake, is known, then

 (79)

where Ts(t) is the temperature forcing function for

the surface layer (°C/day or cal/cm3-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

cmputed variable temperature response in each

of the segments was then compared with a
87

Figure 43. Assumed horizontal and vertical dispersion coefficients for Lake 3 model.
88

Table 11. Monthly Heat Flux Values

Month 1965 Qt* 1966 Qt 1967 Qt 1968 Qt 1969 Qt

Five-Year

Average

January

February

March

April

May

June

July

August

September

October

November

December

-329.1

-216.6

 4.6

 275.2

 533.2

 463.0

 270.7

 67.0

 -43.7

-235.8

-200.6

-226.1

-392.0

-139.5

 80.3

 207.0

 413.7

 610.9

 307.0

 69.5

 -55.2

 -64.4

-162.4

-349.1

-214.2

-278.1

 48.4

 318.8

 420.6

 618.6

 303.0

 51.2

 -92.7

-108.7

-254.8

-296.8

-329.2

-244.9

 88.4

 361.1

 353.7

 471.7

 293.3

 198.5

 102.1

 -46.9

-206.7

-370.2

-275.5

-153.6

 11.3

 313.1

 448.1

 475.4

 -

 -

 -

 -

 -

 -

-308.0

-206.54

 46.6

 295.04

 380.8

 464.8

 293.5

 96.55

 -22.375

-113.95

-206.125

-310.55

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 interfaces and its

142 flow interfaces, is extensive so that a

transport regime based on verification against

generalized temperature 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 was the temperature

characterization of Lee (1972). Composite

temperature 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

nearshore 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 dataset for temperature 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 Figures 45 and

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

the flow field were investigated.
89

Figure 44. Regional division of Lake Ontario by Lee (1972).
90

Figure 45. Comparison of computed temperatures from Lake 3 model with data from Lee (1972),

winter and fall.
91

Figure 46. Comparison of computed temperatures from Lake 3 model with data from Lee (1972),

spring.
92

Table 12. Regional to Lake 3 Segmentation Equivalence

Segment of Lake 3

Lee (Region) 1 2 3 4 5

1

2

3

4

5

6

7 

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

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

60

53,55,56

56,58

60,62

54,57

57,61

55,59 

 -

 -

 -

 -

63,64

64,65

 -

 -

 -

 -

 -

66

66,67

 -

As can be seen, no marked effect on the

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

nearshore to main lake resolution, the surface

areas of the segments are so large that the

interfaces for vertical transport are much greater

than those for horizontal and flow transport.

Hence, vertical dispersion terms dominate.

The grid is not fine enough and effects such as

nearshore thermal bar flow restriction and

changing wind pattern effects on horizontal

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

nearshore to mid-lake segments horizontal

dispersion for a period in the spring to simulate

thermal bar, and simulating stratification earlier

in the nearshore 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

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

8.3 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 Chapter 5. The purpose of these runs is

to determine the general direction of the output

with special reference to spatial differences in

phytoplankton biomass. Considerably 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 onto the Lake 3

grid and only those stations for which a

reasonably complete set of data were available
93

are included in the mapping. Figure 47 shows

the location of the data stations and the Lake 3

grid. Some segments have only one station for

comparison while other segments (primarily

nearshore stations) may have three stations for

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

Figure 48 shows very early 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 embayment and as seen reaches

peak values of almost 15 μg chlorophyll/l 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 10 μg/l is computed. This

results primarily from the thermal bar effect

simulated by variations in the horizontal

dispersion (see Figure 43). Gradients across the

lake in the fall are not pronounced and reflect the

general well-mixed character of the lake at that

time.

Figure 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 Figures 48 and 49 is for

zero phytoplankton settling velocity.

8.3.1 Output Data Display

The Lake 3 model generates a substantial output

totaling 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 comprehend 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.

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

preliminary run for a full year of the surface

phytoplankton chlorophyll.
94

Figure 47. Principal IFYGL stations and the Lake 3 grid.
95

Figure 48. Preliminary results of Lake 3 model, 0-4 meters.
96

Figure 49. Preliminary comparison of Lake 3 model output to observed data at two segments for

June.
97

Figure 50. Data output flow diagram used for verification and display purposes.
98

Figure 51. Three-dimensional plot of phytoplankton chlorophyll calculated from Lake 3 model in

June, 0-4 meters.
99

APPENDIX A

LAKE 1 - INPUT AND PROGRAM LISTING

A.1 INPUT INFORMATION

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 Figure 9.

The purpose of this appendix is to sumarize the

input information needed to run the model and to

present, where available, the background

information from which this information was

drawn. Units for the information given are

consistent with those required by the

computational software used, which is named

Water Analysis Simulation (WASP). A listing of

WASP and sample output are included.

A.2 MORPHOMETRY AND HYDRODYNAMIC

REGIME

The information required to describe the

morphometry of Lake Ontario for the Lake 1

model is presented in Table A-1, with the

associated hydrodynamic regime.

The bulk of the morphometric data was

generated using data supplied by the CCIW

(personal communication). The bathymetric

data, supplied on punched cards, were depth

averages for a square grid laid on a polyconic

projection of Lake Ontario, with geographical

areas of 4 km (IJC, 1969). These surface areas

were defined by sides which were actually

curves on the earth’s surface. These data were

processed to yield the necessary information.

The average depth of Lake Ontario is 396 feet

(90 meters). The IJC (1969) 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-1 and A-2

which show the variation of volume and surface

area with depth were generated with the CCIW

supplied data.

Table A-1. Physical Parameters

SegmentInterSegmentRepresentationVolume

(106

 cf)

Depth

(ft)

Intersegment

Areas (ft2)

Inflow

(cfs)

Outflow

(cfs)

Dispersion

Coefficient

(mi2/day)

1

2

3

1-2

2-3

Epilimnion

Hypolimnion

Benthos

1.05 x 107

4.85 x 107

4.80 x 104

55.8

240.5

0.5

1.76 x 1011

0.959 x 1011

43,.500.

188,400.

0.0

43,500.

188,500.

0.0

See Figure A-4

0.0
100

Figure A-1. Variation of volume with depth.

Figure A-2. Variation of area with depth.
101

The hydrodynamic regime is inputted into

WASP, as opposed to a modeling configuration

which would incorporate a hydrodynamic

submodel internally. This differentiation is most

significant for more spatially defined models than

Lake 1. The hydrodynamics must, therefore, be

obtained from a database, the literature, or an

external hydrodynamic submodel. The mean

annual outflow for Lake Ontario is 232,000 cfs

(IJC, 1969). Inflow was set equal to outflow and

was split proportionally to the segment depths,

as shown in Table A-1.

The development of the thermocline was

simulated by temporally varying the vertical

dispersion coefficient as shown in Figure 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 dispersion coefficient

was set at 0.0003 mi2

/day (90 cm2/sec) for

periods of non-stratification, which is in the range

reported in Hydroscience, Inc. (1973).

A.3 PHYSICAL EXOGENOUS VARIABLES

Water Temperature - The temporal variation of

water temperature used in Lake 1 is shown in

Figure A-4. Since the segments represent

completely mixed volumes, representative

temperatures were chosen. Data averaged from

CCIW (1969) were used as temperature input.

Only stations with depths over 50 meters were

considered. These were viewed as

representative of main lake conditions.

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 overland radiation to overlake

radiation ratio proposed by Richards and Lowen

(1965). The overland radiation data, lake-to-land

ratio and adjusted temporal lake radiation

function used is shown in Figure A-5.

Photoperiod - That fraction of the day from

sunrise to sunset is variable during the year, and

is a function of latitudinal position. Sunsetsunrise data from Newspaper Enterprise

Association (1974) were used to obtain the

functions shown in Figure A-6. The Lake Ontario

photoperiod function was obtained by

interpolating to the equivalent of 43°40'N

latitude.

Background Water Clarity - The variation of the

extinction coefficient with time is a function of the

phytoplankton concentration in the water body.

Given no phytoplankton, there 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 Figure A-7. The extinction coefficient

was taken to be:

after Beeton (1958).

A.4 BIOLOGICAL AND CHEMICAL

EXOGENOUS VARIABLES

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 situations with the “downstream”

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

Figure A-5. Temporal variation of solar radiation.
104

Figure A-6. Temporal variation of photoperiod.
105

Figure A-7. Determination of background water clarity - Ke.
106

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 phosphorus

species used in the nutrient systems (non-living

organic, ammonia and nitrate, nitrogen, nonliving organic and available phosphorus) was not

available. Only total nutrient information was

given. Species 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.

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

A.5 CHEMICAL AND BIOLOGICAL

ENDOGENOUS VARIABLES

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

the present configuration of Lake 1.

NOTE:

Source listings for Lake 2 model and sample

input and output are available from USEPA,

Grosse Ile Laboratory, Grosse Ile, Michigan.

This program is currently operational on the

Control Data Corporation 6600 at Courant

Institute of New York University.

A version will be available on Optimal systems,

Incorporated, in the near future.

Table A-2. Lake 1 - Mass Loadings and Initial Conditions

System Boundary Concentration (mg/L)

Initial Conditions

(mg/l)

Phytoplankton

Zooplankton - Herbivorous

Nitrogen

 Non-Living Organic

 Ammonia-N

 Nitrate-N

Phosphorus

 Non-Living Organic

 Available P

Zooplankton - Carnivorous

Zooplankton - Upper Trophic Level 1

Zooplankton - Upper Trophic Level 2

1.0

.05

.443

.0150

.250

.0454

.0143

.01

By-passed system

By-passed system

2.0

.005

.090

.0150

.250

.005

.0143

.005

I
107

Table A-3. System Parameters

System - Parameter Symbol Value Units

Phytoplankton Chlorophyll a

 1. Saturated Growth Rate (	o)

 2. Saturating Light Intensity

 3. Michaelis Constant for Nitrogen

 4. Michaelis Constant for Phosphorus

 5. Endogenous Respiration Rate (	20)

 6. Carbon-to-Chlorophyll Ratio

 7. Sink Velocity

 8. Preference Structure for Inorganic Nitrogen

K1T

K1C

IS

KMN

KMP

K2T

K2C

CCHL

SVEL

PNH3

1.066

0.580

350.

0.025

0.002

1.08

0.10

0.05

0.10

0.5

None

day-1

ly/day

mg N/l

mg P/l

None

day-1

mg C/g Chl a

m/day

None

Non-Living Organic Nitrogen

 1. Organic N-NH3 Hydrolysis Rate

 2. Nitrogen-to-Chlorophyll a Ratio

 3. Sink Coefficient

Ammonia Nitrogen

 1. NH3-NO3 Nitrification Rate

Nitrate-Nitrogen

 None

K34T

K34C

NCHL

K33

K45T

K45C

0.00175

0

0.010

0.001

0.002

0

day °C-1

day-1

mg N/mg Chl a

day-1

day °C-1

day-1

Non-Living Organic Phosphorus

 1. Organic P-Inorganic P Conversion Rate

 2. Phosphorus-to-Chlorophyll a Ratio

 3. Sink Coefficient

K67T

K67C

PCHL

K66

0.007

0

0.001

0.001

day °C-1

day-1

mg P/μg Chl a

day-1
108

Table A-3. System Parameters (Continued)

System - Parameter Symbol Value Units

Available Phosphorus

 1. Sinking Coefficient

Herbivorous Zooplankton-Carbon

 1. Grazing Rate

 2. Michaelis Constant for Phytoplankton

 3. Conversion Efficiency

 4. Endogenous Respiration Rate

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

Upper Trophic Level 2 - Carbon

 1. Grazing Rate

 2. Conversion Efficiency

 3. Endogenous Respiration Rate

K77

CGC

CGT

KMP

AZP

K4T

K4C

CGT8

CGC8

A8Z

K8T

K8

CGT9

CGC9

A98

K9T

K9

CGT0

CGC0

A109

K10T

K10

K9

0

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

0.06

0

0.60

0.002

0

day-1

l/mg C/day °C

l/mg C/day

μg Chl a/l

None

day °C-1

day-1

l/mg C/day °C

l/mg C/day

None

day °C-1

day-1

l/mg C/day °C

l/mg C/day

None

day °C-1

day-1

1/mg C/day °C

1 mg C/day °C

None

day °C-1

day-1
109

REFERENCES

Beeton, A.M. 1958. Relationship Between

Secchi Disc Readings and Light Penetration in

Lake Huron. Amer. Fish. Soc. Transac., 87:73-

79.

Beeton, A.M. and D.C. Chandler. 1963. The St.

Lawrence Great Lakes. In: D.C. Frey (Ed.),

Limnology in North America. University of

Wisconsin Press, Milwaukee, Wisconsin.

Bellman, R. 1960. Introduction to Matrix

Analysis. McGraw-Hill Book Company,

Incorporated, New York, New York.

Birkhoff, G. and R.S. Varga. 1958. Reactor

Criticability and Nonnegative Matrices. J. SIAM,

6(4):354-377.

Bonham, Carter, G., J.H. Thomas, and D.

Lackner. 1973. A Numerical Model of Steady

Wind-Driven Currents in Lake Ontario and the

Rochester Embayment Based on Shallow-Lake

Theory, International Field Year for the Great

Lakes Report Number 7, Rochester Embayment

Project. Department of Geological Sciences,

University of Rochester, Rochester, New York.

Burns, N.M. and C. Ross. 1972. Project Hypo.

Canada Centre for Inland Waters Publication

Number 6. Report to the U.S. Environmental

Protection Agency, Washington, D.C. EPA-TS05-71-208-24.

Burns, N.M. and A.E. Pashley. 1974. In Situ

Measurements of the Settling Velocity Profile of

Particulate Organic Carbon in Lake Ontario. J.

Fish. Res. Bd. Canada, 31(3):291-297.

Canada Centre for Inland Waters (CCIW).

Personal Communication on Descriptive

Limnology. Burlington, Ontario, Canada.

Canada Centre for Inland Waters (CCIW). 1969.

Limnological Data Reports, Lake Ontario, 1966-

1969. Canadian Oceanographic Data Centre,

Burlington, Ontario, Canada.

Canadian Hydrographic Service. 1970.

Bathymetric Chart (881) of Lake Ontario.

Department of Energy, Mines and Resources,

Ottawa, Ontario, Canada.

Canale, R.P., S. Nachiappan, D.J. Hineman, and

H.E. Allen. 1973. A Dynamic Model for

Phytoplankton Production in Grand Traverse

Bay. In: N.A. Rukavina, J.S. Seddon, and P.

Casey (Eds.), Proceedings of the 16th

Conference on Great Lakes Research, pp.

21-23. Braun-Brumfield Publishers, Ann

Arbor, Michigan.

Caperon, J. and J. Meyer. 1972. NitrogenLimited Growth of Marine Phytoplankton. II.

Uptake Kinetics and Their Role in NutrientLimited Growth of Phytoplankton. Deep Sea

Res., 19:619-632.

Casey, D.J. 1965. Lake Ontario Environmental

Summary. Report to the Lake Ontario Basin

Office, U.S. Environmental Protection Agency,

Rochester, New York.
110

Di Toro, D.M., D.J. O’Connor, and R. Thomann.

1971. A Dynamic Model of the Phytoplankton

Population in the Sacramento-San Joaquin

Delta. In: Advances in Chemistry, Number 106,

pp. 131-180. American Chemical Society,

Washington, D.C.

Di Toro, D.M., D.J. O’Connor, R.V. Thomann,

and J.L. Mancini. 1973. Preliminary

Phytoplankton, Zooplankton Nutrient Model of

Western Lake Erie. In: B.C. Patten (Ed.),

Systems Analysis and Simulation in Ecology,

Volume 3. Academic Press, Incorporated, New

York, New York.

Di Toro, D.M. 1974a. An Evaluation of

Phytoplankton Kinetic Behavior in Flask

Experiments. In preparation.

Di Toro, D.M. 1974b. Combining Chemical

Equilibrium and Water Quality Models.

Manhattan College, Riverdale, New York.

 

Epply, R.W. 1972. Temperature and

Phytoplankton Growth in the Sea. Fish. Bullet.,

70(4):1063-1085.

Glooschenko, W., J.E. Moore, and R.A.

Vollenweider. 1972. The Seasonal Cycle of

Pheo-Pigments in Lake Ontario With Particular

Emphasis on the Role of Zooplankton Grazing.

Limnol. Oceanogr., 17(4):597-605.

Glooschenko, W.A., J.E. Moore, M. Munawar,

and R.A. Vollenweider. 1974. Primary

Production in Lakes Ontario and Erie: A

Comparative Study. J. Fish. Res. Bd. Canada,

31(3):253-263.

Goering, J.J., D.N. Nelson, and J.A. Carter.

1973. Silicic Acid Uptake by Natural Populations

of Marine Phytoplankton. Deep Sea Res.,

20(9):777-789.

Great Lakes Basin Commission. 1971.

Appendix II - Levels and Flows. In: Great Lakes

Framework Study, Draft Number 3. Great Lakes

Basin Commission, Ann Arbor, Michigan.

Helgeson, H.C. 1967. Thermodynamics of

Complex Dissociation in Aqueous Solution at

Elevated Temperatures. J. Phys. Chem.,

71:3121-3136.

Hendrey, G.R. and E. Welch. 1973. The Effects

of Nutrient Availability and Light Intensity on the

Growth Kinetics of Natural Phytoplankton

Communities. Presented at the 36th Conference

of the American Society of Limnology and

Oceanography, Salt Lake City, Utah. June 1973.

Hubbard, J.E. Personal communication.

Department of Geology and Earth Science, State

University College, Brockport, New York.

Hutchinson, G.E. 1957. A Treatise on

Limnology, Volume 1 - Geography, Pysics and

Chemistry. John Wiley and Sons, Incorporated,

New York, New York.

Hutchinson, G.E. 1967. A Treatise on

Limnology, Volume 2 - Introduction to Lake

Biology and the Limnoplankton. John Wiley and

Sons, Incorporated, New York, New York.

Hydroscience, Incorporated. 1973. Limnological

Systems Analysis of the Great Lakes. Report to

the Great Lakes Basin Commission, Westwood,

New Jersey. 474 pp.

International Joint Commission. 1969.

International Lake Erie and 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

Rivers, Volumes 1, 2, and 3. Windsor, Ontario,

Canada.

Kern, D.M. 1960. The Hydration of Carbon

Dioxide. J. Chem., 37(1):14-22.

Kibby, N.V. 1971. Effects of Temperature on

the Feeding Behavior of Daphnia rosea. Limnol.

Oceanogr., 16(3):580-581.
111

Kierstead, H. and L.B. Slobodkin. 1953. The

Size of Water Masses Containing Plankton

Blooms. J. Mar. Res., 12:141-147.

Kramer, J.R. 1964. Theoretical Model for the

Chemical Composition of Fresh Water With

Application to the Great Lakes. Great Lakes

Research Division, Publication Number 11, pp.

149-160. The University of Michigan, Ann Arbor,

Michigan.

Kramer, J.R. 1967. Equilibrium Models and

Composition of the Great Lakes. In: Equilibrium

Concepts in Natural Water Systems, Advances

in Chemistry, Series Number 67, pp. 243-254.

American Chemical Society, Washington, D.C.

Lee, A.H. 1972. Regional Characterizations of

the Thermal Properties of Lake Ontario. In:

D.V. Anderson and J.S. Seddon (Eds.),

Proceedings of the 15th Conference on Great

Lakes Research, pp. 625-634. BraunBrumfield Publishers, Ann Arbor, Michigan.

McNaught, D.C. and M. Buzzard. 1973.

Zooplankton Production in Lake Ontario as

Influenced by Environmental Perturbations. In:

First Annual Reports of the EPA/IFYGL Projects,

pp. 28-70. National Environmental Research

Center, Office of Research and Development,

U.S. Environmental Protection Agency, Corvallis,

Oregon. EPA-600/3-73-021, 336 pp.

Murphy, G.M. 1960. Ordinary, Differential

Equations and Their Solutions. D. Van Nostrand

Corporation, Princeton, New Jersey.

Murthy, C.R., G. Kullenberg, H. Westerberg, and

K.C. Miners. 1974. Large Scale Diffusion

Studies. In: International Field Year for the

Great Lakes Bulletin Number 10, pp. 22-49.

U.S. Department of Commerce, National

Oceanic and Atmospheric Administration,

Rockville, Maryland.

Newspaper Enterprise Association. 1974. The

World Almanac and Book of Facts. Newspaper

Enterprise Association, New York, New York.

O’Brien and Gere Engineers, Incorporated.

1971. Onondaga Lake Study. Water Quality

Office, U.S. Environmental Protection Agency,

Washington, D.C. 11060 FAE 04/71, 460 pp.

Paasche, E. 1973. Silicon and the Ecology of

Marine Plankton Diatoms. I. Thalassiosira

pseudomona (Cyclotella nana) Grown in a

Chemostat With Silicate as Limiting Nutrient.

Marine Biol., 19:117-126.

Patalas, K. 1969. Composition and Horizontal

Distribution of Crustacean Plankton in Lake

Ontario. J. Fish. Res. Bd. Canada, 26(8):2135-

2164.

Phillips, D.W. and J.A.W. McCulloch. 1972.

The Climate of the Great Lakes Basin.

Climatological Study Number 20, Environment

Canada, Toronto, Ontario, Canada.

Platt, T. and B. Irwin. 1973. Preliminary

Productivity Measurements in St. Margaret’s

Bay. Technical Report Number 203. Fisheries

Research Board of Canada, Ontario, Canada.

Porter, K.G. 1973. Selective Grazing and

Differential Digestion of Algae by Zooplankton.

Nature, 244(5412)179-180.

Richards, T.L. and P. Lowen. 1965. A

Preliminary Investigation of Solar Radiation Over

the Great Lakes as Compared to Adjacent Land

Areas. In: Proceedings of the 8th Conference

on Great Lakes Research, pp. 278-282.

Braun-Brumfield Publishers, Ann Arbor,

Michigan.

Riley, G.A.; H. Stommel, and D.F. Bumpus.

1949. Quantitative Ecology of the Plankton of

the Western North Atlantic. Bull. Bingham

Oceanogr. Coll., 12(3):1-169.

Rodgers, G.K. and G.K. Sato. 1971. Energy

Budget for Lake Ontario. Canadian

Meteorological Research Report, Atmospheric

Environment Service, Ontario, Canada.
112

Shapiro, N. 1962. Analysis by Migration in the

Presence of Chemical Reaction. Report Number

P-2596, Rand Corporation, Santa Monica,

California.

Shapley, M. and L. Cutler. 1970. Rand’s

Chemical Composition Program: A Manual.

Report Number R-495-PR, Rand Corporation,

Santa Monica, California.

Simons, T.J. 1971. Development of Numerical

Models of Lake Ontario. In: D.V. Anderson

and J.S. Seddon (Eds.), Proceedings of the

14th Conference on Great Lakes Research,

pp. 654-669. Braun-Brumfield Publishers,

Ann Arbor, Michigan. 

Simons, T.J. 1972. Development of Numerical

Models of Lake Ontario, Part 2. In: D.V.

Anderson and J.S. Seddon (Eds.),

Proceedings of the 15th Conference on Great

Lakes Research, pp. 655-672. BraunBrumfield Publishers, Ann Arbor, Michigan.

Simons, T.J. 1973. Development of a

Three-Dimensional Numerical Model of the

Great Lakes. Scientific Series Number 12.

Canada Centre for Inland Waters, Burlington,

Ontario, Canada. 26 pp.

Steele, J.H. 1956. Plant Production on Fladen

Ground. J. Mar. Biol. Assoc., 35:1-33.

Stumm, W. and J.J. Morgan. 1970. Aquatic

Chemistry. John Wiley and Sons, Incorporated,

New York, New York.

Thomann, R.V., D.M. Di Toro, D.J. O’Connor,

and R.P. Winfield. 1973. Mathematical

Modeling of Eutrophication of Large Lakes. In:

N.A. Thomas (Ed.), First Annual Report of the

EPA IFYGL Projects. U.S. Environmental

Protection Agency, National Environmental

Research Center, Corvallis, Oregon. EPA660/3-73-024. 336 pp.

Thomann, R.V., D.M. Di Toro, and D.J.

O’Connor. 1974. Preliminary Model of Potomac

Estuary Phytoplankton. J. Environ. Eng., ASCE,

100(SA3):699-715.

Trussell, R.R. and J.F. Thomas. 1971. A

Discussion of the Chemical Character of Water

Mixtures. J. Amer. Water Works Assoc.,

63(1):49-51.

Wagman, D.C., W.H. Evans, V.B. Parker, I.

Halow, S.M. Bailey, and R.H. Schumm. 1968.

Selected Values of Chemical Thermodynamic

Properties. Technical Note 270-3, National

Bureau of Standards, Washington, D.C.

Watson, N.H.F. and G.F. Carpenter. 1974.

Seasonal Abundance of Crustacean

Zooplankton and Net Plankton Biomass of Lake

Huron, Erie and Ontario. J. Fish. Res. Bd.

Canada, 31(3):309-317.

Yentsch, C.S. 1962. Marine Plankton. In: R.A.

Lewin (Ed.), Physiology and Biochemistry of

Algae, pp. 771-797. Academic Press, New York,

New York.
113