Application of a 3D hydrodynamic–biological model for seasonal and spatial dynamics of water quality and phytoplankton in Lake Erie

Journal of Great Lakes Research 37 (2011) 41–53

Contents lists available at ScienceDirect

Journal of Great Lakes Research j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / j g l r

Application of a 3D hydrodynamic–biological model for seasonal and spatial dynamics of water quality and phytoplankton in Lake Erie Luis F. Leon a,f, Ralph E.H. Smith a,⁎, Matthew R. Hipsey b,c, Serghei A. Bocaniov a, Scott N. Higgins a,d, Robert E. Hecky a,e, Jason P. Antenucci b, Jorg A. Imberger b, Stephanie J. Guildford a,e a

University of Waterloo, Department of Biology, Waterloo, Ont., Canada N2L 3G1 University of Western Australia, Centre for Water Research, Crawley, WA, 6009, Australia University of Western Australia, School of Earth and Environment, Crawley, WA, 6009, Australia d University of Wisconsin, Madison Center for Limnology, Madison, WI 53706, USA e University of Minnesota, Large Lakes Observatory, Duluth, MN 55812, USA f National Water Research Institute, Environment Canada, Burlington, Ont., Canada L7R 4A6 b c

a r t i c l e

i n f o

Article history: Received 10 December 2009 Accepted 14 September 2010 Available online 28 January 2011 Communicated by Ram Yerubandi Index words: Lake Erie Three dimensional modeling Nutrients Phytoplankton

a b s t r a c t In large lakes, temporal variability is compounded by strong spatial variability associated with mesoscale physical processes such as upwelling and basin-scale circulation. Here we explore the ability of a three dimensional model (ELCOM–CAEDYM) to capture temporal and spatial variability of phytoplankton and nutrients in Lake Erie. We emphasized the east basin of the lake, where an invasion by dreissenid mussels has given special importance to the question of spatial (particularly nearshore–offshore) variability and many comparative observations were available. We found that the model, which did not include any simulation of the mussels or of smaller diffuse nutrient sources, could capture the major features of the temperature, nutrient and phytoplankton variations. Within basin variability was large compared to among-basin variability, especially but not exclusively in the western regions. Consistent with observations in years prior to, but not after, the mussel invasion the model predicted generally higher phytoplankton concentrations in the nearshore than the offshore zones. The results suggest that the elevated phytoplankton abundance commonly observed in the nearshore of large lakes in the absence of dreissenid mussels does not have to depend on localized nutrient inputs but can be explained by the favourable light, temperature and nutrient environment in the shallower and energetic nearshore zone. The model is currently being extended to allow simulation of the effects of dreissenid mussels. © 2010 International Association for Great Lakes Research. Published by Elsevier B.V. All rights reserved.

Introduction Large lakes, such as the Laurentian Great Lakes, share many physical similarities with coastal marine ecosystems (Rao and Schwab, 2007). Physical processes contribute to differentiation of water masses between nearshore and offshore regions while upwelling and other manifestations of episodic energy inputs transport and mix materials and organisms on a variety of time scales. In combination with often localized inputs of nutrients and suspended materials from the catchment, the physical processes generate pronounced variations of chemical and biological properties in space and time (e.g. Budd et al., 2001; Millie et al., 2002). The resulting variability is a challenge in testing and predicting biological and chemical responses to important influences such as nutrient loading or altered food web structure (e.g. Charlton et al., 1999).

⁎ Corresponding author. Tel.: + 1 519 888 4567. E-mail address: [email protected] (R.E.H. Smith).

Lake Erie is a good example of the complexities created by the physically-driven, three-dimensional variability of large lakes as they interact with major anthropogenic and other perturbations. It has a history of eutrophication and subsequent reversal based on control of external P loading (Charlton et al., 1999). More recently it has undergone major food web alterations, notably due to colonization by dreissenid mussels but also through establishment of non-native zooplankton and fish. Some symptoms of eutrophication, notably hypolimnetic hypoxia, have not continued the expected course of improvement, or may even have worsened, but trends and causes are still debatable (e.g. Conroy et al., 2005) given the complex and nonlinear interactions that occur within the lake ecosystem. The lake has pronounced large scale patterns in morphometry and external nutrient inputs. Most of the external nutrient loading enters the shallow and relatively well-mixed western basin (Dolan and McGunagle, 2005) which is also frequently the site of problem blooms of phytoplankton. The large central basin is deep enough to form a seasonal thermocline but the hypolimnion is thin and prone to hypoxia while nutrient and plankton concentrations are highly variable, particularly in the area adjacent to the west basin (Charlton

0380-1330/$ – see front matter © 2010 International Association for Great Lakes Research. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.jglr.2010.12.007

42

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

et al., 1999; Guildford et al., 2005; Smith et al., 2005). The east basin is the point of maximum depth for the lake but also features extensive shallow nearshore areas. It has been proposed that the potentially high impact of the benthic dreissenid mussels in the shallower nearshore waters has re-configured the spatial distribution of phytoplankton and nutrients in the east basin and even promoted problematic growth of benthic algae (Hecky et al., 2004, Higgins et al., 2008). The management of hypolimnetic hypoxia, nuisance blooms of cyanobacteria, dreissenid mussel impacts, and other issues of concern requires an understanding of physical, chemical, and biological processes all of which are dynamically coupled through space and time. The application of numerical and computer simulation models that can account for these dynamic features would be of immense utility for management and ecosystem theory. Quantification and prediction of invasive mussel impacts, for example, really depends on a credible representation of the physical environment that dictates their access to particles in the water column (e.g. Boegman et al., 2008). As another example, observations of seasonal oxygen depletion in the central basin often reveal strong spatial and temporal variations in its progress related to variable advection, mixing and internal wave activity (Rao et al., 2008). Models that can account for such physical processes together with key ecological processes are an important tool for elucidating the mechanisms involved and for designing appropriate monitoring and management programs. Ecosystem modeling in large lakes, and specifically the Laurentian Great Lakes, has a long history. Earlier models were highly aggregated and did not aim to capture mixing, exchange and circulation processes in a mechanistic manner, and would not capture within-basin processes at all. More recently, two dimensional models (e.g. Boegman et al., 2008; Zhang et al., 2008) have been used to analyze chemical and biological, as well as physical, dynamics. Three dimensional models have been used to simulate and analyze physical properties (e.g. Beletsky and Schwab, 2001; Leon et al., 2005), study advection of larval fish (Beletsky et al., 2007) and generate onedimensional portrayals of nutrient and plankton dynamics (Chen et al., 2002). Published examples of three dimensional modeling of physical, chemical and biological processes appear to be limited so far to Leon et al. (2006), who demonstrated the ability of a coupled hydrodynamic–biological model (ELCOM–CAEDYM) to capture the major features of oxygen and phosphorus dynamics in Lake Erie, and Schwab et al. (2009), who used vertically averaged dynamics from the Princeton Ocean Model to analyze spatial and temporal patterns of phosphorus in Lake Erie. Here we report on further application and assessment of ELCOM– CAEDYM (ELCD for short) for modeling nutrients and phytoplankton in Lake Erie. ELCD has been successfully applied to model water quality, plankton, and benthic–pelagic interaction processes in a variety of aquatic systems (e.g. Romero et al., 2004; Robson and Hamilton, 2004; Spillman et al., 2008). Our objectives here were, first, to use the extensive data available for Lake Erie in 2002 to make the first assessment of how well ELCD could represent highly dynamic physical (e.g. temperature, circulation), chemical (standard nutrients) and biological (phytoplankton biomass) properties within one of the world's largest lakes. Second, we utilize the model to illustrate how physical and ecological processes may be expected to influence the spatial dynamics of phytoplankton production. Third, we utilize the model to simulate the expected zonation of physical, chemical, and biological properties across the nearshore–offshore gradient and compare these estimates to measured values pre- and post-dreissenid invasion. Since our model does not currently include dreissenid filtration, our hypothesis was that the model should provide reasonable representations of the nearshore–offshore gradients as they existed during pre-invasion, but not post-invasion, years. Differences in chemical and biological values between modeled and measured values during post-invasion would then provide an

