e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 available at www.sciencedirect.com journal homepage: www.elsevier.com/locate/ecolmodel Modelling spatial heterogeneity of phytoplankton in Lake Mangueira, a large shallow subtropical lake in South Brazil Carlos R. Fragoso Jr. a,∗ , David M.L. Motta Marques a , Walter Collischonn a , Carlos E.M. Tucci a , Egbert H. van Nes b a Federal University of Rio Grande do Sul, Hydraulic Research Institute, CP 15029 Porto Alegre, RS, Brazil Wageningen University, Aquatic Ecology and Water Quality Management Group, Department of Environmental Sciences, P.O. Box 47, 6700 AA Wageningen, The Netherlands b a r t i c l e i n f o Article history: a b s t r a c t We present a model describing phytoplankton growth in Lake Mangueira, a large subtropical Received 22 October 2007 lake in the Taim Hydrological System in South Brazil (817 km2 , average depth 2 m). The hori- Received in revised form zontal 2D model consists of three modules: (a) a detailed hydrodynamic module for shallow 7 August 2008 water, which deals with wind-driven quantitative flows and water level, (b) a nutrient mod- Accepted 21 August 2008 ule, which deals with nutrient transport mechanisms and some conversion processes and Published on line 25 September 2008 (c) a biological module, which describes phytoplankton growth in a simple way. We solved the partial differential equations numerically by applying an efficient semi-implicit finite Keywords: differences method to a regular grid. Hydrodynamic parameters were calibrated to continu- Subtropical shallow lake ous measurements of the water level at two different locations of the lake. An independent Wetland validation data set showed a good fit of the hydrodynamic module (R2 ≥ 0.92). The nutrient Phytoplankton patchiness and biological modules were parameterized using literature data and verified by compar- Spatial heterogeneity ing simulated phytoplankton patterns with remote sensing data from satellite images and Complex model field data of chlorophyll a. Moreover, a sensitivity analyses showed which parameters had the largest influence on the simulated phytoplankton biomass. The model could identify zones with a higher potential for eutrophication. It has shown to be a first step towards a management tool for prediction of the trophic state in subtropical lakes, estuaries and reservoirs. © 2008 Elsevier B.V. All rights reserved. 1. Introduction During the last 200 years, many lakes have suffered from eutrophication, implying an increase of both the nutrient loading and organic matter (Wetzel, 1996). This degradation process usually resulted in an increase of water turbidity due to blooms of cyanobacteria or green algae and subsequently the disappearance of submerged aquatic macrophytes (Moss, 1998). ∗ An aspect that has often been neglected in freshwater systems is the fact that phytoplankton is horizontally often not evenly distributed in space. Though the occurrence of phytoplankton patchiness is known for a long time in marine systems (e.g. Platt et al., 1970; Steele, 1978; Steele and Henderson, 1992), phytoplankton in shallow lakes is often assumed to be homogeneous. However, there are various mechanisms that may cause horizontal heterogeneity in shallow lakes. For example, grazing by aggregated zooplank- Corresponding author. Tel.: +55 51 33086654; fax: +55 51 33087509. E-mail addresses: crubertofj@hotmail.com (C.R. Fragoso Jr.), dmm@iph.ufrgs.br (D.M.L.M. Marques), collischonn@iph.ufrgs.br (W. Collischonn), tucci@iph.ufrgs.br (C.E.M. Tucci), egbert.vannes@wur.nl (E.H. van Nes). 0304-3800/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2008.08.004 126 e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 Fig. 1 – Taim Hydrological System situated in the South of Brazil. The meteorological stations in the north, middle and south of Lake Mangueira are denominated as TAMAN, TAMAC and TAMAS, respectively. ton and other organisms may cause spatial heterogeneity in phytoplankton (Scheffer and De Boer, 1995). Submerged macrophytes beds may be another mechanism by reduction of wave resuspension and allopathic effect on the algal community (Van den Berg et al., 1998). For large shallow lakes, wind can be a dominant factor leading to both spatial and temporal heterogeneity of phytoplankton (Carrick et al., 1993), either indirectly by affecting the local nutrient concentration due to resuspended particles, or directly by resuspending algae from the sediment (Scheffer, 1998). In the management of large lakes, prediction of distributed phytoplankton can assist the manager to decide on an optimal course of actions, such as biomanipulation and regulation of recreation or potable water supply (Reynolds, 1999). However, it is difficult to measure the spatial distribution of phytoplankton. Mathematical modelling of phytoplankton population can be an important alternative methodology in improving our knowledge regarding the physical, chemistry and biological processes related to the phytoplankton ecology (Scheffer, 1998; Edwards and Brindley, 1999; Mukhopadhyay and Bhattacharyya, 2006). There is already a large variety of phytoplankton models. The simplest models are based on a steady-state or on the assumption of complete mixing (Schindler, 1975; Smith, 1980; Thoman and Segna, 1980). Phytoplankton models based on more complex vertical 1D hydrodynamic processes gave a more realistic representation of the stratification and mixing processes in deep lakes (Imberger and Patterson, 1990; Hamilton et al., 1995a,b; Imberger, 1995). However, the vertical 1D assumption might be restrictive, especially in large shallow lakes that are poorly stratified and often characterized by a significant differences between pelagic and shore zones. In those cases, a horizontal 2D model with a complete description of the hydrodynamic and ecological processes can offer more insight in the factors determining local water quality. Currently, computational power is not limiting the development of 2D and 3D models anymore, and these models are more frequently applied. There is a large diversity of 2D and 3D hydrodynamic models, of which most are designed to study deep-ocean circulation or coastal, estuarine and lagunal zones (Blumberg and Mellor, 1987; Casulli, 1990). However, only few of them are coupled with biological components (Lord et al., 1994; Rajar and Cetina, 1997; Bonnet and Wessen, 2001). Over the past decade there has been a concerted effort to increase the realism of ecosystem models that describe plankton production as biological indicator of eutrophication. Most of this effort has been expended on the description of phytoplankton in temperate lakes; thus, multi-nutrient, photo acclimation models are now not uncommon (e.g. Olsen and Willen, 1980; Edmondson and Lehman, 1981; Sas, 1989; Fasham et al., 2006; Mitra and Flynn, 2007; Mitra et al., 2007). In subtropical lakes, eutrophication has been intensively studied, but only focusing on measurement of changes in nutrient concentrations (e.g. Matveev and Matveeva, 2005; Kamenir et al., 2007). Here, we present a phytoplankton model coupled with a horizontal 2D hydrodynamic model for the large subtropical shallow Lake Mangueira (South Brazil) focusing on spatial heterogeneity of phytoplankton. The hydrodynamic parameters were calibrated using continuous water level measurements in two stations in Lake Mangueira and the Taim wetland. The parameters describing algae growth were based on literature values. The generated spatial patterns of chlorophyll a concentration were verified both with a field data set and with a cloud-free satellite image provided by Terra Moderate Resolution Imaging Spectroradiometer (MODIS) with spatial 127 e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 The hydrodynamic model is based on the shallow water equations derived from Navier–Stokes, which describe dynamically a horizontal two-dimensional flow: Fig. 2 – Simplified representation of the interactions involving the state variables (double circle) and the processes (rectangle). resolution of 1.0 km. Additionally, we carried out a sensitivity analysis of the biological parameters. 2. Material and methods 2.1. Study area The Taim Ecological Reserve was established to protect the Taim Hydrological System, a heterogeneous and productive ecosystem in southern Brazil (Fig. 1), harboring an exceptional high biodiversity. The reserve encompasses a variety of habitats such as beaches, dunes, forests, grasslands, lakes and wetlands (Garcia et al., 2006). The system studied is part of this reserve and includes the Taim wetland and Lake Mangueira covering a total area of about 950 km2 , of which 86% is occupied by the lake (Fig. 1). The average depth of the lake is approximately 2 m and its trophic state ranges from oligotrophic to mesotrophic. The mesotrophic conditions occur in the spring and summer when it suffers from a notable water withdrawal to irrigation of rice crops (approximately 2 L ha−1 s−1 during 100 days), as well as a high input of nutrients loading from its watershed. The system is located in a subtropical climate region. Hydrological and meteorological variables have being monitored at hourly frequency since 2001 by a federal conservation program. These variables include water level, precipitation, temperature, solar radiation, wind velocity and direction in three stations in Lake Mangueira (Fig. 1). 2.2. Model description The present model was structured in three modules: (a) hydrodynamic module, which simulates the variables that quantify the water flow (velocity and water level), (b) transport and conversions of nutrients, and (c) a biological module, which deals with growth and loss of phytoplankton biomass. An overview of the modeled processes is given in Fig. 2. ∂[(h + )u] ∂[(h + )v] ∂ + + =0 ∂t ∂x ∂y (1) ∂u ∂u ∂u ∂ +u +v = −g − u + x + Ah ∇ 2 u + fv ∂t ∂x ∂y ∂x (2) ∂v ∂v ∂ ∂v +u +v = −g − v + y + Ah ∇ 2 v − fu ∂t ∂x ∂y ∂y (3) where u(x, y, t) and v(x, y, t) are the water velocity components in the horizontal x and y directions; t is time; (x, y, t) is the water surface elevation relative to the undisturbed water surface; g is the gravitational acceleration; h(x, y) is the water depth measured from the undisturbed water surface; f is the parameter of Coriolis;  x and  y are the wind stresses in the x and y directions; ∇ = ∂/∂x · i + ∂/∂x · j is a vector operator in the plane x–y;  Ah is the coefficient of horizontal eddy viscosity; and = (g u2 + v2 )/Cz (Daily and Harlerman, 1966) where Cz is the Chezy friction coefficient. Usually, the wind stresses in the x and y directions are written as a function of wind velocity (Wu, 1982): x = CD · Wx · ||W|| (4) y = CD · Wy · ||W|| (5) where CD is the wind friction coefficient; Wx and Wy are the wind velocity components (m s−1 ) in the x and y directions, respectively. Wind velocity is measured at 10 m from water surface; ||W|| = Wx2 + Wy2 is the norm of wind velocity vector. An efficient numerical semi-implicit Eulerian–Lagrangian finite differences scheme was used in order to assure stability, convergence and accuracy (Casulli, 1990; Casulli and Cheng, 1990; Casulli and Cattani, 1994). The nutrient module considers the advection and diffusion of each substance, inlet and outlet loading, sedimentation and resuspension through following equation: ∂(HC) ∂(uCH) ∂(vCH) ∂ + + = ∂t ∂x ∂y ∂x  Kh ∂(HC) ∂x  + + source or sink ∂ ∂y  Kh ∂(HC) ∂y  (6) where C is the average concentration in the water column; H =  + h is the total depth; and Kh is the horizontal scalar diffusivity assumed as 0.1 m2 day−1 (Chapra, 1997). Eq. (6) was applied to model total phosphorus, total nitrogen and phytoplankton. All these equations are solved dynamically, using a simple numerical semi-implicit central finite differences scheme (Gross et al., 1999a,b) (Fig. 2). Thus, the mass balances involving phytoplankton and nutrients can 128 e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 be written as ∂(Ha) ∂(uHa) ∂(vHa) ∂ + + = eff Ha + ∂t ∂x ∂y ∂x + ∂ ∂y  Kh  Kh ∂(Ha) ∂y ∂(Ha) ∂x   + inlet/outlet (7a) ∂(Hn) ∂(uHn) ∂(vHn) ∂ + + = −ana eff Ha + ∂t ∂x ∂y ∂x + ∂ ∂y  Kh ∂(Hn) ∂y   ∂(Hn) Kh ∂x T = Gmax · TT−20  + inlet/outlet (7b) ∂(Hp) ∂(uHp) ∂(vHp) ∂ + + = −apa eff Ha−kphos p+ ∂t ∂x ∂y ∂x ∂ + ∂y  ∂(Hp) Kh ∂y   Kh where P (T, N, I) is the primary production rate as a function of temperature (T), nutrients (N), and light (I); L is the loss rate due respiration, excretion and grazing by zooplankton and a is the chlorophyll a concentration. The temperature effect on primary production was assumed to be an exponential function that is widely used in phytoplankton models (e.g. Eppley, 1972; Canale and Vogel, 1974), which presents only two parameters to calibrate. ∂(Hp) ∂x  where T is the growth rate (day−1 ) at temperature T (◦ C); Gmax is the maximum growth rate algae at 20 ◦ C; and  T is the temperature effect coefficient. We used a commonly used Monod saturating function to model nutrient limitation. In our case, involving multiple nutrients, there are several ways in which the nutrient limitation term could be refined. We used Liebig’s law, where the most limiting nutrient controls phytoplankton growth rate (Lucas, 1997): + inlet/outlet (7c) where a, n and p are chlorophyll a, total nitrogen and total phosphorus concentrations, respectively; ana is the N/Chla ratio equal to 8 mg N mg Chla−1 ; apa is the P/Chla ratio equal to 1.5 mg N mg Chla−1 , inlet/outlet represents the balance between all inlets and outlets in a control volume ∂x ∂y ∂z; and kphos is the settling coefficient of the phosphorus which can be estimated by (Chapra, 1997): N = min ∂(Hn) ∂(uHn) ∂(vHn) ∂ + = −ana eff Ha + + ∂t ∂x ∂y ∂x + ∂ ∂y  Kh ∂(Hn) ∂y  Kh ∂(Hn) ∂x  Four important assumptions were made in this scheme: (a) fixed stoichiometric conversions were applied to transfers between nutrients and phytoplankton; (b) there is no loss of mass due to degradation processes, (c) total nitrogen in the water was considered as a conservative substance, we thus discarded processes like denitrification, N fixation, sedimentation and resuspension and (d) total phosphorus was assumed to fixate to the sediment, resuspension was neglected. The phytoplankton growth (primary production) and loss processes are represented through effective growth rate (Lucas, 1997). The effective growth rate itself is not a simple constant, but varies in response to environmental factors such as temperature, nutrients, respiration, excretion and grazing by zooplankton: ∂(Hp) ∂(uHp) ∂(vHp) + + = −apa eff Ha − kphos p ∂t ∂x ∂y ∂(Hp) Kh ∂x + inlet/outlet  ∂ + ∂y   2.718fp (e−˛1 − e−˛2 ) ke ( + h) ˛1 = fPAR Ia + inlet/outlet ˛2 =  n p , n + kN p + kP (11) (12) where, (8) ∂ + ∂x  where N is the growth rate due nutrients uptake (day−1 ); and kN and kP are the half-saturation for nitrogen and phosphorus uptake, respectively. The dependence of the growth rate of phytoplankton on light was approached by an optimum function (Steele, 1965), incorporating light inhibition a high light levels. L =  (10) ∂(Hp) Kh ∂y  (9) e−ke (+h) Is fPAR Ia Is (13) (14) where L is the phytoplankton growth rate on light dependence (day−1 ); fp is the photoperiod; fPAR is the fraction of photosynthetically active radiation (PAR); Ia is the light level (kJ m−2 day−1 ); IS is the optimal light level (kJ m−2 day−1 ), and ke is the light attenuation coefficient (m−1 ). We assumed that light attenuation coefficient can be related with the amount of chlorophyll a concentration in the water (Riley, 1956): ke = k′e + 0.0088a + 0.054a2/3 (15) where k′e is the light attenuation (m−1 ) due other factors rather than phytoplankton, such as particle-free water and color, nonvolatile suspended solids and detritus. The processes contribute to the loss rate of phytoplankton considered were respiration, excretion and grazing by zooplankton and others organisms. They are usually modeled as a single first-order decay being respiration and excretion depending of temperature, as in L = R + G (16) 129 e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 Table 1 – Hydrodynamic and biological parameters description and its values range Parameter Description Unit Values range Hydrodynamic 1 Ah 2 CD 3 CZ Horizontal eddy viscosity coefficient Wind friction coefficient Chezy coefficient m1/2 s−1 – – Biological 1 Gmax 2 IS 3 k′e 4 T 5 R 6 kP 7 kN 8 kre 9 kgz Maximum growth rate algae Optimum light intensity for the algae growth Light attenuation coefficient in the water Temperature effect coefficient Respiration and excretion effect coefficient Half-saturation for uptake phosphorus Half-saturation for uptake nitrogen Respiration and excretion rate Zooplankton grazing rate day−1 cal cm−2 dia−1 m−1 – – mg P m−3 mg N m−3 day−1 day−1 and R = kre · RT−20 (17) G = kgz (18) where L is the total phytoplankton loss rate; R is the loss rate of phytoplankton by respiration and excretion; G is the death of algae due to grazing by zooplankton and other organisms; kre is the respiration and excretion rate; kgz is the grazing rate by zooplankton and other organisms and  R is a coefficient modelling the temperature effect. Thereby, nine parameters control the variation of the effective growth rate of phytoplankton. These biological parameters and their respective values range are listed in Table 1. 2.3. Calibration and validation of hydrodynamic module The basic part of the model is the hydrodynamic module. Thus, an accurate prediction of the hydrodynamic conditions can identify how the phytoplankton is being transported and where zones with high potential of eutrophication and phytoplankton blooms are located. The hydrodynamic module was calibrated by tuning the model parameters within their observed literature ranges (Table 1). Nonetheless, the hydraulic resistance caused by presence of emerged macrophytes in Taim Wetland was represented by a smaller Chezy’s resistance factor than used in other lake areas (Wu et al., 1999). Calibration and validation of the hydrodynamical parameters was done using two different time-series of water level and wind produced for two locations in Lake Mangueira (north and south). We used a period of 26 days for calibration, starting 10 July 2002 at 4:00 p.m., and 15 days for validation, starting 01 January 2003 at 0:00 a.m. The reason of that choice was the availability of continuous data during these periods. The valibration period differed significantly from the calibration period as in this summer period withdrawal of water to rice crops took place. We assumed a constant pumping rate to represent water withdrawal to rice crops, whereas in reality there were unknown daily fluctuations. 5–15 2e−6–4e−6 50–70 1.5–3.0 100–400 0.25–0.65 1.02–1.14 1.02–1.14 1–5 5–20 0.05–0.25 0.10–0.20 Reference White (1974) Wu (1982) Chow (1959) Jørgensen (1994) Schladow and Hamilton (1997) Schladow and Hamilton (1997) Eppley (1972) Schladow and Hamilton (1997) Lucas (1997) Lucas (1997) Chapra (1997) Chapra (1997) The coefficient of determination (R2 ) between model results and field data was computed to measure the model performance. 2.4. Phytoplankton simulation For the parameters of the phytoplankton module we used the average values for the literature range given in Table 1. To evaluate its performance we simulated another period of 86 days, starting 22 December 2002 at 00:00 h (summer). Solar radiation and water temperature data were taken from TAMAN meteorological station, situated in northern part of Lake Mangueira. Photosynthetically active radiation (PAR) at the Taim wetland was assigned as 20% of the total radiation, in order to represent the indirect effect of the emergent macrophytes on phytoplankton growth rate according experimental studies of emerged vegetation stands in situ. At the lake areas, we assumed that the percentage of PAR was 50% of the total solar radiation (Janse, 2005). The resulting phytoplankton patterns were compared with satellite images from MODIS, which provides improved chlorophyll a measurement capabilities over previous satellite sensors. For instance, MODIS can better measure the concentration of chlorophyll a associated with a given phytoplankton bloom. Unfortunately, there were no detailed chlorophyll a and nutrient data available for the same period. Therefore, we compared only the median simulated values with field data from another period (2001 and 2002). 2.5. Sensitivity analysis of the biological module To determine which biological parameters had the strongest effect on the calculated phytoplankton biomass, a simple sensitivity analysis was applied. We considered only biological parameters which are directly responsible for effective growth rate (d−1 ) of phytoplankton. For this analysis we use the same conditions and input variables that were used for hydrodynamic calibration (26 days). For each parameter the model ran twice: first with the minimum value and then with the maximum value. All other parameters were kept at their default values. For each of these runs, the minimum, maximum and mean daily chlorophyll a concentrations were calculated for the simulation period. 130 e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 Fig. 3 – Time series of wind velocity and direction on Lake Mangueira and water levels fitted at the north and south of Lake Mangueira into calibration and validation periods (simple line: observed, dotted line: calculated). 3. Results 3.1. Calibration and validation of hydrodynamic module The simulated and observed values of water levels at two stations of Lake Mangueira during the calibration and validation period are shown in Fig. 3. The model was able to reproduce the water level well in both extremities of Lake Mangueira. Wind-induced currents can be considered the dominant factor controlling transport of substances and phytoplankton in Lake Mangueira, producing advective movement of superficial water masses in a downwind direction. For instance, a southwest wind, with magnitude approximately greater than 4 m s−1 , can causes a significant transport of water mass and substances from south to north of Lake Mangueira, leading to a almost instantaneous increase of the water level in the northeastern parts and, hence the decrease of water level in southwestern areas. Apart from wind effects, also the water balance between precipitation and evaporation is an important factor determining the water level in Lake Mangueira. Note that in the validation period the water withdrawal to rice crops was an important factor in the water balance. The assumption of a constant pumping rate turned out to be a reasonable approach, which it might have led to the slightly lower coefficients of determination (R2 ). In summary, the fit and verification of the model have showed that a semi-implicit Eulerian–Lagrangian difference finite scheme presents good results, respecting stability, convergence and precision principles as it was previously seen e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 Fig. 4 – Relative effect of each biological parameter on the minimum, mean and maximum daily values of chlorophyll a considering its (a) maximum and (b) minimum parameters values from ranges presented in Table 1. in others papers (e.g. Casulli and Cheng, 1990; Casulli and Cattani, 1994). 3.2. Sensitivity analysis The results of the sensitivity analysis are summarized in Fig. 4. The first panel corresponds to relative effect of the maximum value assigned in the range for a particular parameter (Fig. 4a) and the second one using its minimum value (Fig. 4b). The histogram shows that the model outcome is rather sensitive to various parameters. The parameters with a strongest effect in the model outcome are related to the effect of temperature (Gmax ,  T and  R ) and to loss processes, such as respiration and zooplankton grazing (kra and kgz ). The parameters related to light penetration (Is and k′e ) and phosphorous/nitrogen uptake (kP , kN ) have a weaker effect on the effective growth rate, indicating that the phytoplankton growth in subtropical conditions (summer-autumn period) was not strongly limited by light and nutrients. 3.3. Phytoplankton simulation in the Lake Mangueira The model was used to determinate the spatial distribution of chlorophyll a and to identify locations with higher growth 131 and phytoplankton biomass in Lake Mangueira. Fig. 5 shows the spatial distribution of phytoplankton biomass for different times during the simulation period. Specifically, in Lake Mangueira there is a strong gradient of phytoplankton productivity from littoral to pelagic zone (Fig. 5). Moreover, the model outcome suggests that there is a significant transport of phytoplankton and nutrients from littoral to pelagic zones through hydrodynamic processes. This transport was intensified by several large sand bank formations that are formed perpendicular to the shore line of the lake, carrying nutrients and phytoplankton from the shallow to deepest zones. In addition, it was possible to identify zones with the highest productivity. There is a trend of phytoplankton aggregation in the southwest and northeast areas as dominant wind directions coincide with longitudinal direction of the Lake Mangueira (Fig. 6). The clear water in the Taim wetland, north of Lake Mangueira, was caused by shading of emergent macrophytes, modeled as a fixed reduction of PAR. After 1200 h of simulation (50 days), the daily balance between the total primary production and loss was negative. That means that daily losses such as respiration, excretion and grazing by zooplankton exceeded the primary production in the photoperiod, leading to a significant reduction of the chlorophyll a concentration for the whole system (Fig. 5d and e). We verified the modeled spatial distribution of chlorophyll a with those estimated by remote sensing (Fig. 7). The simulated patterns had a reasonably good similarity with the evaluated patterns by remote sensing (Fig. 7a and b). In both figures, large phytoplankton aggregations can be observed in both southern and northern parts of Lake Mangueira, as well as in littoral zones. Unfortunately we did not have independent data of phytoplankton in the simulation period. Therefore we could only compare the median values of simulated and observed chlorophyll a, total nitrogen and total phosphorus for three points of Lake Mangueira (Fig. 8). We thus assumed that the median values were comparable between the years. The fit of these variables was reasonable, considering that we did not calibrate the biological parameters of the phytoplankton module. We also fitted the model without considering spatial processes on the median chlorophyll a data of all stations (results not shown), but obviously it is then impossible to model the differences between the stations. Interestingly, this nonspatial simulation resulted in a poorer fit, as the simulated median chlorophyll a values were systematically overestimated (ca. 28 mg m−3 ). That indicates that in the model hydrodynamic processes (i.e. advection and diffusion) had a rather strong effect on the average pelagical chlorophyll a values, probably due to the exchange with the littoral zones. 4. Discussion Recognition of the importance of spatial and temporal scales is a relatively recent issue in ecological research on aquatic food webs (Bertolo et al., 1999; Woodward and Hildrew, 2002; Bell et al., 2003; Mehner et al., 2005). Among other things, the 132 e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 Fig. 5 – Phytoplankton dry weight concentration fields in ␮g l−1 , for the whole system at different times: (a) 14 days; (b) 28 days; (c) 43 days; (d) 57 days; (e) 71 days; (f) 86 days. The color bar indicates the phytoplankton biomass values. A wind sleeve, in each frame, indicates the direction and the intensity of the wind. The border between Taim wetland and Lake Mangueira is showed as well. observational or analytical resolution necessary for identifying spatial and temporal heterogeneity in the distributions of populations is an important issue (Dungan et al., 2002). Most ecological systems exhibit heterogeneity and patchiness on a broad range of scales, and this patchiness is fundamental to population dynamics, community organization and stability. Therefore, ecological investigations require an explicit determination of spatial scales (Levin, 1992; Hölker and Breckling, 2002), and it is essential to incorporate spatial heterogeneity in ecological models to improve understanding of ecological processes and patterns (Hastings, 1990; Jørgensen et al., 2008). In models of lake ecosystems, horizontal spatial heterogene- e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 133 Fig. 6 – The wind rose, showing the distribution of wind speeds, and the frequency of the varying wind directions. ity of phytoplankton and hydrodynamic processes are often neglected. Our model analysis showed that it is important to consider such spatial heterogeneity in large lakes, as the water quality is expected to differ significantly between the shores and the pelagic zones. Especially for prediction of the water quality (including the variability due to wind) at the littoral zones of such lake, incorporation of spatial explicit processes is essential. Such information is often important for recreationists and lake managers. Also for more detailed studies of the growth and competition of phytoplankton species, detailed information of water movement is very important (Huisman et al., 1999). Fig. 7 – (a) MODIS-derived chlorophyll a image of 1 km spatial resolution on 8 February 2003 and (b) chlorophyll a concentration field simulated in Lake Mangueira on 8 February 2003. 134 e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 Fig. 8 – Comparison between box-plot diagrams corresponding to 37 samples (2001 and 2002) and the median values of chlorophyll a, total nitrogen and total phosphorus simulated by model at the three sampling locations of Lake Mangueira (TAMAS = south, TAMAC = centre and TAMAN = north). The hydrodynamic module, using a semi-implicit Eulerian–Lagrangian finite differences scheme, showed excellent results during hydrodynamic calibration and validation (i.e. coefficients of determination above 0.92). The scheme, which is rarely used in ecological researches (e.g. Lucas and Cloern, 2002; Robson and Hamilton, 2004; Romero et al., 2004a; Romero et al., 2004b; Spillman et al., 2007), allowed us to use a larger step time than in other numerical schemes, assuring stability, convergence and precision (Casulli and Cheng, 1990; Casulli and Cattani, 1994). Computationally, the resulting algorithm is suitable for the simulations of complex two or three-dimensional flow using fine spatial resolution and relatively large time steps. The present formulation also is fully vectorizable and allows for the simulation of flooding and drying of tidal flats. The model describes the most important hydrodynamic and the main biological processes of the phytoplankton in an integrated way, in order to help to understand the role of phytoplankton heterogeneity in a shallow lake. However, as the model focuses on hydrodynamic processes many biological processes were simplified. For instance, important microbial processes such as nitrification, resuspension, mineralization of detritus and interactions between water-sediment were neglected or strongly simplified. The modelling of phytoplankton growth was also kept very simple. All species of phytoplankton were lumped without distinction between functional groups, such as cyanobacteria, diatoms and other small edible algae. Therefore, their different characteristics could not be taken account. Other aquatic organisms, such as phytoplankton in the sediment, zooplankton, zoobenthos, macrophytes and fishes, were not modeled dynamically, limiting the interactions between trophic groups. The interactions with emergent macrophytes in wetland area were also only indirectly modeled by a fixed reduction of PAR. Others more realistic mechanisms can also be implemented for wetland areas, such as moderate effect of wind-induced resuspension and distinct nutrient processes in the water column. The current model could be extended to provide a more complete description of aquatic food-web, but the complexity involved in a full food-web model may be large, whereas a simpler model allows for better understanding of the main processes that determine the results (Van Nes and Scheffer, 2005). As it is good to combine different approaches of different complexity (Van Nes and Scheffer, 2005), we are currently developing a more complex version of this model which aggregates most neglected biological processes (Fragoso et al., 2007). Our sensitivity analysis, analyzing the summer-autumn period, indicated that phytoplankton growth in this subtropical shallow lake is particularly sensitive to temperature and loss processes like respiration and grazing by zooplankton and other organisms. This is in line with model analyses and empirical results in temperate lakes (Barko and Smart, 1981; Spencer, 1986; Scheffer et al., 1993). Of course, subtropical lakes do not freeze in winter, but they can also be subject to relatively strong temperature variations during the season (e.g. in Lake Mangueira, water temperature can take values between 8 ◦ C and 27 ◦ C). In view of the key role of phytoplankton in subtropical lakes, the relatively strong effect of temperature on phytoplankton biomass indicates that climate e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 changes can influence the tropic state of the subtropical lakes as was observed for temperate lakes (Scheffer et al., 2001; Van Leeuwen et al., 2007). Remote sensing data can help us to calibrate and verify a distributed heterogeneously model outcome. However, the application of satellite remote sensing for lake water is constrained by the need for high spatial resolution image data and thus remains limited by spectral resolution capabilities. Furthermore, it is difficult to quantify chlorophyll a in waters characterized by high and heterogeneous suspended sediment concentrations (SSC). The SSC dominates the spectral reflectance, masking the spectral influence from other components in broad spectral band systems, making chlorophyll a determination from remote sensing imagery difficult. Also in Lake Mangueira, SSC can be a dominate component in the water column, as resuspension can be significant in this shallow windy lake. Apart of results from remote sensing, we observed that an independent data set of the spatial distribution of chlorophyll a and other states variables is important for a better verification of the modeled phytoplankton results. The lack of spatial and temporal distributed data for the Lake Mangueira made it impossible to compare simulated and observed values in a detailed way. However, the good fit in the median values of nutrients and phytoplankton indicated that the model is a promising step towards a management tool for subtropical ecosystems. Acknowledgments This work was supported by the Dutch Research Council of Aquatic Ecology and Water Quality Management Group, Department of Environmental Sciences (Wageningen University). We acknowledge with sincere appreciation the assistance provided by the National Institute of Environmental (IBAMA, ESEC TAIM) and the National Council for Scientific and Technological Development (CNPq, Long Term Ecological Research Program) of Brazil, which fully supported this research financially. We also thank an anonymous reviewer for valuable comments. references Barko, J.W., Smart, R.M., 1981. Comparative influences of light and temperature on the growth and metabolism of selected submersed fresh water macrophytes. Ecol. Monogr. 51, 219–236. Bell, T., Neill, W.E., Schluter, D., 2003. The effect of temporal scale on the outcome of trophic cascade experiments. Oecologia 134, 578–586. Bertolo, A., Lacroix, G., Lescher-Moutoue, F., 1999. Scaling food chains in aquatic mesocosms: do the effects of depth override the effects of planktivory? Oecologia 121, 55–65. Blumberg, A., Mellor, G., 1987. A description of the three-dimensional coastal ocean circulation model. In: Heaps, N. (Ed.), Three-dimensional Coastal Ocean Model. AGU, Washington, DC. Bonnet, M.P., Wessen, K., 2001. ELMO, a 3D water quality model for nutrients and chlorophyll: first application on a lacustrine ecosystem. Ecol. Model. 141, 19–33. 135 Canale, R.P., Vogel, A.H., 1974. Effects of temperature on phytoplankton growth. J. Environ. Eng. Div.-ASCE 100, 231– 241. Carrick, H.J., Aldridge, F.J., Schelske, C.L., 1993. Wind influences phytoplankton biomass and composition in a shallow, Productive Lake. Limnol. Oceanogr. 38, 1179– 1192. Casulli, V., 1990. Semi-implicit finite-difference methods for the 2-dimensional shallow-water equations. J. Comput. Phys. 86, 56–74. Casulli, V., Cheng, R.T., 1990. Stability analysis of Eulerian–Lagrangian methods for the one-dimensional shallow-water equations. Appl. Math. Model. 14, 122– 131. Casulli, V., Cattani, E., 1994. Stability accuracy and efficiency of a semiimplicit method for 3Dimensional shallow-water flow. Comput. Math. Appl. 27, 99–112. Chapra, S.C., 1997. Surface water-quality modeling. McGraw-Hill Series in Water Resources and Environmental Engineering, 844 pp. Chow, V.T., 1959. Open Channel Hydraulics, McGraw-Hill, New York, 680 pp. Daily, J.W., Harlerman, D.R.F., 1966. Fluid Dynamics. AddisonWesley. Dungan, J.L., Perry, J.N., Dale, M.R.T., Legendre, P., Citron-Pousty, S., Fortin, M.J., Jakomulska, A., Miriti, M., Rosenberg, M.S., 2002. A balanced view of scale in spatial statistical analysis. Ecography 25, 626–640. Edmondson, W.T., Lehman, J.T., 1981. The effect of changes in the nutrient income on the condition of Lake Washington. Limnol. Oceanogr. 26, 1–29. Edwards, A.M., Brindley, J., 1999. Zooplankton mortality and the dynamical behaviour of plankton population models. Bull. Math. Biol. 61, 303–339. Eppley, R.W., 1972. Temperature and phytoplankton growth in sea. Fish. Bull. (Washington, DC) 70, 1063–1085. Fasham, M.J.R., Flynn, K.J., Pondaven, P., Anderson, T.R., Boyd, P.W., 2006. Development of a robust marine ecosystem model to predict the role of iron in biogeochemical cycles: a comparison of results for iron-replete and iron-limited areas, and the SOIREE iron-enrichment experiment. Deep-Sea Research. Part I. Oceanographic Research Papers 53, pp. 333–366. Fragoso Jr., C.R., Ferreira, T.F., Motta Marques, D., Collischonn, W., Van Nes, E.H., Scheffer, M., 2007. A complex computational system to cascading trophic interactions evaluation and alternative steady states in subtropical and tropical ecosystems. In: ABRH (Ed.), Proceedings of the 11th International Conference on Diffuse Pollution and the 1st Joint Meeting of the IWA Diffuse Pollution and Urban Drainage Specialist Groups. Minas Gerais, Brazil. Garcia, A.M., Hoeinghaus, D.J., Vieira, J.P., Winemiller, K.O., Marques, D., Bemvenuti, M.A., 2006. Preliminary examination of food web structure of Nicola Lake (Taim Hydrological System, south Brazil) using dual C and N stable isotope analyses. Neotrop. Ichthyol. 4, 279–284. Gross, E.S., Koseff, J.R., Monismith, S.G., 1999a. Evaluation of advective schemes for estuarine salinity simulations. J. Hydraulic Eng.-ASCE 125, 32–46. Gross, E.S., Koseff, J.R., Monismith, S.G., 1999b. Three-dimensional salinity simulations of south San Francisco Bay. J. Hydraulic Eng.-ASCE 125, 1199–1209. Hastings, A., 1990. Spatial heterogeneity and ecological models. Ecology 71, 426–428. Hamilton, D., Schladow, S., Zic, I., 1995. Modelling Artificial Destratification of Prospect and Nepean Reservoirs: Final Report. WP 922 DH, UWA, Centre for Water Research. Hamilton, D.P., Hocking, G.C., Patterson, J., 1995. Criteria for selection of spatial dimension-ality in the application of one 136 e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 and two dimensional water quality models. In: T.U.O.N. Modsim’ 95 (Ed.). Water Res. Ecol. International Congress on Modelling and Simulation. Hölker, F., Breckling, B., 2002. Scales, hierarchies and emergent properties in ecological models: conceptual explanations. In: Hölker, F. (Ed.), Scales, hierarchies and emergent properties in ecological models. Theorie in der Ökologie 6. Peter Lang, Frankfurt, pp. 7–27. Huisman, J., Van Oostveen, P., Weissing, F.J., 1999. Critical depth and critical turbulence: two different mechanisms for the development of phytoplankton blooms. Limnol. Oceanogr. 44, 1781–1787. Imberger, J., 1995. Flux Paths in a Stratified Lake: A Review IU-TAM Symposium on Physical Limonology. Broome, Australia. Imberger, J., Patterson, J.C., 1990. Physical limnology. Adv. Appl. Mech. 27, 303–475. Janse, J.H., 2005. Model Studies on the Eutrophication of Shallow Lakes and Ditches. Wageningen University, Wageningen, 376 pp. Jørgensen, S.E., 1994. Fundamentals of Ecological Modelling. Developments in Environmental Modelling, 19, 2nd edition. Elsevier, Amsterdam, p. 628. Jørgensen, S.E., Fath, B.D., Grant, W.E., Legovic, T., Nielsen, S.N., 2008. New initiative for thematic issues: an invitation. Ecol. Model. 215, 273–275. Kamenir, Y., Dubinsky, Z., Alster, A., Zohary, T., 2007. Stable patterns in size structure of a phytoplankton species of Lake Kinneret. Hydrobiologia 578, 79–86. Levin, S.A., 1992. The problem of pattern and scale in ecology. Ecology 73, 1943–1967. Lord, D., Imberger, J., Pattiaratchi, C., 1994. Management of coastal waters in Western Australia, the use of integrated models. In: Yokosuka, J. (Ed.), Internationial Conference on Hydro-technical Engineering for Port and Harbour Construction. Lucas, L.V., 1997. A Numerical Investigation of Coupled Hydrodynamics and Phytoplankton Dynamics in Shallow Estuaries. University of Stanford. Lucas, L.V., Cloern, J.E., 2002. Effects of tidal shallowing and deepening on phytoplankton production dynamics: a modeling study. Estuaries 25, 497–507. Matveev, V.F., Matveeva, L.K., 2005. Seasonal succession and long-term stability of a pelagic community in a productive reservoir. Mar. Freshwater Res. 56, 1137–1149. Mehner, T., Holker, F., Kasprzak, P., 2005. Spatial and temporal heterogeneity of trophic variables in a deep lake as reflected by repeated singular samplings. Oikos 108, 401–409. Mitra, A., Flynn, K.J., 2007. Importance of interactions between food quality, quantity, and gut transit time on consumer feeding, growth, and trophic dynamics. Am. Nat. 169, 632–646. Mitra, A., Flynn, K.J., Fasham, M.J.R., 2007. Accounting for grazing dynamics in nitrogen–phytoplankton–zooplankton models. Limnol. Oceanogr. 52, 649–661. Moss, B., 1998. Shallow lakes: biomanipulation and eutrophication. Scope Newlett. 29, 45. Mukhopadhyay, B., Bhattacharyya, R., 2006. Modelling phytoplankton allelopathy in a nutrient-plankton model with spatial heterogeneity. Ecol. Model. 198, 163–173. Olsen, P., Willen, E., 1980. Phytoplankton Response to Sewage Reduction in Vattern, a Large Oligotrophic Lake in Central Sweden. Arch. Hydrobiol. 89, 171–188. Platt, T., Dickie, L.M., Trites, R.W., 1970. Spatial Heterogeneity of Phytoplankton in a near-Shore Environment. J. Fish. Res. Board Can. 27, 1453–1465. Rajar, R., Cetina, M., 1997. Hydrodynamic and water quality modelling: an experience. Ecol. Model. 101, 195–207. Reynolds, C.S., 1999. Modelling phytoplankton dynamics and its application to lake management. Hydrobiologia 396, 123–131. Riley, G.A., 1956. Oceanography of Long Island Sound, 1952–1954. II. Physical oceanography, Bulletin of the Bingham Oceanographic Collection XV, 15–46. Robson, B.J., Hamilton, D.P., 2004. Three-dimensional modelling of a Microcystis bloom event in the Swan River estuary, Western Australia. Ecol. Model. 174, 203–222. Romero, J.R., Antenucci, J.P., Imberger, J., 2004a. One- and three-dimensional biogeochemical simulations of two differing reservoirs. Ecol. Model. 174, 143– 160. Romero, J.R., Hipsey, M.R., Antenucci, J.P., Hamilton, D., 2004b. Computational Aquatic Ecosystem Dynamics Model: CAEDYM v2. 1 Science Manual. Centre for Water Research, University of Western Australia, Nedlands, WA, Australia. Sas, H., 1989. Lake Restoration by Reduction of Nutrient Loading: Expectations, Experiences, Extrapolations. Academia Verlag Richarz, St. Augustin, 1–497 pp. Scheffer, M., 1998. Ecology of Shallow Lakes. Population and Community Biology. Chapman and Hall, London, 0– 357 pp. Scheffer, M., De Boer, R.J., 1995. Implications of spatial heterogeneity for the paradox of enrichment. Ecology 76, 2270–2277. Scheffer, M., Bakema, A.H., Wortelboer, F.G., 1993. MEGAPLANT—a simulation model of the dynamics of submerged plants. Aquat. Bot. 45, 341–356. Scheffer, M., Straile, D., Van Nes, E.H., Hosper, H., 2001. Climatic warming causes regime shifts in lake food webs. Limnol. Oceanogr. 46, 1780–1783. Schindler, D.W., 1975. Modelling the eutrophication process. J. Fish. Res. Board Can. 32, 1673–1674. Schladow, S.G., Hamilton, D.P., 1997. Prediction of water quality in lakes and reservoirs. 2. Model calibration, sensitivity analysis and application. Ecol. Model. 96, 111–123. Smith, R.A., 1980. The theoretical basis for estimating phytoplankton production and specific growth-rate from chlorophyll, light and temperature data. Ecol. Model. 10, 243–264. Spencer, D.F., 1986. Early growth of Potamogeton pectinatus L. in response to temperature and irradiance: morphology and pigment composition. Aquat. Bot. 26, 1–8. Spillman, C.M., Imberger, J., Hamilton, D.P., Hipsey, M.R., Romero, J.R., 2007. Modelling the effects of Po River discharge, internal nutrient cycling and hydrodynamics on biogeochemistry of the Northern Adriatic Sea. J. Mar. Syst. 68, 167–200. Steele, J.H., 1965. Notes on some theoretical problems in production ecology. In: Goldman, C.R. (Ed.), Primary Production in Aquatic Environments. University of California Press, Berkeley, CA. Steele, J.H., 1978. Spatial Pattern in Plankton Communities. Publisher, Plenum Press, New York, 470 pp. Steele, J.H., Henderson, E.W., 1992. A simple model for plankton patchiness. J. Plankton Res. 14, 1397– 1403. Thoman, R.V., Segna, J.S., 1980. Dynamic phytoplankton-phosphorus model of Lake Ontario: ten-year verification and simulations. In: Loehr, C., Martin, C.S., Rast, W. (Eds.), Phosphorus Management Strategies for Lakes. Ann Arbor Science Publishers, Ann Arbor, MI, pp. 153– 190. Van den Berg, M.S., Coops, H., Meijer, M.L., Scheffer, M., Simons, J., 1998. Clear water associated with a dense Char a vegetation in the shallow and turbid Lake Veluwemeer, the Netherlands. In: Jeppesen, E., Søndergaard, M., Søndergaard, M., Kristoffersen, K. (Eds.), Structuring Role of Submerged Macrophytes in Lakes. Springer-Verlag, New York, pp. 339–352. Van Leeuwen, E., Lacerot, G., Van Nes, E.H., Hemerik, L., Scheffer, M., 2007. Reduced top–down control of phytoplankton in e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 125–137 warmer climates can be explained by continuous fish reproduction. Ecol. Model. 206, 205–212. Van Nes, E.H., Scheffer, M., 2005. A strategy to improve the contribution of complex simulation models to ecological theory. Ecol. Model. 185, 153–164. Wetzel, R.G., 1996. Limnology. W.B. Saunders Co., Philadelphia. White, F.M., 1974. Viscous Fluid Flow, McGraw-Hill, New York, 32 pp. 137 Woodward, G., Hildrew, A.G., 2002. Food web structure in riverine landscapes. Freshwater Biol. 47, 777–798. Wu, F.C., Shen, H.W., Chou, Y.J., 1999. Variation of roughness coefficients for unsubmerged and submerged vegetation. J. Hydraulic Eng.-ASCE 125, 934–942. Wu, J., 1982. Wind-stress coefficients over sea-surface from breeze to hurricane. J. Geophys. Res.-Oceans Atmos. 87, 9704–9706.