- Email: [email protected]
Interfacing Physical and Biochemical Models of Large Lakes
INTERFACING PHYSICAL AND BIOCHEMICAL MODELS OF LARGE LAKES
T.J. Simons Canada Centre for Inland Waters P.O. Box 5050 Burlington, Ontario L7R 4A6
ABSTRACT Comprehensive models for water quality control and management of large lakes have two identifiable components, namely a physical submodel and a package of biochemical routines. The physical model component supplies all essential information on water transports in the basin and describes the distribution of temperature and radiation required for a simulation of photosynthesis and other biological processes. The biochemical submodel is concerned with the various chemical and biological interactions resulting in a local increase of one biochemical variable at the expense of others. The two major components of the water quality model are tied together by a computational framework based essentially on a generalized advection-diffusion equation which expresses the principle of mass conservation as applied to each model variable. An essential property of such interdisciplinary model is a modular structure allowing for modifications pertaining to specific hydrodynamic or biological processes without affecting the overall computer program . This paper discusses the different components of current water quality models of the North American Great Lakes with particular emphasis on the transport mechanisms acting in such lakes .
the use of modern-day high-speed computers. Prerequisites for such investigations are obviously an adequate physical lake model, a realistic biochemical model, and a practical transport framework to couple these submodels. This paper will discuss some aspects of this interfacing problem followed by a brief description of the pertinent submodels . MODEL FRAMEWORK
I NTRODUCTI ON
The purpose of this effort is to construct a computational framework to serve as a basis of interaction and synthesis for an interdisciplinary group of scientists involved with water quality models. The underlying mathematical principles are that the models are prognostic in character, such that all variables are computed forward in time from specified initial values, and all changes of model variables are computed from mass conservation principles, regardless of the type of natural process to be simulated. From a computer-technical viewpoint, the conditions are that scientific concepts regarding specific hydrodynamic or biological processes can be incorporated with no more than a superficial knowledge of submodels pertaining to other disciplines. Thus the model is constructed in a modular fashion, characterized by a hierarchy of routines linked directly or by magnetic tape.
Large natural water bodies such as the Great Lakes show considerable horizontal variations of chemical and biological parameters. Simulation models of aquatic eco-systems tend to concentrate on vertical variations and often assume that the lake or reservoir is horizontally well-mixed. Such models ca n, in principle, be used to simulate the conditions at various points in a large body of water as a function of local physical parameters and external sources. However, individual zones of a lake are not independent from neighboring ones but are coupled by horizontal transport and dispersion mechanisms. It has been long recognized that upward mixing of nutrients and downward movements of plant material are essential elements of a plankton ecosystem. An investigation of similar processes acting in the horizontal has received less attention and, indeed, would be practically impossible without
As shown in Figure 1, the complete water quality model is visualized as a combination of three major components, namely, the physical model, the biochemical model, and the synthesis model. The physical model consists of a package of modeling routines concerned with the hydrodynamic and thermodynamic properties of large water bodies, whereas the biochemical model simulates the various chemical, biological, and geological processes in the lake. The synthesis model is responsible for the coupling of physical and biological components and coordinates the operation of the total model. In essence, the synthesis model is a transport model linking different compartments of the water body and various natural processes within each of these compartments. Mathematically, the synthesis model is formulated on the basis of the so-called advection-diffusion equation, which expresses the principle of mass conservation as applied to all
495
model variables.
hydrodynamic results are interpolated and integrated to obtain the flow of water perpendicular to any arbitrary boundary separating t he volume elements of the biochemical model. Again this operation is subject to certain restrictions, in particular, the conservation of water mass contained in each element . This is not a trivial problem when some form of interpolation is involved.
Let us consider a volume element within the lake bounded by imaginary surfaces which are permeable with regard to the flow of water and suspended or dissolved material. At any time the total mass of a given chemical substance in this element can be expressed as the product of the mean concentration of this substance and the volume of the element. Now this mass can change by essentially three processes. The first one consists of Without entering into the details of the advection, transports by organized currents and vertical water diffusion, and source calculations, the following motions through the various faces of the volume summarizes the basic principles. The advective element, usually referred to as advection. The flux of material across a permeable interface is second process results from quasi-random fluid computed by multiplication of the flow of water perpendicular to the interface with the appropriate motions which tend to mix the properties of adjacent volume elements in a manner analogous to concentration of the respective material. A molecular diffusion. The third contribution comes typical assumption is that the concentration of from external sources or the internal exchange the volume element from which the outflow takes place, represents a good approxi mation to the between individual chemical substances appearing as a source of one element and a sink of others. concentration at the interface between two The advection-diffusion equation states that the adjacent elements. On the other hand, the time rate of change of mass of a given substance, contributions from diffusive mi xing processes are contained in the element considered, can be simulated on the basis of the gradient hypothesis. computed by accumulating the above effects. Once Thus the diffusive flux through an interface is the time rate of change is known, it is a simple taken to be proportional to the difference in matter to compute the future concentrations from concentrations of adjacent elements and its the present one, at least for a reasonably short direction is from high to low concentration. The time step during which the environmental conditions constant of proportionality, usually called the may be considered quasi-invariant. After sub-grid-scale diffusion coefficient, may be based completing this computation for all the volume on dye-experiments in lakes or oceans. A elements of the model, the advection, diffusion, detailed investigation of this numerical procedure was made by Lam and Simons(l), on the basis of and sources can be calculated for these new conditions and the cycle can be repeated until the one year of physical and chemical data for Lake Erie. First a hydrodynamic model was run desired prediction is completed. In practice, water quality models of this type may divide the to compute the water circulations throughout the year from meteorological records for shore water basin in more than one hundred compartments, and they will use time steps of one day or less stations . The model output was then combined with to predict seasonal and yearly changes in a lake. chemical observations taken at approximately monthly intervals by ship cruises on the lake. Using the principles outlined above, the The synthesis model is coupled with the physical distribution of a quasi-conservative substance, lake model through the advection-diffusion in this case chloride, was predicted from one mechanism and with the biochemical models through cruise to the next and compared with observations. the source terms in the governing equation. In turn, some of the biological models will draw After a satisfactory calibrat~Qn of this model directly upon the physical model for such was obtained, Lam and Jaquet( ) next considered parameters as temperature and radiation which a chemical component interacting with the lake strongly affect biological processes. Since the bottom, thus including a source term at the lower feedback from the biochemical model components interface of the model. Before further elaborating on the physical model can be neglected to a first on these biochemical aspects it is useful to approximation, the most practical interfacing with evaluate the present physical model. the physical model will consist of computer disks or magnetic tapes. Thus the physical input for the PHYSICAL MODEL water quality model is obtained first, partly from observations and partly from hydrodynamic and The purpose of a physical model is to supply water thermodynamic models, and stored for subsequent circulations as well as temperature and radiation use by the synthesis model. Since the required data essential for computing algal growth and resolution in time and space is in general other biological processes. Among the applications different for hydrodynamic versus water quality of temperature information, a very important factor modeling, the first step is to make the necessary is the effect of stratification on vertical mi xing. adjustments. The hydrodynamic models to be In the winter the lakes show a fairly uniform considered here have a typical horizontal vertical temperature profile, but by early summer resolution of 5 km and a ti me step of a few minutes, most lakes exhibit a well-defined thermocline much smaller than the corresponding values for which effectively prohibits mixing between the the biological models mentioned above. Thus the warm upper layers and the cold bottom water. Since output from the physical models is first averaged algal growth is confined to ~he photosynthetic over time. This averaging period is not arbitrary, upper layers, the upper part of the lake tends but as we will see subsequently, it depends on the to get depleted fro m nutrients. When the autumn geographic position of the lake . Next the cooling starts at the lake surface, the lake is
496
thoroughly mixed by sinking of co1der water and rising of warmer water. This seasonal cycle of temperature stratification and associated variations in vertical mixing have been mode1ed with considerable success by a large number of investigators. Similarly, very accurate atmospheric radiation models have been developed and are available in the meteorological literature. In our present studies we have decided to by-pass this step in the mode1ing program, and instead, this model input is obtained directly from observations. The most ambitious program to document the detailed spatial and temporal distribution of physical lake properties was the 1972 International Field Year on Lake Ontario. This data base allows a full description of surface radiation and lake temperatures throughout the year. Furthermore, the vast amount of current meter observations taken in the course of this field program, were used to validate hydrodynamic models. Thus, in conjunction with the field experiments, a comprehensive modeling program was initiated to simulate transports of water and material in the lake.
the second relates to time variation of local currents. One of the most important effects controlling the flow pattern in a lake is the bottom topography . Indeed, when one considers water transports integrated between the surface and the bottom of the lake, such transports tend to run parallel to the bottom contours. Furthermore, the Great Lakes have a typical elongated shape and show a relatively smooth transition in depth from the shore zones to the interior. Now if the wind blows across the lake, in particular, along its length axis, the water in the nearshore area will be carried along faster than the deep water because there is a greater mass of water per unit surface in the interior of the lake. Actually, the deep water may start moving against the wind because the lake is bounded and thus the transport integrated over any cross section must be necessarily small. It turns out that such current reversals are likely to occur along the bottom contour corresponding to the mean depth . In view of this, we choose this particular depth contour as the boundary between shore elements and deep water segments in our water quality model. Figure 3 illustrates this lake response as computed by our hydrodynamical model of Lake Ontario after a strong westerly wind. Since the northern half of the lake is much shallower than the southern half, the water transports are generally eastward along the northern shore.
The hydrodynamic model verificatio~ ba~ been summarized in two papers by Simons(3) (4). The model consists of four layers and has a horizontal gridmesh of 5 km. At each gridpoint and for each layer the model computes currents, vertical velocities, and temperatures from specified initial conditions and for given windstresses acting at the water surface. In view of our earlier remarks, the temperature is continually updated by recourse to observations taken by weekly ship cruises. Although the model appeared reasonably successful in predicting temperature changes from one week to the next, the accumulation of errors would soon result in poor agreement with observations, since the model does not incorporate all mechanisms contributing to seasonal temperature changes. On the other hand, the temperature is an excellent indicator of the model performance with regard to transport calculations, since the redistribution of heat in the lake over a period of a week or so is essentially controlled by water transports just like any other dissolved or suspended material. Figure 2 shows an example of such comparisons. Given two ship cruises spaced two weeks apart, one can derive the observed temperature changes for any given depth as a function of horizontal coordinates, and a similar distribution can be obtained from the model output . It is quite apparent that the model results display the essential features of the observed distribution pattern. Areas of cooling, indicated by dashed lines, can be related to upwelling zones where co1der water from the lake bottom moves to the surface . This process is caused by winds moving the local surface water away from the shore and it is the primary agent responsible for shortterm temperature changes in lakes or oceans .
The time variation of local currents in a large lake is largely controlled by the effect of the earth's rotation. Thus water particles tend to move in a straight line in an absolute frame of reference, but to an observer attached to the rotating earth it appears as a spiraling motion. If one considers the speed and direction of the current vector observed or computed in a given point as a function of time, the result looks like Figure 4. This figure displays observations taken by a string of current meters at four depths and corresponding model results for four layers. The clockwise rotation of the current vector is quite apparent, the period being about 17 hours for Lake Ontario. At the same time it can be seen that the amplitude varies ~ith the same periodicity, thus implying a net displacement over a number of days. It will now be apparent why the averaging period for the hydrodynamical model output is not arbitrary. Obviously, meaningful averaged water transports are obtained by choosing the averaging period equal to this rotational period. The most accurate procedure is to apply a digital filter to the original time series. BIOCHEMICAL MODEL Among the variety of water quality models, the more interesting biochemical routines are the recently developed plankton models of large lakes. The nutrient-plankton interactions in lakes and oceans have been simulated in studies of aquatic ecosystem dynamics, but usually without considering spatial interactions. In the last couple of years, studies of eutrophication and water quality problems have tried to include a full description of plankton interactions as well as a complete physical submodel. In particular,
Among the results of hydrodynamic modeling studies, there are two that may be mentioned here because they are of direct relevance to water quality modeling. The first one concerns the horizontal pattern of water transports in large lakes whereas
497
the Lake Washington model by Chen and Orlob(5) and the Lake Ontario model by Thoman et al. (6) may be mentioned in the context of large lake mana gement models. The biochemical submodel deals with kinetic interactions between chemical and biological system variables and its mathematical representation is characterized by a system of coupled ordinary differential equations. Each equation of the system prescribes the time rate of change of one system variable by accumulating the effects of its interaction with all remaining variables. Central to the aquatic plankton model are the algal or phytoplankton components and three major conversion processes between system variables, namely, photosynthesis, grazing, and regeneration of nutrients. Photosynthesis is the uptake of nutrients by phytoplankton resulting in a growth of algal components at the expense of chemical components. Grazing represents a loss of phytoplankton and results in the growth of zooplankton and higher trophic model components. Regeneration of nutrients is the recycling of chemical components through excretion by zooplankton and the oxidation of dead organic material. The photosynthetic activity and the resulting growth of algal populations depend on light, temperature, and the availability of nutrients. In a comprehensive water quality model the first two parameters are supplied by the physical submodel. Commonly considered nutrients are ammonia and nitrate nitrogen and soluble reactive phosphate. Effects of multiple nutrients are usually incorporated by multiplication of their individual effects. Phytoplankton mortality is not due only to predation but endogeneous respiration and other causes of natural death contribute to plankton losses. Similarly, the growth of zooplankton is not related to grazing in a straightforward manner since the simulated growth of a population must include births, increase of individual weight, and food assimilation. A typical formulation of grazing activity appears as the product of the biomass of the prey and the predator. The regeneration of nutrients constitutes the feedback loop from plankton to nutrient components in the biochemical model. Dead cell material and excreta of indigested food decompose by bacterial activity, resulting in the release of nutrients. The more complex models explicitly consider dead organic components as a link between plankton death and nutrient regeneration. Decomposition of organic material may be represented by a three-step first order reaction from organic nitrogen to ammonia, nitrite, and nitrate, and a first order decay from organic phosphorus to orthophosphate. In turn, a simulation of this oxidation process may lead to the inclusion of the dissolved oxygen balance in the model. It is clear that any operational water quality model should be extremely flexible to accommodate any new developments in the simulation of such complex phenomena.
REFERENCES (1)
Lam, D.C.L. and T.J. Simons, Numerical Computations of Advective and Diffusive Transports of Chloride in Lake Erie for the Year 1970. Canada Centre for Inland Waters, Burlington, Ontario, (1974), 17 pp.
(2)
Lam, D.C.L. and J.M. Jaquet, Computation of Physical Transport and Regeneration of Phosphorus in Lake Erie in the Fall of 1970 . Canada Centre for Inland Waters, Burlington, Ontario, (1974), 20 pp.
(3)
Simons, T.J., Verification of Numerical Models of Lake Ontario, Part 1, Circulation in Spring and Early Summer. J. Phys. Oceanogr., Vol. 4, (1974), pp 507-523.
(4)
Simons, T.J., Verification of Numerical Models of Lake Ontario, Part 2, Stratified Circulations and Temperature Changes. J. Phys. Oceanogr., Vol . 5, (1975).
(5)
Chen, C.W. and G.T. Orlob, Ecologic Simulation for Aquatic Environments. Water Resources Engineers, Inc., Walnut Creek, California, (1972), 156 pp.
(6) Thomann, R. V., D.M. DiToro, R.P. Winfield and D.J. O'Connor, Mathematical Modeling of Phytoplankton in Lake Ontario. Manhattan College, Bronx, New York, (1974),177 pp.
., ATE R
Ph~sical
Q U A liT Y
MOO E l
Synthes i 5 Model
'·lodel
Bi ochemi ca 1 Model
Photosynthes i s
IRadiation
n
Tempera ture
Wate r Transpor t s
Advection +
c::::::;'>
Diffusion
+ Sources
Mixing and
!Dispersion
n
IRiver Discharges
I Figu re 1.
498
Predation
~
Nutrient Regene r a t 10n
Sediment Exchange
Surface Aeration
Final State
Major conponents of cl typical wate r quality model f or large lakes .
LAKE ONTARIO OBSERVED TEMPERATURE CHANGES QC AUGUST 2 - 15. 15, 1972 ; 20 - 40 METRES
2 )
).\ J ,,'
.--//,
, ~. t
/
\
1
,
,
,
•
~8
~~
~
l(~~~\~/ ,"
/'
/-\"~_ _ _~), ,
0
1
,,
.'
I'
-"
"
499
i " '2 .
,
'
-
~
,
(
"''."
... . y
,;
Figure 2.
·6
-
.'....
,-<~-' ~~ ---0"'~'
if
?)
·8
<.
,
r
.
"l .•••
~
,
G'
o
\ ../
~-
r
,:
'1
~ _... '.--.' - ",: ' .'
-
"
.I,
11
· ::-:'~ ..... : (·::;"·:'t. ;,rtfr ,
"
' •.0; , .• 'd'
';" 1',1
3 ;: ' >.:
' " . _.r·· .
:t0~-' ~ <--1 ', "\ ,..
, .'1, .
' (~, ~.
J/ ~o ~
>
'! :!
,' 1 { .
COMPUTED TEMPERATURE CHANGES QC AUGUST 2 - 15, 15. 1972 ; 20 - 40 METRES
,, :
'-61
( /:
)
"
it ','~ (I
'<-'
~/
"~ (~~ 7~~ c~~~- j( "" , " \, ' ~ ~ VI , -' - \ ( \ , "'r{' ..".:..
LAKE ONTARIO
,
)
~
..
."'.'~
•
~
,__ _ _
(.
, ,:, J'(/~ ( \ . ,.' (
'/ '.·" ,,;.,
.
~-- . ~~~ ~ .; . . ':.Jj )
2)
'r, '
"
'". -3 '
'3 '
- "~'"
:
~
'
.. : l
) ':
- -,'
:
.. .. o.: :.'. ·5 r·,
0
.........~ ... \" ,.' ~ '"
~,~ ~~
,
I
'
i
,
Horizontal pattern of temperature changes obtained from two ship cruises on Lake Ontario, spaced two weeks apart, and corresponding changes computed by a three-dimensional hydrodynamic model. model . Negative temperature changes, indicated by dashed lines, point at upwelling areas where colder water from the lake bottom moves to the surface surface..
• r
LAKE ONTARIO
... ... . ".
COMPUTED WATER TRANSPORTS AUGUST 10 , 1972 ; STRAT IFIED
•
~ ~
~~
•
I
r,
I '
(
..
:l ",1
..
'
;
.... .. '
.... .. . ~ ~
,
'
,
.
-- : :
~
--
.
. . "L'to' '0
... ..
10 M 2/S
500
LAKE ONTARIO COMPUTED WATER TRANSPORTS 12 , 1972 ; STRAT IFIED AUGUST 12,1972;
.. :,J
0' , •
•
•
•
•
•
,
/
•
I
~
\1~
.'
Figure 3.
"
.
Water transports, integrated between lake surface and bottom, as computed by a four-layer hydrodynamic model \vesterly wind and two days later. The immediate effect of the storm is that the iimmediately mmediately after a strong westerly shallow water along the northern shore moves in the direction of the wind while the deep water returns against the wind.
".
30
'0
'.
·
';; u
-.."'...Eu 0 W
'"'"w
lL CL
'0
20 2+
10 10
+
---Tr\---
~
~ 30
"
m
11'1
(\_ -; (
! \
'
I
.
'
11'1 1ft .
__
Obl.r v. d Obs.rv.d
,0",
\\
j'
\\
/\\.
~
/
le
I~
_ . ... 50.IL --', :SO M .
-'-"--T-----r-.- ..r~ r ~ ____ _ _ _ _ _ _-j~---=_----+--------+_--------___!
,/ ,.
I
I
_ _ _" ,
--- -- -.
Ob •• rvld Obl.r wld
,""', ( . \
m
..-- 10 M.
\
\\
/t~ \
(f) en
N
z
0o
5
~ t-
u U
'"'"a:ow
--j f--------f~__+--_il_----_t_f___j__
- .. - - -
0
.!_2_ - - l __ _ _ _ 3: f-L--J-UN-E--7-2_ -.-_-_-~4--.--2-3-_-JU-!-.E-._-'..L-7-2~~~~_-t:~~~~242~4-'--JJUNE U-N-E-7-2---I'--.-=25==J=U=N=E=:~=====!====--'-L-=~=~~=====1~~:=~2-=7==JU=N=E:_- .:7:-:2:-_---L.--'-----j N
JUNE
72
27
JUNE
72
30
501
I:.I
...u
·.
-!: -..
20 - - .
E
I\\. "..
u
~ 0o
'"'""-ww CL
Cl> III
10 10
in'"\" " "Cc", "
",
-
--
\ " \/ "
-
/ ........ /"
I- -
-----i - - - - - - - - + - - -I-.f-- -.-----.-.-----. --.-.--I
/ ' \\-yt'--!
I
Computed Comput.d
Computed
.- - ~: :-:. .'
,~,. r'' -. '. J( -·---- +---~-I --L.i F'''<~x-,c~r) iI , ------.J ,J -~ --- . I ,~ J ' / tI ",'- 'J I .!
f-,- ---
.:' ' " , ,____
___________
/-C,
--I
""
, ' -",
;
,~
" , Jk
,
, _, " I
'\,
"
_.--,0/ .~ .../ c ~
', '
N
~~. / ".\,- -.- ---
D---2 / ' ,,,/
w
z o
... ...a:
5s
~
~,~:r'''+:-=-'~--:-:/
--- +
t-
'"o
---;r--------'-
//
u U
E
- - _ . _ - .,L
I
--1 --- --
/
N NL-~-I
___________L______~LL / " / ____~____~ /1 / __________JI@ L----------~~-L------~------~--~~---------It I hi
Figure 4. 4. Time series of speed and direction of observed currents for a vertical string of current meters in Lake Ontario Onta r io co rres ponding results from from a four-layer four-l ayer hydrodynamic model mode l of Lake Ontario. Ontario . Due to a storm on 22 June 1972, and corresponding the currents first increase in magnitude magnit ude and subsequently display di sp lay a spiraling motion motion related to the earth's ea r th's rotation.
T.J. Simons Canada Centre for Inland Waters P.O. Box 5050 Burlington, Ontario L7R 4A6
ABSTRACT Comprehensive models for water quality control and management of large lakes have two identifiable components, namely a physical submodel and a package of biochemical routines. The physical model component supplies all essential information on water transports in the basin and describes the distribution of temperature and radiation required for a simulation of photosynthesis and other biological processes. The biochemical submodel is concerned with the various chemical and biological interactions resulting in a local increase of one biochemical variable at the expense of others. The two major components of the water quality model are tied together by a computational framework based essentially on a generalized advection-diffusion equation which expresses the principle of mass conservation as applied to each model variable. An essential property of such interdisciplinary model is a modular structure allowing for modifications pertaining to specific hydrodynamic or biological processes without affecting the overall computer program . This paper discusses the different components of current water quality models of the North American Great Lakes with particular emphasis on the transport mechanisms acting in such lakes .
the use of modern-day high-speed computers. Prerequisites for such investigations are obviously an adequate physical lake model, a realistic biochemical model, and a practical transport framework to couple these submodels. This paper will discuss some aspects of this interfacing problem followed by a brief description of the pertinent submodels . MODEL FRAMEWORK
I NTRODUCTI ON
The purpose of this effort is to construct a computational framework to serve as a basis of interaction and synthesis for an interdisciplinary group of scientists involved with water quality models. The underlying mathematical principles are that the models are prognostic in character, such that all variables are computed forward in time from specified initial values, and all changes of model variables are computed from mass conservation principles, regardless of the type of natural process to be simulated. From a computer-technical viewpoint, the conditions are that scientific concepts regarding specific hydrodynamic or biological processes can be incorporated with no more than a superficial knowledge of submodels pertaining to other disciplines. Thus the model is constructed in a modular fashion, characterized by a hierarchy of routines linked directly or by magnetic tape.
Large natural water bodies such as the Great Lakes show considerable horizontal variations of chemical and biological parameters. Simulation models of aquatic eco-systems tend to concentrate on vertical variations and often assume that the lake or reservoir is horizontally well-mixed. Such models ca n, in principle, be used to simulate the conditions at various points in a large body of water as a function of local physical parameters and external sources. However, individual zones of a lake are not independent from neighboring ones but are coupled by horizontal transport and dispersion mechanisms. It has been long recognized that upward mixing of nutrients and downward movements of plant material are essential elements of a plankton ecosystem. An investigation of similar processes acting in the horizontal has received less attention and, indeed, would be practically impossible without
As shown in Figure 1, the complete water quality model is visualized as a combination of three major components, namely, the physical model, the biochemical model, and the synthesis model. The physical model consists of a package of modeling routines concerned with the hydrodynamic and thermodynamic properties of large water bodies, whereas the biochemical model simulates the various chemical, biological, and geological processes in the lake. The synthesis model is responsible for the coupling of physical and biological components and coordinates the operation of the total model. In essence, the synthesis model is a transport model linking different compartments of the water body and various natural processes within each of these compartments. Mathematically, the synthesis model is formulated on the basis of the so-called advection-diffusion equation, which expresses the principle of mass conservation as applied to all
495
model variables.
hydrodynamic results are interpolated and integrated to obtain the flow of water perpendicular to any arbitrary boundary separating t he volume elements of the biochemical model. Again this operation is subject to certain restrictions, in particular, the conservation of water mass contained in each element . This is not a trivial problem when some form of interpolation is involved.
Let us consider a volume element within the lake bounded by imaginary surfaces which are permeable with regard to the flow of water and suspended or dissolved material. At any time the total mass of a given chemical substance in this element can be expressed as the product of the mean concentration of this substance and the volume of the element. Now this mass can change by essentially three processes. The first one consists of Without entering into the details of the advection, transports by organized currents and vertical water diffusion, and source calculations, the following motions through the various faces of the volume summarizes the basic principles. The advective element, usually referred to as advection. The flux of material across a permeable interface is second process results from quasi-random fluid computed by multiplication of the flow of water perpendicular to the interface with the appropriate motions which tend to mix the properties of adjacent volume elements in a manner analogous to concentration of the respective material. A molecular diffusion. The third contribution comes typical assumption is that the concentration of from external sources or the internal exchange the volume element from which the outflow takes place, represents a good approxi mation to the between individual chemical substances appearing as a source of one element and a sink of others. concentration at the interface between two The advection-diffusion equation states that the adjacent elements. On the other hand, the time rate of change of mass of a given substance, contributions from diffusive mi xing processes are contained in the element considered, can be simulated on the basis of the gradient hypothesis. computed by accumulating the above effects. Once Thus the diffusive flux through an interface is the time rate of change is known, it is a simple taken to be proportional to the difference in matter to compute the future concentrations from concentrations of adjacent elements and its the present one, at least for a reasonably short direction is from high to low concentration. The time step during which the environmental conditions constant of proportionality, usually called the may be considered quasi-invariant. After sub-grid-scale diffusion coefficient, may be based completing this computation for all the volume on dye-experiments in lakes or oceans. A elements of the model, the advection, diffusion, detailed investigation of this numerical procedure was made by Lam and Simons(l), on the basis of and sources can be calculated for these new conditions and the cycle can be repeated until the one year of physical and chemical data for Lake Erie. First a hydrodynamic model was run desired prediction is completed. In practice, water quality models of this type may divide the to compute the water circulations throughout the year from meteorological records for shore water basin in more than one hundred compartments, and they will use time steps of one day or less stations . The model output was then combined with to predict seasonal and yearly changes in a lake. chemical observations taken at approximately monthly intervals by ship cruises on the lake. Using the principles outlined above, the The synthesis model is coupled with the physical distribution of a quasi-conservative substance, lake model through the advection-diffusion in this case chloride, was predicted from one mechanism and with the biochemical models through cruise to the next and compared with observations. the source terms in the governing equation. In turn, some of the biological models will draw After a satisfactory calibrat~Qn of this model directly upon the physical model for such was obtained, Lam and Jaquet( ) next considered parameters as temperature and radiation which a chemical component interacting with the lake strongly affect biological processes. Since the bottom, thus including a source term at the lower feedback from the biochemical model components interface of the model. Before further elaborating on the physical model can be neglected to a first on these biochemical aspects it is useful to approximation, the most practical interfacing with evaluate the present physical model. the physical model will consist of computer disks or magnetic tapes. Thus the physical input for the PHYSICAL MODEL water quality model is obtained first, partly from observations and partly from hydrodynamic and The purpose of a physical model is to supply water thermodynamic models, and stored for subsequent circulations as well as temperature and radiation use by the synthesis model. Since the required data essential for computing algal growth and resolution in time and space is in general other biological processes. Among the applications different for hydrodynamic versus water quality of temperature information, a very important factor modeling, the first step is to make the necessary is the effect of stratification on vertical mi xing. adjustments. The hydrodynamic models to be In the winter the lakes show a fairly uniform considered here have a typical horizontal vertical temperature profile, but by early summer resolution of 5 km and a ti me step of a few minutes, most lakes exhibit a well-defined thermocline much smaller than the corresponding values for which effectively prohibits mixing between the the biological models mentioned above. Thus the warm upper layers and the cold bottom water. Since output from the physical models is first averaged algal growth is confined to ~he photosynthetic over time. This averaging period is not arbitrary, upper layers, the upper part of the lake tends but as we will see subsequently, it depends on the to get depleted fro m nutrients. When the autumn geographic position of the lake . Next the cooling starts at the lake surface, the lake is
496
thoroughly mixed by sinking of co1der water and rising of warmer water. This seasonal cycle of temperature stratification and associated variations in vertical mixing have been mode1ed with considerable success by a large number of investigators. Similarly, very accurate atmospheric radiation models have been developed and are available in the meteorological literature. In our present studies we have decided to by-pass this step in the mode1ing program, and instead, this model input is obtained directly from observations. The most ambitious program to document the detailed spatial and temporal distribution of physical lake properties was the 1972 International Field Year on Lake Ontario. This data base allows a full description of surface radiation and lake temperatures throughout the year. Furthermore, the vast amount of current meter observations taken in the course of this field program, were used to validate hydrodynamic models. Thus, in conjunction with the field experiments, a comprehensive modeling program was initiated to simulate transports of water and material in the lake.
the second relates to time variation of local currents. One of the most important effects controlling the flow pattern in a lake is the bottom topography . Indeed, when one considers water transports integrated between the surface and the bottom of the lake, such transports tend to run parallel to the bottom contours. Furthermore, the Great Lakes have a typical elongated shape and show a relatively smooth transition in depth from the shore zones to the interior. Now if the wind blows across the lake, in particular, along its length axis, the water in the nearshore area will be carried along faster than the deep water because there is a greater mass of water per unit surface in the interior of the lake. Actually, the deep water may start moving against the wind because the lake is bounded and thus the transport integrated over any cross section must be necessarily small. It turns out that such current reversals are likely to occur along the bottom contour corresponding to the mean depth . In view of this, we choose this particular depth contour as the boundary between shore elements and deep water segments in our water quality model. Figure 3 illustrates this lake response as computed by our hydrodynamical model of Lake Ontario after a strong westerly wind. Since the northern half of the lake is much shallower than the southern half, the water transports are generally eastward along the northern shore.
The hydrodynamic model verificatio~ ba~ been summarized in two papers by Simons(3) (4). The model consists of four layers and has a horizontal gridmesh of 5 km. At each gridpoint and for each layer the model computes currents, vertical velocities, and temperatures from specified initial conditions and for given windstresses acting at the water surface. In view of our earlier remarks, the temperature is continually updated by recourse to observations taken by weekly ship cruises. Although the model appeared reasonably successful in predicting temperature changes from one week to the next, the accumulation of errors would soon result in poor agreement with observations, since the model does not incorporate all mechanisms contributing to seasonal temperature changes. On the other hand, the temperature is an excellent indicator of the model performance with regard to transport calculations, since the redistribution of heat in the lake over a period of a week or so is essentially controlled by water transports just like any other dissolved or suspended material. Figure 2 shows an example of such comparisons. Given two ship cruises spaced two weeks apart, one can derive the observed temperature changes for any given depth as a function of horizontal coordinates, and a similar distribution can be obtained from the model output . It is quite apparent that the model results display the essential features of the observed distribution pattern. Areas of cooling, indicated by dashed lines, can be related to upwelling zones where co1der water from the lake bottom moves to the surface . This process is caused by winds moving the local surface water away from the shore and it is the primary agent responsible for shortterm temperature changes in lakes or oceans .
The time variation of local currents in a large lake is largely controlled by the effect of the earth's rotation. Thus water particles tend to move in a straight line in an absolute frame of reference, but to an observer attached to the rotating earth it appears as a spiraling motion. If one considers the speed and direction of the current vector observed or computed in a given point as a function of time, the result looks like Figure 4. This figure displays observations taken by a string of current meters at four depths and corresponding model results for four layers. The clockwise rotation of the current vector is quite apparent, the period being about 17 hours for Lake Ontario. At the same time it can be seen that the amplitude varies ~ith the same periodicity, thus implying a net displacement over a number of days. It will now be apparent why the averaging period for the hydrodynamical model output is not arbitrary. Obviously, meaningful averaged water transports are obtained by choosing the averaging period equal to this rotational period. The most accurate procedure is to apply a digital filter to the original time series. BIOCHEMICAL MODEL Among the variety of water quality models, the more interesting biochemical routines are the recently developed plankton models of large lakes. The nutrient-plankton interactions in lakes and oceans have been simulated in studies of aquatic ecosystem dynamics, but usually without considering spatial interactions. In the last couple of years, studies of eutrophication and water quality problems have tried to include a full description of plankton interactions as well as a complete physical submodel. In particular,
Among the results of hydrodynamic modeling studies, there are two that may be mentioned here because they are of direct relevance to water quality modeling. The first one concerns the horizontal pattern of water transports in large lakes whereas
497
the Lake Washington model by Chen and Orlob(5) and the Lake Ontario model by Thoman et al. (6) may be mentioned in the context of large lake mana gement models. The biochemical submodel deals with kinetic interactions between chemical and biological system variables and its mathematical representation is characterized by a system of coupled ordinary differential equations. Each equation of the system prescribes the time rate of change of one system variable by accumulating the effects of its interaction with all remaining variables. Central to the aquatic plankton model are the algal or phytoplankton components and three major conversion processes between system variables, namely, photosynthesis, grazing, and regeneration of nutrients. Photosynthesis is the uptake of nutrients by phytoplankton resulting in a growth of algal components at the expense of chemical components. Grazing represents a loss of phytoplankton and results in the growth of zooplankton and higher trophic model components. Regeneration of nutrients is the recycling of chemical components through excretion by zooplankton and the oxidation of dead organic material. The photosynthetic activity and the resulting growth of algal populations depend on light, temperature, and the availability of nutrients. In a comprehensive water quality model the first two parameters are supplied by the physical submodel. Commonly considered nutrients are ammonia and nitrate nitrogen and soluble reactive phosphate. Effects of multiple nutrients are usually incorporated by multiplication of their individual effects. Phytoplankton mortality is not due only to predation but endogeneous respiration and other causes of natural death contribute to plankton losses. Similarly, the growth of zooplankton is not related to grazing in a straightforward manner since the simulated growth of a population must include births, increase of individual weight, and food assimilation. A typical formulation of grazing activity appears as the product of the biomass of the prey and the predator. The regeneration of nutrients constitutes the feedback loop from plankton to nutrient components in the biochemical model. Dead cell material and excreta of indigested food decompose by bacterial activity, resulting in the release of nutrients. The more complex models explicitly consider dead organic components as a link between plankton death and nutrient regeneration. Decomposition of organic material may be represented by a three-step first order reaction from organic nitrogen to ammonia, nitrite, and nitrate, and a first order decay from organic phosphorus to orthophosphate. In turn, a simulation of this oxidation process may lead to the inclusion of the dissolved oxygen balance in the model. It is clear that any operational water quality model should be extremely flexible to accommodate any new developments in the simulation of such complex phenomena.
REFERENCES (1)
Lam, D.C.L. and T.J. Simons, Numerical Computations of Advective and Diffusive Transports of Chloride in Lake Erie for the Year 1970. Canada Centre for Inland Waters, Burlington, Ontario, (1974), 17 pp.
(2)
Lam, D.C.L. and J.M. Jaquet, Computation of Physical Transport and Regeneration of Phosphorus in Lake Erie in the Fall of 1970 . Canada Centre for Inland Waters, Burlington, Ontario, (1974), 20 pp.
(3)
Simons, T.J., Verification of Numerical Models of Lake Ontario, Part 1, Circulation in Spring and Early Summer. J. Phys. Oceanogr., Vol. 4, (1974), pp 507-523.
(4)
Simons, T.J., Verification of Numerical Models of Lake Ontario, Part 2, Stratified Circulations and Temperature Changes. J. Phys. Oceanogr., Vol . 5, (1975).
(5)
Chen, C.W. and G.T. Orlob, Ecologic Simulation for Aquatic Environments. Water Resources Engineers, Inc., Walnut Creek, California, (1972), 156 pp.
(6) Thomann, R. V., D.M. DiToro, R.P. Winfield and D.J. O'Connor, Mathematical Modeling of Phytoplankton in Lake Ontario. Manhattan College, Bronx, New York, (1974),177 pp.
., ATE R
Ph~sical
Q U A liT Y
MOO E l
Synthes i 5 Model
'·lodel
Bi ochemi ca 1 Model
Photosynthes i s
IRadiation
n
Tempera ture
Wate r Transpor t s
Advection +
c::::::;'>
Diffusion
+ Sources
Mixing and
!Dispersion
n
IRiver Discharges
I Figu re 1.
498
Predation
~
Nutrient Regene r a t 10n
Sediment Exchange
Surface Aeration
Final State
Major conponents of cl typical wate r quality model f or large lakes .
LAKE ONTARIO OBSERVED TEMPERATURE CHANGES QC AUGUST 2 - 15. 15, 1972 ; 20 - 40 METRES
2 )
).\ J ,,'
.--//,
, ~. t
/
\
1
,
,
,
•
~8
~~
~
l(~~~\~/ ,"
/'
/-\"~_ _ _~), ,
0
1
,,
.'
I'
-"
"
499
i " '2 .
,
'
-
~
,
(
"''."
... . y
,;
Figure 2.
·6
-
.'....
,-<~-' ~~ ---0"'~'
if
?)
·8
<.
,
r
.
"l .•••
~
,
G'
o
\ ../
~-
r
,:
'1
~ _... '.--.' - ",: ' .'
-
"
.I,
11
· ::-:'~ ..... : (·::;"·:'t. ;,rtfr ,
"
' •.0; , .• 'd'
';" 1',1
3 ;: ' >.:
' " . _.r·· .
:t0~-' ~ <--1 ', "\ ,..
, .'1, .
' (~, ~.
J/ ~o ~
>
'! :!
,' 1 { .
COMPUTED TEMPERATURE CHANGES QC AUGUST 2 - 15, 15. 1972 ; 20 - 40 METRES
,, :
'-61
( /:
)
"
it ','~ (I
'<-'
~/
"~ (~~ 7~~ c~~~- j( "" , " \, ' ~ ~ VI , -' - \ ( \ , "'r{' ..".:..
LAKE ONTARIO
,
)
~
..
."'.'~
•
~
,__ _ _
(.
, ,:, J'(/~ ( \ . ,.' (
'/ '.·" ,,;.,
.
~-- . ~~~ ~ .; . . ':.Jj )
2)
'r, '
"
'". -3 '
'3 '
- "~'"
:
~
'
.. : l
) ':
- -,'
:
.. .. o.: :.'. ·5 r·,
0
.........~ ... \" ,.' ~ '"
~,~ ~~
,
I
'
i
,
Horizontal pattern of temperature changes obtained from two ship cruises on Lake Ontario, spaced two weeks apart, and corresponding changes computed by a three-dimensional hydrodynamic model. model . Negative temperature changes, indicated by dashed lines, point at upwelling areas where colder water from the lake bottom moves to the surface surface..
• r
LAKE ONTARIO
... ... . ".
COMPUTED WATER TRANSPORTS AUGUST 10 , 1972 ; STRAT IFIED
•
~ ~
~~
•
I
r,
I '
(
..
:l ",1
..
'
;
.... .. '
.... .. . ~ ~
,
'
,
.
-- : :
~
--
.
. . "L'to' '0
... ..
10 M 2/S
500
LAKE ONTARIO COMPUTED WATER TRANSPORTS 12 , 1972 ; STRAT IFIED AUGUST 12,1972;
.. :,J
0' , •
•
•
•
•
•
,
/
•
I
~
\1~
.'
Figure 3.
"
.
Water transports, integrated between lake surface and bottom, as computed by a four-layer hydrodynamic model \vesterly wind and two days later. The immediate effect of the storm is that the iimmediately mmediately after a strong westerly shallow water along the northern shore moves in the direction of the wind while the deep water returns against the wind.
".
30
'0
'.
·
';; u
-.."'...Eu 0 W
'"'"w
lL CL
'0
20 2+
10 10
+
---Tr\---
~
~ 30
"
m
11'1
(\_ -; (
! \
'
I
.
'
11'1 1ft .
__
Obl.r v. d Obs.rv.d
,0",
\\
j'
\\
/\\.
~
/
le
I~
_ . ... 50.IL --', :SO M .
-'-"--T-----r-.- ..r~ r ~ ____ _ _ _ _ _ _-j~---=_----+--------+_--------___!
,/ ,.
I
I
_ _ _" ,
--- -- -.
Ob •• rvld Obl.r wld
,""', ( . \
m
..-- 10 M.
\
\\
/t~ \
(f) en
N
z
0o
5
~ t-
u U
'"'"a:ow
--j f--------f~__+--_il_----_t_f___j__
- .. - - -
0
.!_2_ - - l __ _ _ _ 3: f-L--J-UN-E--7-2_ -.-_-_-~4--.--2-3-_-JU-!-.E-._-'..L-7-2~~~~_-t:~~~~242~4-'--JJUNE U-N-E-7-2---I'--.-=25==J=U=N=E=:~=====!====--'-L-=~=~~=====1~~:=~2-=7==JU=N=E:_- .:7:-:2:-_---L.--'-----j N
JUNE
72
27
JUNE
72
30
501
I:.I
...u
·.
-!: -..
20 - - .
E
I\\. "..
u
~ 0o
'"'""-ww CL
Cl> III
10 10
in'"\" " "Cc", "
",
-
--
\ " \/ "
-
/ ........ /"
I- -
-----i - - - - - - - - + - - -I-.f-- -.-----.-.-----. --.-.--I
/ ' \\-yt'--!
I
Computed Comput.d
Computed
.- - ~: :-:. .'
,~,. r'' -. '. J( -·---- +---~-I --L.i F'''<~x-,c~r) iI , ------.J ,J -~ --- . I ,~ J ' / tI ",'- 'J I .!
f-,- ---
.:' ' " , ,____
___________
/-C,
--I
""
, ' -",
;
,~
" , Jk
,
, _, " I
'\,
"
_.--,0/ .~ .../ c ~
', '
N
~~. / ".\,- -.- ---
D---2 / ' ,,,/
w
z o
... ...a:
5s
~
~,~:r'''+:-=-'~--:-:/
--- +
t-
'"o
---;r--------'-
//
u U
E
- - _ . _ - .,L
I
--1 --- --
/
N NL-~-I
___________L______~LL / " / ____~____~ /1 / __________JI@ L----------~~-L------~------~--~~---------It I hi
Figure 4. 4. Time series of speed and direction of observed currents for a vertical string of current meters in Lake Ontario Onta r io co rres ponding results from from a four-layer four-l ayer hydrodynamic model mode l of Lake Ontario. Ontario . Due to a storm on 22 June 1972, and corresponding the currents first increase in magnitude magnit ude and subsequently display di sp lay a spiraling motion motion related to the earth's ea r th's rotation.