estimate of dreissenid impacts for regions of the lake where predreissenid data are not available. Materials and methods Model structure The three-dimensional (3D) Estuary, Lake and Coastal Ocean Model (ELCOM) solves the unsteady Reynolds–averaged Navier– Stokes equations for heat and momentum transfer across the water surface due to wind and atmospheric thermodynamics in order to simulate the spatial and temporal variations of the physical and transport processes in the water body. It has been described and its performance analyzed in previous publications (e.g. Hodges et al., 2000; Leon et al., 2005; Morillo et al., 2008), to which the reader is referred for further information. The Computational Aquatic Ecosystem Dynamics Model (CAEDYM) interfaces with ELCOM and is documented and described in detail elsewhere (e.g. Hipsey et al., 2006). For this application it was used to simulate the C, N, P, DO, and Si cycles, along with inorganic suspended solids and phytoplankton (App. A). Each simulated variable is subject to advection, mixing and boundary forcing by ELCOM. The process equations and rationale have been presented in full detail in previous publications (e.g. Robson and Hamilton, 2004; Romero et al., 2004; Spillman et al., 2008; Gal et al., 2009) and the complete model structure and operation is explained in online documents freely available from the Centre for Water Research, University of Western Australia. Twelve state variables are required to model the algal biomass, including the dynamically calculated internal nutrient stores of N and P (App. A, B). Five dissolved inorganic nutrients (PO4, NO3, NH4, DIC, and RSi), three dissolved organic (DOC, DON, and DOP) and three particulate detrital organic matter groups (POC, PON, and POP), two inorganic suspended solids size classes (SS1 and SS2), and dissolved oxygen (DO) were simulated. For comparison with data, total nitrogen (TN) and total phosphorus (TP) were derived by totaling the relevant N and P state variables. The model simulates PO4 but operational measurements are conventionally termed soluble reactive P (SRP) and we use this latter terminology in subsequent sections. For primary production, the shortwave radiation (280–2800 nm) intensity at the surface is converted to the photosynthetically active component (PAR) assuming that 45% of the incident spectrum lies between 400 and 700 nm. PAR is assumed to penetrate into the water column according to the Beer–Lambert Law (App. A, Eq. (1)) with the light extinction coefficient dynamically adjusted to account for variability in the concentrations of algal, inorganic and detrital particulates, and dissolved organic carbon levels. Two inorganic particle groups (SS) were included within the simulation (App. A, Eq. (2)) and modeled as a balance between settling (based on Stokes Law) and resuspension. At the sediments, CAEDYM maintains mass balance of C, N, P, Si, DO and SS in both the water column and a single sediment layer. The sediment fluxes of dissolved inorganic and organic nutrients are based on empirical formulations that account for sensitivities to temperature and oxygen (as a proxy for the sediment redox condition). The DO balance includes atmospheric exchange, sediment oxygen demand, decomposition and nitrification demands, photosynthesis and respiration (App. A, Eq. (3)). Sediment oxygen demand (SOD) is a function of water temperature and dissolved oxygen levels but at 20 °C and non-limiting dissolved oxygen concentrations would be 0.6 gO2 m− 2 d− 1, within the range of published SOD estimates for central basin Lake Erie (e.g. Matisoff and Neeson, 2005). Dissolved oxygen dynamics are the subject of forthcoming work but the model appears to generate realistic predictions (Leon et al., 2006) to support the oxygen-dependent nutrient cycling and biological process descriptions that are the focus of this paper.

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

Both the inorganic and organic, as well as dissolved and particulate forms of C, N and P are modeled explicitly along the degradation pathway of POM to DOM to dissolved inorganic matter, but only the balance equations for inorganic forms are shown here (App. A, Eqs. (5)–(7)). The nitrogen cycle includes the processes of denitrification, nitrification and N2 fixation. Decomposition follows first-order kinetics, dependent on substrate concentrations, temperature and DO. The silica cycle is simpler and includes uptake of dissolved Si (RSi) by diatoms into the internal Si (ISi) pool, dissolved sediment fluxes of RSi, diatom mortality directly into the RSi sediment pool, settling of ISi, and resuspension of ISi. This relatively simple representation assumes that diatom frustules rapidly settle into the sediments and the organic silica is rapidly mineralized. The DIC balance (App. A, Eq. (5)) includes atmospheric fluxes of CO2 that are based on the difference between the atmospheric and water column values of pCO2. Phytoplankton description While the hydrodynamics model (ELCOM) has been previously applied in Lake Erie (Leon et al., 2005), the coupled hydrodynamics– biological model (ELCOM–CAEDYM, or ELCD) has not previously been applied to Lake Erie or other large temperate lakes. A key issue in the current application was the configuration of the phytoplankton groups. The model allows for up to 7 different groups defined, as is common in similar ecosystem models (e.g. Bierman and Dolan, 1981; Chen et al., 2002; Arhonditis and Brett, 2005; Mieleitner and Reichert, 2008) by a mixture of taxonomic and functional characteristics. We chose to use five phytoplankton groups with plausible characteristics relative to the major groups known to be important in the lake (e.g. Makarewicz et al., 1999; Barbiero and Tuchman, 2001; Smith et al., 2007). The Early Diatoms group represents early-blooming diatom taxa with high Si requirements and rapid sinking rates. The Late Diatoms group represents diatoms with lower Si requirements and sinking rates that occur later in the seasonal succession; this group could include some of the silicified chrysophyceans as well. A Cyanophyte group was defined to represent mainly the larger and potentially N-fixing taxa that are often associated with warm and stable waters. The Flagellates group was defined to represent cryptophytes and other flagellates typical of cooler waters and/or deeper strata; this group could include some nonmotile forms that have similar preferences and low sinking rates (e.g. some picocyano-

43

bacteria). Finally, the Other Phytoplankton group was defined to represent flagellates and nonmotile forms that are typical of warmer and brighter conditions (e.g. many Chlorophytes and some dinoflagellates, chrysophyceans and haptophytes). Phytoplankton dynamics account for variations in potential growth rates (gross photosynthesis), losses to respiration, excretion and mortality, and the difference between settling and resuspension rates (App. A, Eq. (4)). Growth rate is determined by a maximum potential rate at 20 °C as modified by a temperature response function and the minimum value of expressions for limitation by light, P, N and Si (App. B). The temperature function for growth, f T1, employs constants selected to produce maximum productivity at a temperature TOPT and zero growth at some higher temperature TMAX. Below the temperature (TSTD) the productivity follows a simple exponential relation of the same form as for algal loss rates, Ri (App. B). The internal P and N within the algal groups are modeled as dynamic intracellular stores, and these are able to regulate growth through f(N) and f(P). N and P uptake rates are functions of external and internal nutrient concentrations and algae are given a preference for ammonium over nitrate (App. B). Losses to respiration, natural mortality and excretion (Ri) follow Arrhenius-based, temperature-dependent, scaling (App. B). For each group, respiration comprises a constant fraction (fres) of the total losses. A constant fraction (fdom) of the non-respiratory losses goes to dissolved compounds (excretion) and the balance to particulate detritus. Metabolic loss of nutrients due to mortality and exudation is proportional to the dynamic internal nutrient ratio multiplied by the loss rate. Losses to sedimentation were represented by a constant settling rate, adjusted for water viscosity, specific to each phytoplankton group. Phytoplankton parameters are summarized in App. C. Maximum specific growth rate parameters were set in a similar range to some previous models of the Great Lakes. The cardinal temperatures defining the temperature responses of the groups were chosen to reflect the seasonal occurrence of related groups in Lake Erie (e.g. Makarewicz et al., 1999; Carrick, 2004; Smith et al., 2007). The temperature multiplier for growth mainly determines the steepness of response below TSTA and above TOPT and was similar among groups. Ik (marking the upper limit of the light-limited response region) was set to resemble the measured properties of Great Lakes communities (Fahnenstiel et al., 1989; Smith et al., 2005 and unpub. data; Depew et al., 2006) at times when the different groups would normally be most

Table 1 Initialization values of chemical parameters in early spring for the lake-wide simulation in 2002. Basin

Station

Depth (m)

TP (mgP m− 3)

SRP (mgP m− 3)

NH4 (mgN m− 3)

NO3 (mgN L− 1)

SRSi (mgSi m− 3)

Chl-a (mg m− 3)

East

879 886 934 940 945 946 949 952 953 954 961 962 965 966 968 969 970 971 973 974

62 50 28 49 21 23 22 22 23 24 21 19 14 11 10 8 10 9 6 8

12.9 13.2 12.6 13.7 14.7 12.6 16.6 16.1 16.6 18.3 13.3 22.8 21.5 18.4 21.6 24.8 18.4 18.4 28.0 18.4

8.3 8.5 8.1 8.9 9.6 8.1 10.9 10.6 10.9 12.1 8.6 15.3 14.4 5.1 6.0 6.9 5.1 5.1 7.8 5.1

22 16 19 19 16 11 19 10 9 8 2 35 31 50 38 90 13 10 90 90

300 284 305 307 197 257 221 188 248 225 100 1080 1110 906 1066 1720 906 906 1390 906

