Key Words: computational fluid dynamics, phytoplankton spatio-temporal distributions, velocity field, irradiance
R.D. Hedger 1(rdh@geo.ed.ac.uk), N.R.B. Olsen 2, D.G. George 3, P.M. Atkinson 4, T.J. Malthus 1,
1 Department of Geography, University of Edinburgh, Edinburgh, UK, EH8 9XP
2 Department of Hydraulic and Environmental Engineering, Norwegian University of Science and Technology, N-7034 Trondheim, Norway
3 Institute of Freshwater Ecology, Far Sawrey, Ambleside, Cumbria, UK, LA22 OLP
4 Department of Geography, University of Southampton, Southampton, UK, SO17 1BJ
The relationships between the spatio-temporal distribution of phytoplankton concentration and the environmental conditions of Esthwaite Water (a small eutrophic lake in the English Lake District, U.K.) were examined using a 3-D computational fluid dynamics (CFD) model. The water velocity field was obtained through solving the 3-D Navier Stokes equation for turbulent flow on a finite-volume, unstructured non-orthogonal grid. The spatio-temporal distributions of two types of phytoplankton were modelled: the cyanobacterium Microcystis, and the dinoflagellate Ceratium. Cyanobacterial buoyancy were estimated according to the Kromkamp and Walsby model, and dinoflagellate motility was estimated according to a model that we devised using empirical data from Esthwaite Water and other similar lakes. Circulation patterns of water and phytoplankton, as simulated by the CFD model, were similar to those obtained through field observations. Downwind surface drift currents were initiated by wind stress, with sub-surface return gradient currents initiated near the thermocline. Near-surface accumulations of cyanobacteria were pushed downwind by the surface currents and accumulated at downwelling areas, and near-thermocline accumulations of dinoflagellates were pushed upwind by the sub-surface return currents, and accumulated at upwelling areas. In all cases, the Coriolis force greatly influenced patterns, causing a clockwise deflection of water flow and phytoplankton accumulation. Through the use of the CFD model, it was possible to conclude that the horizontal and vertical phytoplankton distributions resulted from the interaction between the vertical motility of the phytoplankton (dependent on the light environment) and the velocity vectors at the depths at which the phytoplankton accumulated (dependent upon wind stress and basin morphometry).
Phytoplankton concentration is spatially heterogenous within many lakes. This heterogeneity is often exhibited as a distinct pattern in the horizontal and vertical domains ('pattern' being defined as non-random spatial variation). This pattern is caused by the interaction between the phenomenological characteristics of phytoplankton species (such as buoyancy) and the velocity field (advection, turbulence) and light field properties of the lake. For instance, at basin-wide scales, horizontal pattern may be caused by positively buoyant species accumulating in downwelling areas at the downwind shore of a lake, and negatively buoyant species accumulating in upwelling areas at the upwind shore of a lake (Heaney, 1976; George and Edwards, 1976; George and Heaney, 1978; George, 1993). Likewise, vertical pattern may be caused by a species moving to a depth of optimum light intensity within the water column, either through its own motility (Heaney and Talling, 1980) or alterations in its buoyancy (Reynolds, 1984).
The aim of this study is to determine if current knowledge of how lakes function can be used in dynamic modelling to simulate phenomena that have been measured. This will involve analysis of the effect of the interaction between specific characteristics of individual species and lake environmental properties, particularly irradiance and current vectors.
Esthwaite Water is a small lake in the English Lake District, approximately 2.5 km long, with a mean depth of 6.4 m and a maximum depth of approximately 15 m.
The lake was chosen as a study area because it has suitable morphometric and phytoplankton phenomenological characteristics for interactions to result in pattern, and it is highly sampled (because of its near proximity to the Institute of Freshwater Ecology, Far Sawrey).
The lake has relatively shallow bays along its east and west sides (Figure 1), which often result in great spatial variation in flow properties across the lake (occasionally, water in the shallows may be near-stationary compared to water in the central trough). Stratification occurs from May until early October, with a thermocline occurring at approximately 6 m. Subsurface return flows occur just above this thermocline (George and Heaney, 1978).
Figure 1. Bathymetry of Esthwaite Water.
The lake is eutrophic, with the dominant phytoplankton being cyanobacteria (including Microcystis, Aphanizomenon and Cryptomonas) and dinoflagellates (including Ceratium hirundinella). Temporal variation in phytoplankton concentration follows a pattern of exponential growth in May through to July, relatively constant levels in August, and rapid decline in September and October (Heaney and Talling, 1980). Since net growth of the population is slow it can be inferred that basin-wide variation is the result of the interaction between phytoplankton characteristics and the velocity field rather than nutrient gradients (Heaney, 1976).
In a study of the spatial distributions of phytoplankton in 1973 (George and Heaney,1978) it was found that spatial distributions were species-dependent. On the 2 August 1973, there was an upwind accumulation of the dinoflagellate Ceratium. It was proposed that the dinoflagellates had avoided the surface, had been entrained in a sub-surface return current, and been upwelled upwind. On the 9 August 1973, there was a downwind accumulation of the the cyanobacterium Microcystis (9 August 1973). It was proposed that the cyanobacteria had stayed at the surface (being positively buoyant), had been entrained in a downwind surface drift current, and had accumulated over downwelling areas downwind. Enough data on existent lake environmental properites exist for establishing the boundary conditions for running a CFD model. It is the primary objective of this paper to simulate these two distributions using CFD modelling.
The SSIIM2.0 (Sediment Simulation In Intakes with Multiblock option) CFD model (Olsen, 1998) was used to estimate the velocity field properties by solving the Navier-Stokes equations for turbulent flow in a general three-dimensional geometry to obtain the water velocity:
( 1 )
where U is the Reynolds-averaged velocity (convection), u is the fluctuating velocity (turbulence), x is a space variable, P is pressure, rho is water density, deltaij is the Kronecker delta term (which is 1 if i = j and 0 if i <> j), and t is time. Terms on the left hand side of the equation are the transient, convective and Coriolis terms. Terms on the right hand side of the equation are the convective and Reynold's stress term.
The Reynold's stress term was approximated by the Boussinesq approximation:
( 2 )
where vT is the eddy viscosity (modelled by the k-epsilon model):
( 3 )
where cmicro is a constant (0.09), k is the turbulent kinetic energy, and epsilon is the rate of dissipation of k.
The equations were discretised with a control-volume approach. The Semi-Implict Method for Pressure Linked Equations (SIMPLE) was used to find the pressure, and the water surface elevation. The second-order upwind scheme was used in the discretisation of the convective terms. The numerical methods are further described by Patankar (1980), Melaaen (1992) and Olsen (1997). Wind-induced shear stress on the lake surface was estimated as a function of wind velocity:
( 4 )
where c10, is the friction coefficient (1.1x10-3, after Bengtsson, 1973), rhoair is the density of air (1.2 kg m-3), and Uair is the wind velocity. The shear stress from the wind was added to the Navier-Stokes equations.
The convection-diffusion equation for algal concentration, c, was solved as:
( 5 )
The source term was a combination of sinking velocity of the diatoms and diatom production due to light.
1. Cyanobacterial buoyancy
Rising velocity was estimated from a modification of Stoke's law:
( 6 )
where v is the vertical velocity (positive if the phytoplankton is rising, negative if it is sinking), g is the gravitational acceleration, r is the colony radius, rho c is the phytoplankton density, rho is the water, A is the proportion of the cell volume relative to the colony volume, phi is the form resistance, and n is the water viscosity. Density was estimated as a function of cell carbohydrate concentration according to the phenomenological method of Kromkamp and Walsby (1990) (see also, Howard et al., 1995, Belov and Giles, 1997):
( 7 )
where rho c
( 8 )
where Io is the subsurface irradiance, and epsilon is the vertical extinction coefficient.
2. Dinoflagellate motility
Since dinoflagellate swimming velocities are typically two orders of magnitude less than those of horizontal water currents, the effect of swimming in the horizontal dimension was negligible so only swimming velocities in the vertical dimension were estimated. Dinoflagellates were parameterised to swim to optimum irradiance levels according to the following function:
( 9 )
where v is the vertical velocity of the dinoflagellates (positive if the dinoflagellates are beneath the depth of the optimum light levels, negative if the dinoflagellates are a lesser depth than the optimum light levels), vmax is the maximum velocity, Iopt is the optimum irradiance (150 micromolar photons m-2 s-1, after Heaney and Furnass, 1980), Iz is the irradiance at the depth of the dinoflagellates, and Iran is an irradiance range parameter.
Wind velocity data were obtained from a weather station 5 km to the south-south-east of Esthwaite Water on the east side of Lake Windermere at Haws Wood (B. N. G. 338500 490500). On 2 August 1973, wind speed varied between 2.5 m s-1 and 5 m s-1 (a mean of 2.75 m s-1) with the dominant direction being from the south-south-west (heading towards 25o). On 9 August 1973, wind speed varied between 2.5 m s-1 and 7 m s-1 (a mean of 3.5 m s-1), again with the dominant direction being from the south-south-west. These wind directions are typical for Esthwaite Water, resulting from the channelling effects of the Windermere and Esthwaite valleys.
Data on irradiance were unavailable so this was estimated according to the method outlined by Kirk (1996). This estimate was decreased according to prevalent cloud conditions (measured at the Institute of Freshwater Ecology, Far Sawrey (O.S. B.N.G. 339000 495500)). On the 2 August 1973, the sky was totally overcast, resulting in a low ambient light environment; on 9 August 1973, there was only one octar of cloud, resulting in a high ambient light environment. The lake was strongly stratified throughout this time (18o C at the surface, 1o C at 12 m depth), with a thermocline at 5-6 m.
No data were available on nitrogen and phosphorous concentrations. However, given that the dominant cause of spatial variation in phytoplankton concentration in Esthwaite Water is wind-redistribution rather than nutrient gradients, it is doubtful whether nutrient inputs and outputs needed to be included. For the short amount of time that was to be simulated (12 hrs) it was assumed that there was no great change in mean phytoplankton concentration (because a steady-state had been achieved for the Summer), so nutrient concentration was set to limit growth.
The two dominant types of phytoplankton that were identified from cell counts of point samples acquired on 31 July 1973 and 7 August 1973 were, respectively, dinoflagellates, (predominantly Ceratium hirundinella) and cyanobacteria (predominantly Microcystis). There were approximately 15 species of phytoplankton identified in each of the samples but for both dates, the dominant phytoplankton type was approximately an order of magnitude greater in number than the other phytoplankton types. Transect samples, acquired on the 2 August 1973 and 9 August 1973 were of chlorophyll-a concentration which were converted to phytoplankton concentration using a typical ratio (Reynolds 1984).
No data were available on the initial horizontal or vertical distribution of either of these phytoplankton types. The initial vertical distribution of dinoflagellates was parameterised with a maximum at 4 m, with decreases to near zero at the surface and a depth of 8 m. This distribution was chosen to satisfy requirements for light and oxygen. For the levels of sub-surface irradiance that were estimated to exist at the start of the model's simulation (around 600 micromolar photons m-2 s-1), levels of approximately 150 micromolar photons m-2 s-1) would correspond to a depth of approximately 3 to 4 m for Esthwaite Water, assuming that there were no phytoplankton at lesser depths in the column.
The initial vertical distribution of cyanobacteria was parameterised with a maximum at the surface and a minimum at 8 m depth. This distribution was simulated according to the estimated change in irradiance levels, buoyancy and rising velocity that had occurred over the previous week (with the assumption that the initial mean cyanobacteria population in the upper 8 m of the water column was 7 g m-3 (chlorophyll-a concentration of 80 mg m-3).
1. Dinoflagellates dominant with greatest concentration at thermocline (2 August 1973).
Surface currents were deflected clockwise from the wind direction (Figure 2a) by the Coriolis force. This caused downwelling at the downwind area of the lake, particularly in the east and south (Figure 2b) and upwelling at the upwind area of the lake, particularly in the west and north. Vertical current speeds were of the order of 0.0002 m s-1). There was a progressive rotation in current direction with depth due to a decreased current speed as the thermocline was approached (Figure 3). Immediately above the thermocline, currents moved towards the south at slightly under 0.005 m s-1). Beneath the thermocline, current continued to be deflected clockwise.
(a) 
(b) 
Figure 2. Simulated currents of Esthwaite (15:00 2 August 1973): (a) surface current vectors; (b) vertical velocity at surface.
Figure 3. Current vectors for Esthwaite Water (15:00 2 August 1973) as a function of depth: values at the end of each arrow are the depth (m), values on x and y axes are u and v velocity component (m s-1).
Most velocity vectors near to the surface were similar to those measured. The rate of rotation with depth was underestimated, however, and the effect of the density gradient at the thermocline on flows was negligible.
Measured and simulated surface dinoflagellates distributions were similar (Figure 4), although the simulation estimated an east-west trend that was not observed in the field.
(a) 
(b) 
Figure 4. Surface dinoflagellate concentration (15:00 2 August 1973): (a) measured; (b) simulated.
2. Cyanobacteria dominant with greatest concentration at surface (9 August 1973).
For this example, surface currents were deflected to clockwise from the wind direction (Figure 5a) by the Coriolis force, but the deflection was to a lesser extent than on 2 August 1973 because of the greater wind speed. Downwelling again occurred at the downcurrent area of the lake in the east (Figure 5b) and upwelling again occurred at the upwind area of the lake in the west, but the magnitude was greater than previously (0.0004 m s-1) because of the increased horizontal displacement by wind.
(a) 
(b) 
Figure 5. Simulated current of Esthwaite (15:00 9 August 1973): (a) surface current vectors; (b) vertical velocity at surface.
The cyanobacteria increased in density but remained positively buoyant throughout most of the simulation. The rate of density increase was inversely proportional to distance from the surface because light intensity, and therefore photosynthesis, decreased with depth.
The positively buoyant cyanobacteria accumulated over downwelling areas, resulting in a similar surface distribution to that measured (Figure 6) although with a greater magnitude of variation.
(a) 
(b) 
Figure 6. Surface cyanobacterial concentration (15:00 9 August 1973): (a) measured; (b) simulated.
Any discrepancy between the simulated and measured distributions may be attributed to simulation errors where the hydrodynamics had not been modelled with enough accuracy, or to sampling errors where the sampling scheme had been unable to resolve existing variation or the locations of the sample observations were not known with satisfactory accuracy. The horizontal sample transects undertaken by George and Heaney (1978) covered most of the areas where the discrepancies occurred (particularly in the dinoflagellate simulation) so that horizontal discrepancies may be attributed to simulation uncertainty. Furthermore, although there will have been error in determining the exact positions of the observations, this will have been small relative to the scale over which the gradient was observed.
1. Dinoflagellates dominant with greatest concentration at thermocline (2 August 1973).
The discrepancies between the simulated and measured horizontal distributions may be attributed to the model over-estimating the rate of downwards momentum transfer, and therefore under-estimating the rate of angular deflection with increasing depth. The simulated vertical velocity profile (Figure 3) had flows occurring at a depth of approximately 3 m that have hardly been deflected from flows at the surface. With a peak in dinoflagellate concentration occurring at a depth of 3 m, some dinoflagellates would be entrained in downwind drift currents.
2. Cyanobacteria dominant with greatest concentration at surface (9 August 1973)
Some of the discrepancies between the simulated and measured horizontal distributions may be related to the sampling regime which in the north missed the area of simulated high concentration, and near to the southern outlet missed the areas of simulated low and high concentrations on the west and east shores, respectively. Sampling the simulated distribution with the regime used on 9 August 1973 and interpolating a new distribution would create one that was more similar to that measured. Additionally, use of fine spatial resolution remote sensing data (George 1993) has shown patterns consistent with those simulated , so the simulated variation is not atypical for this lake. However, the position of the simulated peak concentration several hundred metres north of that measured suggests that the eastward component of the near surface velocity vector was too small relative to the northward component in the north-east part of the lake. The large underestimation of surface concentration in the north-west bay may likewise be attributed to inaccurate estimation of the velocity field - a problem that has been encountered in past CFD modelling of Esthwaite Water (Falconer et al., 1991).
This study has shown that even without comprehensive data on the environment of Esthwaite Water, through the use of simple deterministic models and approximate estimates of boundary conditions it was possible to simulate the synoptic features of the spatial distribution of phytoplankton. Such features included: an upwind accumulation of surface-avoiding dinoflagellates, and a downwind accumulation of positively buoyant cyanobacteria. Discrepancies between measured and simulated distributions showed the sensitivity of the lake system to environmental conditions. Simulations (not presented here) using slightly different external boundary conditions (for instance, different estimates of irradiance, different initial vertical phytoplankton profiles) or different internal boundary conditions (for instance, different turbulence models) produced different, but generally similar, spatial distributions. With more comprehensive data on the initial boundary conditions, it is expected that the system can be described with greater accuracy.
This project was supported by funding from the Natural Environment Research Council and the Institute of Freshwater Ecology.
Belov, A.P., Giles, J.D., 1997, Dynamical model of buoyant cyanobacteria, Hydrobiologia, 349, 87-97.
Bengtsson, L., 1973, Wind stress on small lakes, Teniska Hxgskolan, Lund, Sweden.
Falconer, R.A., George, D.G., Hall, P., 1991, Three-dimensional numerical modelling of wind driven circulation in a shallow homogenous lake, Journal of Hydrology, 124, 59-79.
George, D.G., 1993, Physical and chemical scales of pattern in freshwater lakes and reservoirs, The Science of the Total Environment, 135, 1-15.
George, D.G., Edwards, R.W., 1976, The effect of wind on the distribution of chlorophyll-a and crustacean plankton in a shallow eutrophic reservoir, Journal of Applied Ecology, 13, 667-690.
George, D.G., Heaney, S.I., 1978, Factors influencing the spatial distribution of phytoplankton in a small productive lake, Journal of Ecology, 66, 135-155.
Heaney, S.I., 1976, Temporal and spatial distribution of the dinoflagellate Ceratium hirundinella O.F. Muller within a productive lake, Freshwater Biology, 6, 531-542.
Heaney, S.I., Furnass, T.I., 1980, Laboratory models of diel vertical migration in the dinoflagellate Ceratium hirundinella, Freshwater Biology, 10, 163-170.
Heaney, S.I., Talling, J.F., 1980, Dynamic aspects of dinoflagellate distribution patterns in a small productive lake, Journal of Ecology, 68, 75-94.
Howard, A., Kirby, M.J., Kneale, P.E., McDonald, A.T., 1995, Modelling the growth of cyanobacteria (GrowScum), Hydrological Processes, 9, 809-821.
Kirk, T.O., 1994, Light and Photosynthsis in Aquatic Ecosystems, Cambridge University Press, Cambridge, pp509.
Kromkamp, J, Walsby, A.E., 1990, A computer model of buoyancy and vertical migration in cyanobacteria. Journal of Plankton Research, 12, 1,161-183.
Melaaen, M. C., 1992, Calculation of fluid flows with staggered and nonstaggered curvilinear nonorthogonal grids - the theory, Numerical Heat Transfer, Part B, vol. 21, pp 1- 19.
Olsen, N.R.B., 1997, Compuational Fluid Dynamics in Hydraulic and Sedimentation Engineering, Division of Hydraulic and Environmental Engineering, The Norwegian University of Science and Technology, pp 57.
Patankar, 1980, Numerical Heat Transfer and Fluid Flow , Hemisphere Publishing Corporation, Taylor and Francis, New York.
Reynolds, C.S., 1984, The Ecology of Freshwater Phytoplankton, Cambridge University Press, Cambridge, pp384.