870 870 960 940 820 900 890 1200 920 1200 570 2240 2120 1630 1790 2300 1090 1700 2600 1440

0.1 0.3 0.7 0.1 0.3 0.2 0.1 0.4 0.4 0.4 2.2 0.7 0.8 1.5 1.2 0.4 0.9 0.4 0.8 0.8

Central

West

44

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

important. The assigned values were also consistent with Reynolds (1984). The C:Chl a ratio for most groups was made similar to values suggested by comparisons of POC and Chl a in Lake Erie (Guildford et al., 2005; Smith et al., 2007 and Smith unpub. data), except that a higher value was used for Flagellates as suggested by Istvanovics et al. (1994). Nitrogen utilization kinetics were based on algal culture studies (App. C) and on previous successful applications of EC elsewhere (Robson and Hamilton, 2004; Romero et al., 2004; Gal et al., 2009). Compared to a relevant review of phytoplankton nutrient kinetics (Lehman et al., 1975) our parameters differed less among groups but, depending on the exact assumptions made to interconvert units, would support similar rates of N uptake and N-limited growth under most conditions. Since 1995, inorganic N is usually N0.2 mg L− 1 in Lake Erie (Charlton et al., 1999 and Smith unpub. data) so we did not expect the N kinetics to have a strong influence on phytoplankton dynamics. For P, half saturation constants were made the same for all groups except Cyanobacteria, which were assigned a somewhat higher value consistent with literature reports (App. C). Minimum and maximum intracellular P quotas were similar to previous successful applications of ELCD (e.g. Romero et al., 2004) and were comparable to values in several algal studies (App. C), again allowing for uncertainty introduced through unit conversions. Maximum P uptake rates were chosen to reflect the range suggested by Smith and Kalff (1983) and Lehman et al. (1975). Half-saturation constants and constant internal Si quotas for diatoms were representative of the lower end of the range reported in studies of marine and freshwater diatoms (App. C), with Late Diatoms assumed to be more efficient and more lightly silicified than Early Diatoms. With a constant internal Si quota, the maximum Si uptake rate scales linearly with growth rate and is not specified as an additional parameter, but would take values as high as 60 (Late Diatoms) to 110 mgSi mgChla− 1 d− 1 (Early Diatoms). These values are comparatively high and, together with the KSi and kiSi values make the model diatoms quite efficient at Si utilization. Settling rates specified for each group were comparable to several previous modeling studies (App. C) and within the considerable range reported in the literature. Other loss rates (respiration, excretion and mortality) could vary from about 5 to 50% of biomass per day among groups and with changing temperature, given the chosen parameters. Under most conditions, losses would be in the range 6 to 25% of biomass per day. Respiration alone is mostly in the range 5–15% of

light-saturated photosynthesis, a plausible range considering existing knowledge of dark algal respiration (Geider and Osborne, 1989). These ranges are similar to some previous Great Lakes models that, like ours, lump respiration, mortality and all other losses except sinking and advection together (e.g. Bierman and Dolan, 1981). We assumed that Flagellates and Others would be mostly edible to zooplankton and therefore have higher total loss rates (kra) than other groups and that POM and DOM production (notionally related to grazing as well as other mortality processes) would comprise a larger share of the total losses compared to respiration alone (i.e., lower fRes). Model setup and validation The bathymetry was set up as described in Leon et al. (2005) on a 2 km grid. The meteorological data (air temperatures, wind speed/ direction, solar radiation, relative humidity and cloud cover) for 2002 were obtained from the archives of the National Water Research Institute (NWRI) and the NOAA National Data Buoy Center. Different values for the forcing data were used to take into account the spatial variability of meteorological conditions across the lake. The lake was divided in three sections (west, central and east basins, Fig. 1) where hourly input was provided. App. D shows the hourly meteorological data used in the 2002 simulation. Model runs used 40 vertical layers, with variable thickness across the depth, and a time step of 5 min was specified bearing in mind numerical stability (CFL) constraints. To simulate the whole lake for 190 days (mid-April to mid-October) required 112 h of computing time on a standard desk-top computer and values for state variables were stored as 1 hour averages. Data on in-lake concentrations of phytoplankton Chl a, light attenuation, and water chemistry for model initialization and validation were taken from an Environment Canada–University of Waterloo study that sampled nearshore and offshore stations, mainly in the east basin, monthly from April to September (Depew et al., 2006, Smith unpub. data). Many of these data also reside in the ‘Star Database’, housed by the NWRI (National Water Research Institute) of Environment Canada. Observations from the U.S. E.P.A. spring survey cruise (Barbiero and Tuchman, 2001 describe stations and methods) provided the additional data necessary to describe the initial conditions in the lake at the beginning of the simulation period. Additional observations, mostly for central and west basins, were made available from the U.S. E.P.A. Lake Erie Trophic Status project (Matisoff and Ciborowski, 2005), and provided data in June, July,

Fig. 1. Map of Lake Erie showing the stations used for lake wide initialization (Table 1) plus (east basin only) stations used for simulated nearshore–offshore comparison study (App. E).

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

August and September. Values for major state variables (Table 1) at stations sampled in spring-time (April) were used to generate concentration contours and initial values for each model cell, assuming vertical homogeneity at the beginning of the simulation period (Apr. 11, 2002). Fig. 1 shows the initialization stations plus, in the east basin, additional stations that were used in a simulated nearshore–offshore study (see Results). Validation comparisons used stations shown in Fig. 1 plus some additional ones in central and west basins, depending on availability of data in any particular time period. App. E gives a detailed list of stations used for the various comparisons we present. The model includes eleven inflows (Table 2, Fig. 1) and one outflow (Niagara River). Daily flow rates for the Detroit, Niagara and Grand (ON) Rivers were based on either NOAA (Great Lakes Research Laboratory) or Water Survey of Canada (Environment Canada) databases. Daily flows for all Ohio tributaries but the Rocky River were obtained from the database maintained by the Heidelberg College Water Quality Laboratory. For the Rocky River, the daily discharge was assumed to be 0.35 times of that of the Cuyahoga River where 0.35 is the ratio of average annual flows of two rivers (Bolsenga and Herdendorf, 1993; Schwab et al., 2009). Daily flows for Cattaraugus and Buffalo Rivers we estimated based on both U.S. Geological Survey (USGS) data and Great Lakes Environmental Research Laboratory's (GLERL) Large Basin Runoff Model (LBRM) results. The Niagara River is responsible for 88% of the water outflows from the lake while the other sources of outflow are much smaller, e.g. Welland Diversion (3%), evaporation (8%) and consumptive use (1%) (Bolsenga and Herdendorf, 1993). The eleven inflows included here (Table 2) account for 97.5% of all lake inflows (Bolsenga and Herdendorf, 1993). Modeled nutrient loads to Lake Erie included only those from the eleven chosen inflows but these account for almost 87% of TP loads to the lake based on recent assessments (Dolan and McGunagle, 2005; Schwab et al., 2009). This modest underestimation of nutrient loadings should not have a large effect on our results because lake residence time is quite long (2.7 years) and the simulation period (Apr. 11 to Oct. 20) starts after the loading from the smaller tributaries has passed its annual maximum. The in-lake concentrations at the beginning of the simulation thus capture much of the effect of the loading that we do not explicitly simulate. To estimate loadings from the tributaries, we used a combination of government databases and published sources (e.g. Dolan and McGunagle, 2005). To estimate loadings from the Detroit River we used data retrieved from STORET, a U.S. EPA database of water quality monitoring data. The selected station was near Rockwood (station ID # 820017). Daily loadings for each of the Ohio tributaries was calculated as the sum of the tributary load retrieved from the Heidelberg College National Center for Water Quality Research

Table 2 Loadings of suspended solids (SS), total P (TP) and total N (N) from Lake Erie tributaries for the period April 11–October 20, 2002, together with the cumulative flows.

1 2 3 4 5 6 7 8 9 10 11

River

Flow (km3)

SS (Mt)

TP (Mt)

TN (Mt)

Detroit R. Maumee R. Grand R. (Ontario) Sandusky R. Cuyahoga R. Raisin R. Buffalo R. Cattaraugus R. Grand R. (Ohio) Rocky R. Vermilion R. Total

86.64 1.54 0.76 0.39 0.35 0.30 0.30 0.28 0.23 0.12 0.07 90.98

685,000 251,000 38,000 95,000 106,000 34,000 15,000 29,000 30,000 24,000 11,000 1,318,000

1383 592 109 200 148 109 32 60 61 33 19 2746

46,652 11,044 1932 3476 2752 1708 832 1552 581 615 300 71,444

45

database and the load from the direct and indirect industrial and municipal discharges located downstream of the river monitoring site. In the Rocky River, concentrations of the suspended solids, phosphorus and nitrogen were assumed to be 0.65 times of those in the Cuyahoga River where 0.65 represents the ratio of average annual phosphorus concentrations in both rivers (Schwab et al., 2009). Loadings for the Grand River in Ontario were estimated based on the data obtained from the Ontario Ministry of the Environment and Grand River Conservation Authority (Ontario). For the Cattaraugus and Buffalo Rivers, the following approach was used. First, seasonal (April 11 to October 20) total phosphorus (TP) loads were estimated similar to the method described in Schwab et al., 2009. Then, the average TP:TN and TP:SS ratios calculated for all other tributaries (=9 tributaries) with known loads of TP, total nitrogen (TN) and suspended solids (SS) were used to estimate seasonal loads of TN and SS. Finally, daily flows were used to apportion the seasonal loads of TP, TN and SS into daily loads. The model requires estimates for different fractions of the TP and TN loads, just as for in-lake concentrations (Table 1). The TP and TN loads were so-apportioned using the average proportions of the components (PO4, NO3, etc) in the total load as determined from available tributary measurements. Estimated loads of water, suspended solids, total phosphorus and nitrogen are summarized in Table 2. Results Comparisons between model and observations The performance of ELCOM in capturing the variations of temperatures and currents in Lake Erie cannot be fully addressed here, but App. F gives one example of how simulated and observed temperatures compared at the deepest point of the east basin (station 23, Fig. 1). Compared to observations, the seasonal thermocline in the simulations was somewhat more diffuse but the onset and duration of stratification, as well as thickness of the surface mixed layer, was similar between observations and predictions. Surface and deep temperatures were also reproduced quite well. Further analysis of ELCOM's simulation of the observed physical properties of Lake Erie is available elsewhere (Leon et al., 2005) and shows that the model captures the major features of temperature distributions as well as current directions and velocities in nearshore east basin, Lake Erie. Work in other lakes (e.g. Lake Kinneret, Gal et al., 2009) suggests the model also captures circulation patterns with reasonable accuracy but formal tests of this for Lake Erie are objects of ongoing work. ELCD predictions for phytoplankton, nutrients and associated variables were compared to observations across space and time. Model predictions for total epilimnetic phytoplankton biomass (expressed as Chl a) in early August, for example, were in the same range as the observations from sampling expeditions in the same time period (Fig. 2). The spatial plot is specifically for model output for Aug. 5, but observations were collected over a period of days (Aug. 5–8). In making X–Y comparisons (inset, Fig. 2) we used model output for the actual day of observation, averaged over a 24 hour period. Modeled and observed concentrations were significantly (p b 0.01) correlated and both showed elevated concentrations in some areas of the west basin. For total phosphorus (TP) the modeled and observed values were again in the same range and significantly (p b 0.01) correlated in the early August period (Fig. 3). Elevated TP concentrations in west basin were evident in both simulations and observations. Using the available measurements of epilimnetic Chl a and TP concentrations in May, June, July, August and September we found that observations could differ markedly from model predictions for individual stations and dates (Fig. 4a, b). There are many possible reasons for such disparities, especially when the observational data have not all been collected by one agency with completely standard

46

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

Fig. 2. Lake wide horizontal distribution of predicted epilimnetic chlorophyll-a concentrations (mg m− 3) with observed values for Aug 5–8, 2002.

methodology. Predictions and observations were nonetheless significantly correlated overall, as seen best in log–log plots (Fig. 4c, R2 = 0.48, p b 0.0001, n = 277; and 4 d, R2 = 0.40, p b 0.0001, n = 208). The model tended to predict TP higher than the observations when observed TP was low (4 to 8 mgP m− 3). The cause for this bias is not known but may point to a need for more information about forms of P entering the lake (see Discussion). There was no evidence of bias in the relationship between predicted and observed Chl a (Fig. 4c). A similar analysis of predicted and observed soluble reactive silicate concentrations (not shown) indicated a highly significant correlation

across the same body of measurements (R2 = 0.71, p b 0.0001, n = 92) and no apparent bias. Time series comparisons between simulations and observations were made for several stations in east basin, where sampling was relatively frequent in 2002, and gave very similar results in terms of the predicted dynamics and correspondence with available observations. At a station close to the deepest point of east basin (station 938, Fig. 1), for example, the simulated Chl a in the top 5 m was similar to observations (Fig. 5). Observations with a calibrated fluorometer gave better agreement with the simulated spring peak than did

Fig. 3. Lake wide horizontal distribution of predicted epilimnetic TP concentrations (expressed here as mg L− 1) with observed values for August 5–8, 2002.

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

47

Fig. 4. Scatter plot of the observed vs. predicted epilimnetic (top 5 m) concentrations of chlorophyll-a (a, c) and total phosphorus (b, d) in all three basins; (a) and (b) show nontransformed data, while (c) and (d) show logarithmically transformed data.

measurements of extracted Chl a. Chl a profiles measured with the fluorometer (not shown) revealed considerable vertical structure and localized maxima that the sampling for extracted Chl a did not always capture. Among the other stations sampled, there were some observations of extracted Chl a as high as the model predictions for the spring peak but both model and observations suggested that very high values (N4 μg L− 1) were quite ephemeral and spatially localized, and would easily be missed when sampling. Time series comparisons for other important properties showed a variable degree of agreement between simulated and observed values (Fig. 5). Observations of TP, for example, showed less variable values over time than did simulations, while total dissolved P (TDP) observations showed more temporal variation than the simulations. Nitrate, ammonium and vertical attenuation coefficient (Kd) showed comparatively little variation over time in both simulations and observations, while soluble reactive phosphorus (SRP) and soluble reactive silica (SRSi) showed strong dynamics and good agreement between simulations and observations. Even where agreement was weaker (TP, TDP) or dynamics were less pronounced (nitrate, ammonium, Kd), the simulations usually appeared to have the right trend and general range of values, compared to observations. Model predictions of lake behavior The predicted succession among phytoplankton groups is illustrated here for two contrasting locations in the lake. Station 23 is near the deepest point of the oligotrophic east basin (Fig. 1). At a depth of 2 m, which is always within the surface mixed layer, the predicted group succession was from Early Diatoms to a mixture of Late Diatoms and Flagellates and finally to the Other group (Fig. 6a). Predicted biomass of

the Cyanobacteria group was maximal in late summer but in absolute terms was always too low (b0.2 mg Chl a m− 3) for us to show clearly in Fig. 6. At station 970 in the shallower and faster-warming west basin (Fig. 1), the predicted timing of the Early Diatoms peak and subsequent succession was advanced but the sequence was similar to the east basin station (Fig. 6b). The main difference was in the Late Diatoms group, which was relatively more important than at the east basin station. Predicted biomass of Cyanobacteria was again maximal in summer but always low compared to other groups (b0.5 mg Chl a m− 3). The model predicted that temperature was more limiting to the growth rates of all phytoplankton groups than any nutrient at the outset of the simulation. Growth rates for Early Diatoms increased quickly as temperature and light availability increased during spring but then became strongly Si-limited at the time of the predicted peak biomass and subsequent decline. Growth rates of Early Diatoms were severely limited by excessively high temperatures in summer. Growth rates for all other groups were predicted to become strongly P-limited as they approached their seasonal peak biomass, at both the west and east stations. P limitation of growth rates was predicted to persist throughout the simulation period, although it became more moderate in the last 30 days. The predicted dynamics of phytoplankton at the west basin station included very pronounced short-term variations that were not as apparent at the east basin station (Fig. 6) and may have been related to advective motions operating in the more spatially variable west basin (Figs. 2, 3, 7). It is not easy to quantify the role of patch advection in a general, basin-scale, manner as the effects are inevitably specific to the observation point, but we did isolate the advective contributions to daily phytoplankton dynamics for one east basin (station 23) and one west basin (station 970) site for illustrative purposes. A mass

48

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

Fig. 5. Time series output of predicted concentrations of Chl-a, TP, TDP, SRP, SRSi and Kd for the top 5 m together with observations for station 938 (east basin).

balance that included phytoplankton growth, respiration, mortality, settling, resuspension, advection and diffusion was calculated for every day of the simulation season for cells at 2 m depth. Advective processes accounted for almost twice as much of the daily dynamics in the surface mixed layer at the west basin site as at the east, with up to 30% of the daily change in phytoplankton biomass being due to hydrodynamic transport. The predicted seasonal evolution of spatial patterns for Chl a, SRP and TP is shown in Fig. 7 (and cf Figs. 2, 3). In general, the simulations showed highest TP concentrations in west basin and lowest in east basin, but concentrations were far from uniform within basins. TP also tended to decrease from May through June, July and August based on the snapshots shown here. Simulated distributions of Chl a could be quite de-coupled from those of TP. For example, areas of very high TP along the western margins of the west basin in May were quite low in Chl a, despite high concentrations of bioavailable SRP. Moving

eastward along the southern margin of the west basin, predicted TP in the May snapshot remains high but increasing concentrations of Chl a, and decreasing concentrations of SRP, are predicted. In the southwest central basin, high Chl a was predicted in May while TP concentrations were much more moderate and a noticeable depletion of SRP was apparent in the area of high Chl a. In June, highest Chl a was predicted in parts of the northern east basin and not in the west basin where TP was maximal. In these simulations the spatial variation of Chl a within basins can be striking, but the areas of the lake showing the most spatial variation differ among months. In May and August the west and central basins had the strongest spatial heterogeneity, whereas patchiness in June was greatest in the east basin. The figures shown here are only single day snapshots of the spatial distributions but they do exemplify much of the variety in the modeled spatial distributions. Animations are available in App. G to further illustrate the predicted spatial dynamics.

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

Fig. 6. Predicted phytoplankton succession in terms of biomass (treated as chlorophyll a equivalents) at 2 m depth at (a) station 23 in east basin and (b) station 970 in west basin.

A simulated nearshore–offshore study was performed for east basin by computing the average values for a set of nearshore stations and comparing them to the averages for a set of offshore stations (Fig. 1). The comparison showed that nearshore stations were predicted to warm earlier than the offshore (Fig. 8). The simulation period was not long enough to capture the faster cooling that might be expected in the nearshore during the autumn. Nearshore stations were predicted to have higher concentrations of suspended solids

49

(TSS), particularly in the spring. Predicted Chl a increased earlier in the nearshore, was briefly surpassed by the offshore, and was subsequently similar to or higher than offshore concentrations (Fig. 8). The earlier blooming of phytoplankton in the nearshore was reflected in a faster depletion of SRP, but concentrations in nearshore and offshore were similar from July onward. Averages for nearshore and offshore stations by month showed that Chl a was higher nearshore than offshore in all months except June (Table 3). TP was likewise higher nearshore than offshore in most months, but SRP was comparable or lower at the nearshore stations compared to offshore (Table 4). Simulations with east basin tributary inputs set to zero (i.e. a total of eight inflows rather than eleven) gave essentially identical results as exemplified for Chl a in Table 3. The differences in predicted Chl a between zones and months were similar in direction to observations collected in the 1970s, a time when P concentrations and loading were still high and before dreissenid mussels had arrived in the lake (Table 3). The model simulations here correspond to the much lower loadings and spring concentrations of 2002, so to facilitate comparison we have calculated the concentration ratio (as percent) between zones. There were differences in some of the detailed variations but overall the simulations resembled the 1970s observations in showing that nearshore concentrations were usually higher than offshore except in June. The seasonal average ratio of concentrations between zones was remarkably similar between simulated and observed values, with the nearshore averaging about 18% higher in Chl a than the offshore. Historical data to support a similar comparison of TP or SRP could not be found. Discussion Comparison between model and observations Based on snapshot comparisons of spatial patterns and time series observations, ELCD appeared to generate values for all the variables examined that were in the correct range compared to available observations. For some of the key variables (Chl a, TP, SRP and SRSi)

Fig. 7. Model predictions for May 15 (a, d, g), June 15 (b, e, h) and July 15 (c, f, i), showing vertically-averaged values of total P (a, b, c), Chl a (d, e, f) and SRP(g, h, i), all as mg m− 3.

50

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

Fig. 8. Time series of nearshore and offshore epilimnetic (0–5 m) temperatures and concentrations of TSS, Chl a and SRP in the eastern basin of Lake Erie.

the agreement with observations was very close and model results reproduced most of the observed seasonal dynamics. This agreement is partly the result of some informal calibration in which we adjusted certain model parameters, mainly the rate constants for phytoplankton mortality and remineralization of detritus, to improve agreement with time series observations at offshore stations in the east basin. These parameters are difficult to define with reference to the literature because they are simplified summations of complicated ecological processes (e.g. grazing, parasitism, microheterotroph metabolism, etc.) that are not explicit in the model. No tuning was done with respect to observations in the west or central basins, or at nearshore stations in the east basin, so the agreement observed in those areas provides some evidence that the model is capable of making useful predictions outside its range of calibration. Comparisons against other years of observation are planned, and recent implementation of the model (with the same parameter values) in Lake Ontario has shown good agreement between simulations and observations (Leon et al., 2009 and manuscript in prep.). Like Schwab et al. (2009), however, we do not present the model as mechanistically faithful to the complexities of the lake ecosystem. Rather it is a means to portray the implications of the interactions between physical processes (hydrodynamics, temperature and irradiance

variations) and a relatively simple, but defensible, set of assumptions concerning the ecological processes of interest. The major systematic bias that we detected in the ELCD model simulations was an overestimation of TP during the summer stratified season in the east and, especially, central basin. We did not observe a similar bias in west basin, where TP concentrations never fall as low and persistent stratification does not develop. The model of Schwab et al. (2009) also overestimated TP at central and eastern stations with low TP. This may reflect incorrect assumptions about the nature of the TP pool in both models. In Schwab et al. (2009) TP was treated as a single physical entity with a uniform settling velocity. In our more mechanistically-resolved model, we assumed (for lack of observational data) a small (1 mgP m− 3) concentration of sediment-bound P in the lake at the beginning of our simulation, but none was specified in the inflows. The vast majority of P in our simulations was in bioavailable and/or organic forms (dissolved phosphate, labile dissolved and particulate organic P, and phytoplankton internal P). Suspended particles in the Niagara River downstream of Lake Erie were, by contrast, less than 50% bioavailable while mineral (apatite) P comprised over 25% of the total particulate P (Mayer et al., 1991). Stone and English (1993) found apatite to comprise an even larger fraction of particulate P in two small Lake Erie tributaries, especially in

Table 3 Nearshore–offshore comparison of chlorophyll-a (Chl-a, mg m− 3) in the eastern basin of Lake Erie and the difference between nearshore and offshore Chl-a concentrations expressed as the percent of an average nearshore–offshore concentration (%). Month

Apr May Jun Jul Aug Sep Oct Average

Observed1 (1973, 1974 and 1975)

Predicted* (2002) (11 inflows)

Predicted** (2002) (8 inflows)

Nearshore (b 25 m)

Offshore (N 25 m)

%

Nearshore (b 20 m)

Offshore (N20 m)

%

Nearshore (b 20 m)

Offshore (N 20 m)

3.8 7.3 3.6 4.1 4.9 5.1 7.0 5.1

2.7 6.9 3.8 2.7 4.1 4.1 6.3 4.4

33.0 5.9 − 6.1 42.2 16.7 21.4 10.4 17.64

0.8 2.7 3.5 2.2 2.3 2.2 2.1 2.3

0.4 1.7 4.1 2.1 2.3 2.0 1.7 2.0

60.1 46.4 − 15.2 6.7 0.1 9.7 24.2 18.90

0.8 2.8 3.5 2.3 2.2 2.1 2.1 2.2

0.4 1.7 4.0 2.1 2.2 1.9 1.6 2.0

a Observed values are reported as monthly averages for the years of 1973, 1974 and 1975 (EPA, 1980); * results with 11 inflows; ** results with 8 inflows, excluding all direct inflows to east basin.

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

larger particles, while the more available inorganic fraction comprised less than 27% of the sediment-bound P. It seems certain that much of the tributary load of P is not highly bioavailable and that a considerable fraction is in mineral forms that can be expected to settle out during the seasonal stratified period. With inputs and initial concentrations as we specified them, the model would not fully capture such processes. Published data on P fractions in the larger tributaries and the lake itself that could inform better setup of inputs and initial conditions unfortunately seem scarce. Pending collection of such measurements, future work could use the ELCD model to estimate the sensitivity of P dynamics in the lake to different assumptions, within the plausible range, about the functional composition of the TP. Model predictions of phytoplankton succession are challenging to compare against available observations, partly because most of the model groups are somewhat abstract relative to conventional taxonomic groupings, but predictions appeared reasonable where meaningful comparisons were possible. The predicted dominance of Early Diatoms in the spring phytoplankton maximum was consistent with expectations (e.g. Barbiero and Tuchman, 2001; Carrick, 2004; Barbiero et al., 2006; Smith et al., 2007). The timing and magnitude of the bloom, as judged by comparison against Chl a measurements, seemed reasonable (e.g. Fig. 5) and analysis of samples in east basin confirmed that diatoms were dominant during the spring biomass maximum (Smith et al., 2007). The Cyanobacteria were predicted to remain a small (b11%) part of the phytoplankton biomass at the two stations examined, despite the known occurrence of cyanobacterial blooms in west and western central basins (e.g. Budd et al., 2001). Blooms are episodic and spatially limited events, however, and available data support the model predictions of a normally low cyanobacterial contribution to biomass even at the west basin station (Makarewicz et al., 1999; Barbiero and Tuchman, 2001; Ghadouani and Smith, 2005). Our simulated west basin station was not situated in the areas that seem to be most prone to late summer Microcystis blooms and results may have differed if it were placed differently. Model predictions of controls on phytoplankton growth rates were evaluated only at a few reference points in the lake, and always at a depth of 2 m where light availability is expected to be good. Evaluation of growth limiting factors over the whole surface mixed layer would be desirable but is beyond the scope of the current study. It will be pursued in forthcoming work that will use the model to predict lake responses to altered nutrient loading. Even the more limited examination made here has provided insights that can be compared against relevant observations. The model predicted that phytoplankton growth should be limited primarily by temperature and/or light in the early part of the simulation period (April–May), a time when temperature and mean irradiance in the surface layer are typically low (e.g. Guildford et al., 2005). The major decline of the Early Diatom spring bloom was predicted to be accompanied by strong silicate limitation of growth rates in late spring and there is indeed evidence that SRSi concentrations at that time can be strongly depleted (e.g. Guildford et al., 2005 and Fig. 5). Temperatures subsequently increased into the supraoptimal range for Early Diatoms, preventing their re-growth until autumn when temperatures decreased sufficiently that P became the factor most limiting to their growth. All other groups were predicted to become strongly Plimited in their growth rates as they reached and/or passed their seasonal maximum biomass. The prevalence of P limitation was expected because historical responses of the lake to variations of P loading rates support the view that it is a P limited lake (Charlton et al., 1999; Matisoff and Ciborowski, 2005). Experimental assessments of phytoplankton physiology through bioassays, enzyme assays, and stoichiometric measurements in recent (post-1995) times provide more direct evidence for a prevalence of P limitation in summer time (e.g. Guildford et al., 2005; Moon and Carrick, 2007; Rattan, 2009). The model predicted substantial easing of P limitation in the autumn (October), also consistent with direct assays of phytoplankton

51

physiology (Guildford et al., 2005; Rattan, 2009). Detailed comparisons among locations within the lake will form a part of the forthcoming work on model predictions of nutrient loading responses, but the general sequence of temperature, silicate and P limitation described here pertained to both the east and west basin reference points we examined. Scales and patterns of variability Averaged over the typical observation period (April to October) Lake Erie has a west to east gradient of decreasing nutrient and phytoplankton concentrations (e.g. Charlton et al., 1999; Barbiero and Tuchman, 2001), reflecting the fact that most of the nutrient load enters at the western end of the lake. Similar to a recently published three dimensional model for phosphorus (Schwab et al., 2009), our model captured this general longitudinal gradient. Our simulations further showed that TP concentrations were highly dynamic across space and time; within-basin variation often rivaled or exceeded between-basin variation (e.g. Fig. 7). Strong spatial variability has been long been apparent in features observable by remote sensing (e.g. Budd et al., 2001) and models such as ELCD and that of Schwab et al. (2009) seem capable of capturing much of this variability. The west and west-central basins displayed the most spatial variability, consistent with the typical experience from lake-wide, ship-borne, sampling (e.g. Smith et al., 2005). There were also occasions of very strong spatial variability in the central and east basin. For example, strong variability of Chl a concentrations was predicted to occur along the northern coastline in June (Fig. 7d). Such variability could easily be missed by routine offshore sampling, but was captured by our intensive sampling efforts during 2002, which confirmed that June was indeed a time of exceptionally high spatial variability in the east basin (Depew et al., 2006). As noted by Schwab et al. (2009), the ability to characterize this kind of variability is an important step towards better understanding of meso-scale phenomena such as algal blooms. A prominent source of variability in the west basin, as expected, was the influence of the Detroit and Maumee Rivers (Fig. 7). With nutrient concentrations being far higher in the Maumee than the Detroit, it might be expected that the southern and western portions of the basin would have chronically elevated TP and Chl a compared to the northern and eastern parts. Similar to results of Schwab et al. (2009) for TP, the current results portrayed a more complicated situation. The highly variable discharge of the Maumee, which can fall to extremely low values in summer, together with the variable winddriven circulation can produce an almost reverse pattern of distribution of TP at times. The predicted distributions of Chl a were also variable and were often quite different from the TP distributions in the west basin. This model result is at least partly a consequence of the stipulated river loadings, which bring relatively more P (including bioavailable SRP) than Chl a into the lake. Chl a can then develop, depending on the conditions affecting algal dynamics, but with some inevitable lag. Some tributaries (e.g. Grand River, Ontario, R. Smith unpub. data) can at times have very high Chl a concentrations, but we did not have sufficient data to attempt to simulate the effects of riverine inputs of high Chl a water. River Chl a concentrations tend to be highest in the season of lowest discharge, which may limit their direct impact on lake concentrations, but they are certainly one factor (along with improved knowledge of P fractions and bioavailability in tributary inputs) that needs further study if we are to properly understand the variable coupling of TP and Chl a in this area of the lake. Knowledge of the predicted spatial dynamics of phytoplankton growth limitation, and its controlling factors, would also inform our understanding but to extract and synthesize this information is a substantial task that must be deferred to our forthcoming model study of nutrient controls on phytoplankton in the lake. For this spatially complex area of the lake, a more authentic reproduction of events may also require a finer resolution than afforded with the current

52

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

Table 4 Average simulated concentrations of TP and SRP in the top 5 m of nearshore and offshore zones of east basin of Lake Erie. Month

Apr May Jun Jul Aug Sep Oct

TP (mg m− 3)

SRP (mg m− 3)

Offshore

Nearshore

Offshore

Nearshore

13.1 12.9 10.5 8.1 7.1 6.7 7.3

13.5 13.3 10.4 8.3 7.3 6.9 7.1

8.5 7.9 3.5 1.1 0.1 0.1 0.3

8.5 7.3 3.0 1.0 0.1 0.1 0.1

2 km grid and development of a higher resolution (500 m grid) model for the western part of the lake is currently underway. Despite the evident scope for improving the estimates of boundary conditions (including tributary inputs) and spatial resolution of the model, the prediction of a loose coupling between the variations of P and Chl a within basins, and especially west basin, does find support in measurements from the lake. Using the same observational data used here for model validation, we tested for correlations among Chl a and TP measurements using all available stations and dates. The correlations were not significant for any of the basins. For east and central basins, a highly significant (p b 0.001) positive correlation was observed if April values (corresponding to conditions preceding the spring bloom) were excluded but the explanatory power was low (Pearson r ≤ 0.7). Using model predicted TP and Chl a for the same stations and dates, we obtained very similar results with highly significant correlations but low explanatory power (Pearson r ≤ 0.8) for each of the east and central basins. Neither observational nor model predicted values yielded a significant correlation between Chl a and TP in west basin for this more restricted period of observation. Consistent with these results, spatial surveys conducted in 1997 mainly during summer and early fall revealed a significant correlation between Chl a and particulate P in central and east basins but not in west basin (Guildford et al., 2005). While there is good evidence that TP and Chl a concentrations in Lake Erie, including west basin, are coupled over periods of multiple years (e.g. Charlton et al., 1999), the looser relationship within years and basins portrayed by the model does seem to be realistic. The current model resolution may be more appropriate in the east basin, with its simpler bathymetry, larger size and lesser degree of direct tributary influence. The model predictions of systematically higher chlorophyll and TP in nearshore vs. offshore areas of the east basin was consistent with other model results (Schwab et al., 2009) that show a tendency for tributary nutrient inputs to be more or less coast-bound in all three basins. Of the eleven inflows we included only three were in the east basin (Fig. 1). The Grand River (ON) is the largest of them and does, on occasion, produce a recognizable plume of higher TP in our results (e.g. Fig. 3) but it is a relatively small load in the lake budget, usually highly localized and never reaches the south coast. Over 95% of the P load in our model enters the lake in the west or west-central basin and can only reach the east basin by longdistance transport. Model runs with no east basin tributary inputs at all gave virtually the same pattern of nearshore–offshore concentration differences in the east basin. The observation from the model that TP can be systematically elevated in the nearshore of the east basin even with no simulated local inputs suggests that other factors, such as the balance between sedimentation and resuspension, are important influences in this region. The difference in Chl a was even greater, suggesting that the shallower nearshore zone is more conducive to phytoplankton production than the offshore, assuming that phytoplankton mortality does not differ systematically. The elevated concentrations of TP and Chl a in the nearshore predicted by the model are consistent with patterns reported for

many parts of the Great Lakes prior to arrival of dreissenid mussels (Depew et al., 2006) and specifically (in the case of Chl a) for east basin Lake Erie (US EPA, 1980). At least in east basin, the pattern after arrival of the mussels has been different, with Chl a being lower in the nearshore (Depew et al., 2006). While the change in this pattern could be due to mussels (Hecky et al., 2004) it is also possible that P loading controls of various kinds acted to diminish the degree of nearshore enrichment exerted by diffuse nutrient sources, as postulated for Lake Ontario (Hall et al., 2003). Studies of loading changes to Lake Erie have shown, however, that loading controls were most effective on the larger point and tributary sources in the western part of the lake, and that diffuse loading subsequently comprised a larger, not smaller, part of the total load (Dolan and McGunagle, 2005). This suggests that changes in loading should have, if anything, accentuated the historic pattern by preferentially increasing the relative importance of local, diffuse, inputs (i.e. higher nutrient and Chl a in the nearshore than the offshore would be expected). The model results give another line of evidence that the change to lower Chl a nearshore relative to offshore cannot be explained by a decrease of diffuse loading. Even with no direct loading to the east basin, the model still predicted higher Chl a (and TP) in the nearshore. While indirect, this evidence supports the view that mussels may indeed have been critical to the observed reversal of the nearshore–offshore zonation and is the subject of further research and model development. This interpretation is consistent with a recent meta-analysis of dreissenid mussel impacts on lake and river ecosystems, which indicated that phytoplankton biomass was suppressed in littoral habitats by approximately 50% of pre-invasion values (Higgins and Vander Zanden 2010). These results infer that dreissenid mussel impacts in large lake ecosystems may be dramatic within nearshore zones, even if offshore zones show little effect. Thus, offshore monitoring programs and whole lake well-mixed reactor models are unlikely to capture the dramatic changes in nearshore water quality. Successful modeling and management of Great Lake coastal waters will require a greater understanding of the interactions between invasive dreissenids and physical (water circulation, resuspension), chemical (nutrients), and biological (phytoplankton biomass and community structure) parameters. The ELCD model currently uses first order rate constants, modified only by temperature and dissolved oxygen, to describe phytoplankton mortality and nutrient remineralization. If agents such as dreissenid mussels are truly responsible for large scale impacts such as the change in nearshore–offshore patterns, then an important and spatially-variable agent is lacking in the model. Both models (Boegman et al., 2008) and observation (e.g. Nicholls, 1999) show that zooplankton can change the coupling between phytoplankton and phosphorus and, while the current model does simulate transfer of phytoplankton nutrients into detrital pools, it will not capture the spatial dynamics of zooplankton effects accurately. The current model tended to overestimate Chl a more at nearshore than offshore stations in the east basin in summer but not spring, which could reflect both spatial and temporal patterns in the strength of mussel effects on phytoplankton. Work is currently underway to add a portrayal of both dreissenid mussels and zooplankton to ELCD. Recent work with two dimensional models in Lake Erie (Boegman et al., 2008; Zhang et al., 2008) shows that mussel and zooplankton activity can have major impacts on how models portray spatial and temporal variations, and we expect that a three dimensional analysis will provide further insights into the role of the mussels in lake ecosystem dynamics. Supplementary materials related to this article can be found online at doi:10.1016/j.jglr.2010.12.007. Acknowledgements This work was supported primarily by grants to R.E.H.S. from the Natural Sciences and Engineering Research Council of Canada (Strategic

L.F. Leon et al. / Journal of Great Lakes Research 37 (2011) 41–53

Projects and Special Opportunities Programs). It was also supported by a grant to R.E.H. from the US E.P.A. (Lake Erie Trophic Status Study, or LETSS) and in-kind support (e.g. meteorological data buoy) to R.E.H.S. from Environment Canada (National Water Research Institute). D. Rockwell and M. Tuchman of US E.P.A. as well as R. Barbiero (DynCorp Inc., Chicago) kindly assisted in making agency survey data available, and H. Carrick (Penn State University) provided his personal collation of water chemistry measurements stemming from LETSS. M. Charlton and D. Lam (N.W.R.I.) facilitated access to Environment Canada survey data. Our ability to assess the model and derive insights from it would have been far more limited without the generous assistance of all the foregoing individuals, as well as the U. Waterloo colleagues, students and technicians who worked on the U.W. Lake Erie project. References Arhonditis, G.B., Brett, M.T., 2005. Eutrophication model for Lake Washington (USA) Part I. Model description and sensitivity analysis. Ecol. Mod. 187, 140–178. Barbiero, R.P., Tuchman, M.L., 2001. Results from the U.S. EPA's biological open water surveillance program of the Laurentian Great Lakes: I. Introduction and phytoplankton results. J. Great Lakes Res. 27, 134–154. Barbiero, R.P., Rockwell, D.C., Warren, G.J., Tuchman, M.L., 2006. Changes in spring phytoplankton communities and nutrient dynamics in the eastern basin of Lake Erie since the invasions of Dreissena spp. Can. J. Fish. Aquat. Sci. 63, 1549–1563. Beletsky, D., Schwab, D.J., 2001. Modeling circulation and thermal structure in Lake Michigan: annual cycle and interannual variability. J. Geophys. Res. 106 (C9), 19745–19771. Beletsky, D., Mason, D.M., Schwab, D.J., Rutherford, E.S., Janssen, J., Clapp, D.F., Dettmers, J.M., 2007. Biophysical model of larval yellow perch advection and settlement in Lake Michigan. J. Great Lakes Res. 33, 842–866. Bierman Jr., V.J., Dolan, D.M., 1981. Modeling of phytoplankton-nutrient dynamics in Saginaw Bay, Lake Huron. J. Great Lakes Res. 7, 409–439. Boegman, L., Loewen, M.R., Culver, D.A., Hamblin, P.F., Charlton, M.N., 2008. Spatialdynamic modeling of algal biomass in Lake Erie: relative impacts of dreissenid mussels and nutrient loads. J. Environ. Eng. 134, 456–468. Bolsenga, S.J., Herdendorf, C.E., 1993. Lake Erie and Lake Saint Clair handbook. Wayne State Univ. Press. Budd, J.W., Beeton, A.W., Stumpf, R.P., Culver, D.A., Kerfoot, W.C., 2001. Satellite observations of Microcycstis blooms in western Lake Erie. Verh. Internat. Verein. Limnol. 27, 3787–3793. Carrick, H.J., 2004. Algal distribution patterns in Lake Erie: implications for oxygen balances in the eastern basin. J. Great Lakes Res. 30, 133–147. Charlton, M.N., Le Sage, M.N., Milne, J.E., 1999. Lake Erie in transition: the 1990's. In: Munawar, M., Edsall, T., Munawar, I.F. (Eds.), State of Lake Erie, Past Present and Future. Backhuys Publishers, Leiden, Netherlands, pp. 74–114. Chen, C., Ji, R., Schwab, D.J., Beletsky, D., Fahnenstiel, G.L., Jiang, M., Johengen, T.H., Vanderploeg, H., Eadie, B., Budd, J.W., Bundy, M.H., Garnder, W., Cotner, J., Lavrentyev, P.J., 2002. A model study of the coupled biological and physical dynamics in Lake Michigan. Ecol. Mod. 152, 145–168. Conroy, J.D., Kane, D.D., Dolan, D.M., Edward, W.J., Charlton, M.N., Culver, D.A., 2005. Temporal trends in Lake Erie plankton biomass: roles of external phosphorus loading and dreissenid mussels. J. Great Lakes Res. 31 (suppl. 2), 89–110. Depew, D.C., Smith, R.E.H., Guildford, S., 2006. Nearshore–offshore comparison of chlorophyll a and phytoplankton production in the dreissenid-colonized eastern basin of Lake Erie. Can. J. Fish. Aquat. Sci. 63, 1115–1129. Dolan, D.M., McGunagle, K.P., 2005. Lake Erie total phosphorus loading analysis and update: 1996–2002. J. Great Lakes Res. 31 (supp 2), 11–22. Fahnenstiel, G.L., Chandler, J.F., Carrick, H.J., Scavia, D.J., 1989. Photosynthetic characteristics of phytoplankton communities in lakes Huron and Michigan: P–I parameters and end-products. J. Great Lakes Res. 15, 394–407. Gal, G., Hipsey, M.R., Paparov, A., Makler, V., Zohary, T., 2009. Implementation of ecological modeling as an effective management and investigation tool: Lake Kinneret as a case study. Ecol. Mod. 220, 1697–1718. Geider, R.J., Osborne, B.A., 1989. Respiration and microalgal growth: a review of the quantitative relationship between dark respiration and growth. New Phytol. 112, 327–341. Ghadouani, A., Smith, R.E.H., 2005. Phytoplankton distribution in Lake Erie as assessed by a new in situ spectrofluorometric technique. J. Great Lakes Res. 31 (suppl. 2), 154–167. Guildford, S.J., Hecky, R.E., Smith, R.E.H., Taylor, W.D., Charlton, M.N., Barlow-Busch, L., North, R.L., 2005. Phytoplankton nutrient status in Lake Erie in 1997. J. Great Lakes Res. 31 (supp. 2), 72–88. Hall, S.R., Pauliukonis, N.K., Mills, E.L., Rudstam, L.G., Schneider, C.P., Lary, S.J., Arrhenius, F., 2003. A comparison of total phosphorus, chlorophyll a and zooplankton in embayments, nearshore, and offshore habitats of Lake Ontario. J. Great Lakes Res. 29, 54–69. Hecky, R.E., Smith, R.E.H., Barton, D.R., Guildford, S.A., Taylor, W.D., Howell, T.D., Charlton, M.N., 2004. The nearshore shunt: a consequence of ecosystem

53

engineering by dreissenids in the Laurentian Great Lakes. Can. J. Fish. Aquat. Sci. 61, 1285–1293. Higgins, S.N., Malkin, S.Y., Howell, E.T., Guildford, S.J., Campbell, L., Hiriart-Baer, V., Hecky, R.E., 2008. An ecological review of Cladophora glomerata (Chlorophyta) in the Laurentian Great Lakes. J. Phycol. 44, 839–854. Higgins, S.N., VanderZanden, M.J., 2010. What a difference a species makes: a metaanalysis of dreissenid mussel impacts on freshwater ecosystems. Ecol. Monogr. 80, 179–196. Hipsey, M.R., Romero, J.R.R., Antenucci, J.P., Imberger, J., 2006. The computational aquatic ecosystem dynamics model (CAEDYM): a versatile water quality model for coupling with hydrodynamic drivers. Proc. 7th International Conference on Hydroinformatics, Nice, France. Hodges, B.R., Imberger, J., Saggio, A., Winters, K., 2000. Modeling basin scale waves in a stratified lake. Limnol. Oceanogr. 45, 1603–1620. Istvanovics, V., Padisak, J., Pettersson, K., Pierson, D.C., 1994. Growth and phosphorus uptake of summer phytoplankton in Lake Erken (Sweden). J. Plank. Res. 16, 1167–1196. Lehman, J.D., Botkin, D.B., Likens, G.E., 1975. The assumptions and rationale of a computer model of phytoplankton population dynamics. Limnol. Oceanogr. 20, 343–364. Leon, L.F., Imberger, J., Smith, R.E.H., Lam, D., Schertzer, W., 2005. Modeling as a tool for nutrient management in Lake Erie. J. Great Lakes Res. 31 (5), 309–318. Leon, L., Smith, R.E.H., Romero, J., 2006. Lake Erie hypoxia simulations with ELCOM– CAEDYM. Proceedings of 3rd Biennial meeting of Int. Environ. Model. Software Soc., July 2006, Burlington, VT, USA, pp. 1–6. Leon, L.F., Smith, R.E.H., Malkin, S., Depew, D., Hecky, R.E., 2009. Modeling and analysis of Cladophora dynamics and their relationship to local nutrient sources in a nearshore segment of Lake Ontario. Final Report to Ontario Power Generation, Pickering, Ont. Makarewicz, J.C., Lewis, T.W., Bertram, P., 1999. Phytoplankton composition and biomass in the offshore waters of Lake Erie: pre- and post-Dreissena introduction (1983–1993). J. Great Lakes Res. 25, 135–148. Matisoff, G., Ciborowski, J.J.H., 2005. Lake Erie trophic status collaborative study. J. Great Lakes Res. 31 (suppl. 2), 1–10. Matisoff, G., Neeson, T.M., 2005. Oxygen concentration and demand in Lake Erie sediments. J. Great Lakes Res. 31 (suppl. 2), 284–295. Mayer, T., Kuntz, K.W., Moller, A., 1991. Total and bioavailable particulate phosphorus loads from the Niagara River in 1987 and 1988. J. Great Lakes Res. 17, 446–453. Mieleitner, J., Reichert, P., 2008. Modelling functional groups of phytoplankton in three lakes of different trophic state. Ecol. Mod. 211, 279–291. Millie, D.F., Fahnenstiel, G.L., Carrick, H.J., Lohrenz, S.E., Schofield, O.M.E., 2002. Phytoplankton pigments in coastal Lake Michigan: distributions during the spring isothermal period and relation with episodic sediment resuspension. J. Phycol. 38, 639–648. Moon, J.N., Carrick, H.J., 2007. Seasonal variation of phytoplankton nutrient limitation in Lake Erie. Aquat. Microb. Ecol. 48, 61–71. Morillo, S., Imberger, J., Antenucci, J.P., Woods, P.F., 2008. Influence of wind and lake morphometry on the interaction between two rivers entering a stratified lake. J. Hydraul. Eng. 134, 1579–1589. Nicholls, K.H., 1999. Evidence for a trophic cascade effect on north shore western Lake Erie phytoplankton prior to the zebra mussel invasion. J. Great Lakes Res. 25, 942–949. Rao, Y.R., Schwab, D.J., 2007. Transport and mixing between the coastal and offshore waters in the Great Lakes: a review. J. Great Lakes Res. 33, 202–218. Rao, Y.R., Hawley, N., Charlton, M.N., 2008. Physical processes and hypoxia in the central basin of Lake Erie. Limnol. Oceanogr. 53, 2007–2020. Rattan, K. J. Traditional and new fluorometric methods to determine phytoplankton nutrient status for freshwater ecosystems. PhD Thesis, Biology Department, University of Waterloo, Waterloo, Ont. Canada. Reynolds, C.S., 1984. The Ecology of Freshwater Phytoplankton. Cambridge U. Press, Cambridge. Robson, B.J., Hamilton, D.P., 2004. Three-dimensional modelling of a Microcystis bloom event in the Swan River estuary, Western Australia. Ecol. Mod. 174, 203–222. Romero, J.R., Antenucci, J.P., Imberger, J., 2004. One- and three-dimensional biogeochemical simulations of two differing reservoirs. Ecol. Mod. 174, 143–160. Schwab, D.J., Beletsky, D., DePinto, J., Dolan, D.M., 2009. A hydrodynamic approach to modeling phosphorus distribution in Lake Erie. J. Great Lakes Res. 35, 50–60. Smith, R.E.H., Kalff, J., 1983. Competition for phosphorus among co-occurring freshwater phytoplankton. Limnol. Oceanogr. 28, 448–464. Smith, R.E.H., Hiriart-Baer, V., Higgins, S., Guildford, S., Charlton, M., 2005. Planktonic primary production in the offshore waters of dreissenid-infested Lake Erie, 1997. J. Great Lakes Res. 31 (supp. 2), 50–62. Smith, R.E.H., Parrish, C.C., Depew, D.C., Ghadouani, A., 2007. Spatial patterns of seston concentration and biochemical composition between nearshore and offshore waters of a Great Lakes. Freshwat. Biol. 52, 2196–2210. Spillman, C.M., Hamilton, D.P., Hipsey, M.R., Imberger, J., 2008. A spatially resolved model of seasonal variations in phytoplankton and clam (Tapes philippinarum) biomass in Barbamarco Lagoon, Italy. Est. Coast. Shelf Sci. 79, 187–203. Stone, M., English, M.C., 1993. Geochemical composition, phosphorus speciation and mass transport of fine-grained sediment in two Lake Erie tributaries. Hydrobiol. 253, 17–29. U.S. Environmental Protection Agency (US EPA), 1980. Lake Erie nutrient control effectiveness regarding assessment in east basin. Performing Organization S.U.N.Y. Buffalo, N.Y. EPA-600/3-80-067. Zhang, H., Culver, D.A., Boegman, L., 2008. A two-dimensional ecological model of Lake Erie: application to estimate dreissenid impacts on large lake plankton populations. Ecol. Mod. 214, 219–241.