Multi-Year Simulation of Western Lake Erie Hydrodynamics and Biogeochemistry to Evaluate Nutrient Management Scenarios

: During the 1970s, harmful cyanobacteria (HFCB) were common occurrences in western Lake Erie. Remediation strategies reduced total P loads and bloom frequency; however, HFCB have reoccurred since the mid-1990s under increased system stress from climate change. Given these concurrent changes in nutrient loading and climate forcing, there is a need to develop management tools to investigate historical changes in the lake and predict future water quality. Herein, we applied coupled one-dimensional hydrodynamic and biogeochemical models (GLM–AED) to reproduce water quality conditions of western Lake Erie from 1979 through 2015, thereby removing the obstacle of setting and scaling initial conditions in management scenarios. The physical forcing was derived from surface buoys, airports, and land-based stations. Nutrient loads were reconstructed from historical monitoring data. The root-mean-square errors between simulations and observations for water levels (0.36 m), surface water temperature (2.5 ◦ C), and concentrations of total P (0.01 mg L − 1 ), PO 4 (0.01 mg L − 1 ), NH 4 (0.03 mg L − 1 ), NO 3 (0.68 mg L − 1 ), total chlorophyll a (18.74 µ g L − 1 ), chlorophytes (3.94 µ g L − 1 ), cyanobacteria (12.44 µ g L − 1 ), diatoms (3.17 µ g L − 1 ), and cryptophytes (3.18 µ g L − 1 ) were minimized using model-independent parameter estimation, and were within literature ranges from single year three-dimensional simulations. A sensitivity analysis shows that 40% reductions of total P and dissolved reactive P loads would have been necessary to bring blooms under the mild threshold (9600 MTA cyanobacteria biomass) during recent years (2005–2015), consistent with the Annex 4 recommendation. However, these would not likely be achieved by applying best management practices in the Maumee River watershed.


Introduction
Lake Erie, the shallowest and most productive Laurentian Great Lake, has been suffering from eutrophication for the past half century, particularly in the shallow western basin [1]. Eutrophication is primarily driven by excess nutrient loads, particularly phosphorus (P), from agriculture, domestic wastewater, and industrialization [2]. Observations of P limitation in the lower Great Lakes (including western Lake Erie) have resulted in P control policies designed to reverse eutrophication [3].
From the early 1980s to the early 1990s, point-source P abatement programs were implemented as part of the Great Lakes Water Quality Agreement (GLWQA) of 1972 [4]. The 1978 Amendment to the GLWQA set a Lake Erie target total phosphorus (TP) load of 11,000 MTA. As a response, the external P loading entering into the western basin declined, leading to reduced total P concentrations in the water column [1], a decline in phytoplankton biomass [5,6], and improved water quality in the western basin. However, beginning in the mid-1990s, the algal blooms returned [7,8], their recurrence linked to increased spring precipitation flushing increased soluble P from the Maumee River watershed into Lake Erie [9,10]. These changes in nutrient loads are consistent with long-term and predicted tion in the development of nuisance algae blooms in Lake Erie since 1979; (3) evaluate the long-term effectiveness of nutrient load reductions and BMPs on reducing bloom severity. This work is novel, in that it is the longest coupled hydrodynamic and biogeochemical simulation of Lake Erie and the first to employ an automated calibration approach. A full list of acronyms is given in Table 1.

Study Site and Field Data
Lake Erie is the southernmost, shallowest, and smallest by volume of the Great Lakes, and has distinct western, central, and eastern basins. The present study area is the western basin ( Figure 1), which has a surface area of 4837 km 2 , a volume of~20 km 3 , and an average depth of 7.4 m. It has two major tributaries, the Detroit River ( Figure 1 point A), accounting for approximately 90% of the total annual inflow [25], and the Maumee River ( Figure 1 point B), accounting for approximately 47% of the TP (total phosphorus) loads into Lake Erie during 2011-2013 [12]. The TP concentration is 25 times larger in the Maumee River compared to the Detroit River [12]; hence, the Maumee River is the main P source for western Lake Erie. The theoretical hydraulic residence time is~51 days [26], with outflow to the central basin through a rocky chain of islands from Point Pelee, Ontario, to Marblehead, Ohio.

Study Site and Field Data
Lake Erie is the southernmost, shallowest, and smallest by volume of the Great Lakes, and has distinct western, central, and eastern basins. The present study area is the western basin ( Figure 1), which has a surface area of 4837 km 2 , a volume of ~20 km 3 , and an average depth of 7.4 m. It has two major tributaries, the Detroit River ( Figure 1 point A), accounting for approximately 90% of the total annual inflow [25], and the Maumee River ( Figure  1 point B), accounting for approximately 47% of the TP (total phosphorus) loads into Lake Erie during 2011-2013 [12]. The TP concentration is 25 times larger in the Maumee River compared to the Detroit River [12]; hence, the Maumee River is the main P source for western Lake Erie. The theoretical hydraulic residence time is ~51 days [26], with outflow to the central basin through a rocky chain of islands from Point Pelee, Ontario, to Marblehead, Ohio.

Character
Location Source Identifier in Figure 1 Sample Date Observational data for model forcing and calibration ( Figure 1) have been compiled from the Great Lakes Environmental Research Laboratory (GLERL), the National Oceanic and Atmospheric Administration (NOAA), Environment and Climate Change Canada (ECCC), the Environmental Monitoring and Reporting Branch (EMRB), the Ontario Clean Water Agency (OCWA), and published scientific literature (Tables S1 and S2).

Model Description
The 1D hydrodynamic model General Lake Model (GLM v. 3.0.5), coupled with the Aquatic Ecosystem Dynamics module library (AED), was applied in this study [23].
These open-source models have an active development community through the Global Lakes Ecological Observatory Network (GLEON; [23]), and are updated versions of the well-tested and commonly used DYRESM-CAEDYM [31,32]. GLM-AED has been applied to a wide range of water bodies, including simulation of algae/nutrient/oxygen dynamics in temperate lakes [33] and reservoirs [34]. The low computational requirements for a single run allow for extensive calibration and sensitivity analysis.
GLM is a mixed layer model that solves a turbulent kinetic energy balance to simulate the surface heat and momentum budgets and the resultant development of thermal stratification and mixing, including ice and snow cover. The model assumes there is no horizontal variability (horizontally averaged) and adopts a flexible Lagrangian structure, adjusting the vertical layer thicknesses dynamically to resolve the water column structure. AED is a biogeochemical nutrient/plankton/oxygen model that dynamically couples to GLM.
In the present AED setup, 11 state variables were employed to predict seasonal succession of phytoplankton biomass (greens (GREEN), diatoms (DIAT), cryptophytes (CRYPT), and cyanobacteria (CYANO)). These included dissolved oxygen (DO), four dissolved inorganic groups (dissolved reactive phosphorus (PO 4 ), nitrate (NO 3 ), ammonium (NH 4 ), and reactive silica (RSi)), three dissolved organic groups (dissolved organic nitrogen (DON), dissolved organic phosphorus (DOP), and dissolved organic carbon (DOC)), and three particulate detrital organic groups (particulate organic nitrogen (PON), particulate organic phosphorus (POP), and particulate organic carbon (POC)). The phytoplankton growth rates were modelled as a function of the user specified maximum growth rate, photorespiratory loss, ambient temperature, metabolic stress, and the minimum light, N, P, and Si limitation functions [35]. The associated dynamic functions representing algae growth rates (photosynthesis), respiration, excretion, and mortality losses can be found in Hipsey et al. [35].
Photosynthesis is parameterized as the uptake of carbon, and it is determined by a maximum potential growth rate at 20 • C modified by photorespiratory loss, a temperature response function, metabolic stress, and the minimum value of the expressions for limitation by light, P, N, and Si. There is maximum productivity at the optimum growth temperature (T OPT ), zero growth above the maximum growth temperature (T MAX ), and below standard growth temperature (T STD ); the productivity follows Arrhenius scaling. P and N uptake are regulated by external and internal nutrient concentrations. Loss terms, including respiration, natural mortality, and excretion, are modelled with respiration rate coefficients, and the loss rate is divided into two parts: a pure respiratory fraction f res , as well as mortality and excretion. The constant fraction f DOM of mortality and excretion goes to the dissolved organic pool (excretion), and the remainder to the particulate organic pool (detritus). Settling of particles (−0.06 m day −1 ) and mineralization of dissolved organics were modeled as migration with photoinhibition (Rdoc_minerl = 0.001 day −1 ; Rdon_minerl = 0.005 day −1 ; Rdop_minerl = 0.001 day −1 ). Given the difficulty in modelling changes in the associated population dynamics over decadal timescales, we neglect simulation of dreissenid mussels and zooplankton. Following Snortheim et al. [33], we have subsumed the associated mortality and nutrient cycling into the respiration parameter, allowing for regulation without enabling a zooplankton and mussel functional group [35]. This is justified by the efficient internal recycling of nutrients modelled to occur in Lake Erie through predation/excretion [36]. The inclusion of these processes could be the subject of future work.

Initial and Boundary Conditions
The bathymetric profile (area vs. depth; Figure 1b) was manually computed from a 2 km × 2 km × 1 m Lake Erie grid (https://www.ngdc.noaa.gov/mgg/greatlakes/erie. html (last accessed on 2 July 2021)). The model was initialized on 1 May 1979, using available water quality data from the Maumee River, as there were no data available for western Lake Erie (Table 3), and advanced using an hourly timestep. Mean daily meteorological forcing data included air temperature, wind speed, relative humidity, precipitation, shortwave solar radiation, and cloud cover. Daily solar radiation data were disaggregated into sub-daily (rad mode = 0), with longwave radiation computed from cloud cover (lw_type = 'LW_CC'; cloud mode = 1). Mean daily boundary conditions were specified for the Detroit River and Maumee River (flow, temperature, and the water quality state variables; Tables S1 and S2). Exchange to the central basin was specified as an outflow, based on a water balance computed from the Detroit and Maumee River flows, precipitation, evaporation, and observed water levels [37,38]. Default parameters were used for ice and snow (snow_albedo_factor = 1.0; snow_rho_max = 300; snow_rho_min = 50).

Model Calibration
To minimize uncertainty in parameter estimation and avoid user bias in model calibration [20], GLM-AED was calibrated using model-independent parameter estimation and uncertainty analysis (PEST; http://www.pesthomepage.org/ (last accessed on 2 July 2021)). This approach is similar to previous studies that have applied autocalibration methods (Monte Carlo and PEST, respectively) to calibrate the 1D models DYRESM-CAEDYM [39] and Simstrat [24]. To apply PEST to all~60 model parameters would take 10 23 years, and so a sensitivity analysis was employed to determine which parameters required calibration, and the associated calibration ranges. Sensitivity of modelled water temperature, TP, PO 4 , total Chl-a, and the four phytoplankton groups to changes in parameter values were evaluated according to relative sensitivity (RS; [23,40]): where i is the modelled value (output), j is the calibrated parameter value (input), ∆C i is the change in the modelled value, C is is the initial modelled value, ∆β j is the change in the parameter, and β js is the initial parameter value. In this study, to cross-compare RS values, absolute relative sensitivity (ARS) was used, which is the absolute value of RS. If ARS = 0.5, a 10% increase or decrease in the model parameter will cause a 5% change in the modelled state variable. The model parameters with ARS values (Table 4) were calibrated with PEST over parameter ranges based on literature values (Tables S3-S5). PEST optimizes the goodness-of-fit (minimizes root-mean-square error, RMSE) between model output and observation. Datasets for calibration were chosen that most comprehensively covered the entire western basin based on Table 2, Figure 1. In practice, PEST was first applied to calibrate temperature, using GLM, to achieve the smallest RMSE between simulations and observations. Subsequently, PEST was applied to AED-GLM to calibrate nutrients (TP and PO 4 ). Finally, total Chl-a and the four phytoplankton groups were calibrated to ensure the model can reproduce seasonal succession. This iterative approach was favored over attempts to use normalized RMSE, which can give misleading results when RMSE is small [22].  Notes: K r -phytoplankton respiration/metabolic loss rate of 20 • C (d −1 ); P max -maximum phytoplankton growth rate of 20 • C (d −1 ); T std -standard temperature for algal growth ( • C); T opt -optimum temperature for algal growth ( • C); T max -maximum temperature for algal growth ( • C); V T -Arrhenius temp scaling coefficient for growth; K W -extinction coefficient for PAR radiation (m −1 ); Wind_ f actorwind factor; Lw_ f actor-longwave factor; C e -bulk aerodynamic coefficient for latent heat transfer; C h -bulk aerodynamic coefficient for sensible heat transfer;

Phosphorus Loading Reduction Scenarios
The Maumee River delivers significant P to the western basin and has an agricultural watershed that has been the study of the effectiveness of BMPs on nutrient load reductions (e.g., [14,41]). Therefore, to explore the impacts of phosphorus loading reduction scenarios on water quality changes, the nutrient loads from the Maumee River were varied in two management scenarios. In the first scenario, the observed flow rates and nutrient concentrations were scaled according to load reductions realized by implementing BMPs in a published application of the SWAT model to the Maumee River watersheds over 1998-2005 [14]; their BMPs were at a moderate level considered feasible by agricultural specialists. In the second scenario, observed nutrient loads were reduced by 20, 40, and 60% to test if these reductions will achieve the target of limiting PO 4 and TP loads from the Maumee River to 186 MTA and 860 MTA during March to July, and limiting the maximum 30-day average cyanobacteria biomass to 9600 MTA [42].

Temperature, Water Levels, and Ice Cover
Accurate modelling of phosphorus loads and algae growth requires simulation of the water balance [43] and water temperature/ice cover [44], respectively. The time-series of surface water temperature indicated that the model reproduced the annual variations of water temperature, which varied between −5 • C (in winter) and 25 • C (in summer) ( Figure 2). The model also reproduced the seasonal stratification profile, including ephemeral weak stratification ( Figure 3). The surface RMSE was 2.95 • C (1979 to 2015), and vertically averaged RMSE values were 1.39 • C in 1994 and 1.70 • C in 2008 ( Figure S1). These RMSE values were consistent with those from both 1D (1-6 • C; [45]) and 3D (1-3 • C; [46]) models applied to Lake Erie, and Bruce et al. [47] who applied GLM to 32 lakes (RMSE of 1.62 to 1.31 • C through the water column).  Simulated water levels had an RMSE of 0.36 m, in comparison to observations (average of four western basin gauge stations: Buffalo, Cleveland, Port Stanley, and Toledo). These data were daily averages to remove the effects of the 14-hr surface seiche, which was not resolved by the horizontally averaged model. Given that the outflow was determined from a water balance, which included the observed levels, the RMSE was directly attributable to differences between the evaporation and precipitation models employed in the water balance versus GLM. The simulated ice thickness (not shown) had an RMSE of 0.13 m, compared with observations [48] in Sandusky Bay, which is just outside the model domain.

Nutrients
Simulated TP was visually consistent with annually and seasonally averaged observations from four published studies over 1979 through 2015 ( Figure 4). From the early 1980s to the early 1990s, both simulated and observed TP concentrations showed a decreasing trend, coincident with the implementation of point-source P abatement programs ( Figure 4a). Moreover, from the mid-1990s, from May to September, TP concentrations (Figure 4b) increased in accordance with P from the Maumee River watersheds being flushed by increased spring precipitation. Peaks in observed TP concentration occurred in April, consistent with the spring P loads from the Maumee River, and were the major P source for western Lake Erie; the simulations did not capture all peaks, potentially because loads are instantaneously distributed throughout the basin in the horizontally averaged model (Figure 4c,d). The average TP RMSE was 0.01 mg L −1 , comparable to 0.01-0.074 mg L −1 for a 1D model of Lake Ravn [32] and 0.03 mg L −1 for a 3D model of Lake Erie [49].

Nutrients
Simulated TP was visually consistent with annually and seasonally averaged observations from four published studies over 1979 through 2015 ( Figure 4). From the early 1980s to the early 1990s, both simulated and observed TP concentrations showed a decreasing trend, coincident with the implementation of point-source P abatement programs ( Figure 4a). Moreover, from the mid-1990s, from May to September, TP concentrations ( Figure 4b) increased in accordance with P from the Maumee River watersheds being flushed by increased spring precipitation. Peaks in observed TP concentration occurred in April, consistent with the spring P loads from the Maumee River, and were the major P source for western Lake Erie; the simulations did not capture all peaks, potentially be- cause loads are instantaneously distributed throughout the basin in the horizontally averaged model (Figure 4c,d). The average TP RMSE was 0.01 mg L −1 , comparable to 0.01-0.074 mg L −1 for a 1D model of Lake Ravn [32] and 0.03 mg L −1 for a 3D model of Lake Erie [49]. Depth-averaged PO4 concentrations were both simulated and observed to be ~0.02 mg L −1 during 2000 through 2014 (Figure 5a). However, the simulations were, at times, smaller than the observations at the stations located 1 m above the lake bed (Figure 5b), suggesting mineralization of PO4, which has been shown to vary spatially by an order of magnitude in western Lake Erie [26], was not always accurately reproduced with the static release model in AED. Sediment mineralization of PO4 was the only phosphorus boundary condition that was not directly measured. The parameters regulating sediment mineralization ( ℎ _ _ , _ , and _ ; Table 4) were automatically calibrated; therefore, simulation of less PO4 from the sediments suggests that PEST may be compensating for excessive PO4 loads, either from the tributaries or internal cycling.
During the late 1980s through 2010, the PO4 concentration increased, with increases in soluble P from the Maumee River watershed [50], which is in agreement with the simulations, showing PO4 to be maximal during 2010 through 2014. The average RMSE of PO4 was 0.01 mg L −1 , comparable to 0.007-0.061 mg L −1 for simulations of Lake Raven [32] using 1D DYRESM-CAEDYM.  Depth-averaged PO 4 concentrations were both simulated and observed to be~0.02 mg L −1 during 2000 through 2014 (Figure 5a). However, the simulations were, at times, smaller than the observations at the stations located 1 m above the lake bed (Figure 5b), suggesting mineralization of PO 4 , which has been shown to vary spatially by an order of magnitude in western Lake Erie [26], was not always accurately reproduced with the static release model in AED. Sediment mineralization of PO 4 was the only phosphorus boundary condition that was not directly measured. The parameters regulating sediment mineralization (Theta_sed_ f rp, Fsed_ f rp, and Ksed_ f rp; Table 4) were automatically calibrated; therefore, simulation of less PO 4 from the sediments suggests that PEST may be compensating for excessive PO 4 loads, either from the tributaries or internal cycling. The annual and seasonal variations in simulated and sampled NO3 and NH4 were small (Figures 6 and 7). During early summer, the predominant algae were non-N-fixing, and in late summer-if required-they were able to fix nitrogen directly from the atmosphere [51]. The average RMSE for NO3 and NH4 was 0.68 and 0.03 mg L −1 , respectively, comparable to 0.036 mg L −1 for NO3 from 3D simulations of Lake Erie [49]. During the late 1980s through 2010, the PO 4 concentration increased, with increases in soluble P from the Maumee River watershed [50], which is in agreement with the simulations, showing PO 4 to be maximal during 2010 through 2014. The average RMSE of PO 4 was 0.01 mg L −1 , comparable to 0.007-0.061 mg L −1 for simulations of Lake Raven [32] using 1D DYRESM-CAEDYM.
The annual and seasonal variations in simulated and sampled NO 3 and NH 4 were small (Figures 6 and 7). During early summer, the predominant algae were non-N-fixing, and in late summer-if required-they were able to fix nitrogen directly from the atmosphere [51]. The average RMSE for NO 3 and NH 4 was 0.68 and 0.03 mg L −1 , respectively, comparable to 0.036 mg L −1 for NO 3 from 3D simulations of Lake Erie [49].
(black asterisks in Figure 1) for 1986-2015. Data sources: The annual and seasonal variations in simulated and sampled NO3 and NH4 were small (Figures 6 and 7). During early summer, the predominant algae were non-N-fixing, and in late summer-if required-they were able to fix nitrogen directly from the atmosphere [51]. The average RMSE for NO3 and NH4 was 0.68 and 0.03 mg L −1 , respectively, comparable to 0.036 mg L −1 for NO3 from 3D simulations of Lake Erie [49].  Table 2.

Phytoplankton
The model reproduced an increase in total Chl-a from 1994 through 2015 (Figure 8). The Chl-a RMSE between simulations and ECCC observations (17.85 μg L −1 ) was smaller than that between simulations and GLENDA observations (19.21 μg L −1 based on Reavie

Phytoplankton
The model reproduced an increase in total Chl-a from 1994 through 2015 (Figure 8). The Chl-a RMSE between simulations and ECCC observations (17.85 µg L −1 ) was smaller than that between simulations and GLENDA observations (19.21 µg L −1 based on Reavie et al. [30], and 6.25 µg L −1 based on Elbagoury [21]) ( Figure 8). The error was comparable to 10.4 and 12.76 µg L −1 in Lake Ravn in Denmark using 1D DYRESM-CAEDYM [32]. The ECCC data were sampled through the water column across the whole western basin and the vertically averaged simulations from the horizontally averaged model capture these variations. GLENDA samples were measured through the water column during spring and in the epilimnion (from 2 to 4 m depths) during summer at only six stations (ER58, 59, 60, 61, 91, and 92; Figure 1). Particularly, ER59 was located near the Maumee River mouth, along the southwest border in the western basin, resulting in the observations at ER59 during summer being larger than the simulations near the Maumee River plume [52]. study represent early diatoms (optimum temperature of 9.8 °C), with high Si requirements and rapid sinking rates. Like cyanobacteria, the peak of observed diatoms also occurred in 2008 and 2011, which was not reproduced by the model. Peaks in both simulated and observed cyanobacteria were in 2008 and 2011, corresponding to the lowest concentrations of greens and cryptophytes, suggesting the cyanobacteria are out-competing other groups in these years. In general, chlorophytes and diatoms are the preferred food for grazers (zooplankton and fish), with cyanobacteria being less grazed upon.
In addition to the GLENDA cyanobacteria observations, the modeled maximum rolling 30-day average of cyanobacteria concentration was compared to the NOAA bloom severity index (cyanobacteria index, CI; Figure 10). The CI was determined over 10-day intervals from remote sensing, taking the highest cyanobacterial chlorophyll-related index at each pixel available from any of the daily images within a 10-day period. The algal bloom severity was determined from the annual CI [53][54][55], which was the average of the 10-day intervals around the maximum severity of the bloom [56]. Intense algal blooms typically lasted 30-40 days [55]. When the index was in the range of 2-4, blooms were regarded as mild, with severe blooms above 4. Both the highest severity index and simulated maximum 30 [30] and Elbagoury [21]. Data sources: Table 2.  [17] (yellow dots in Figure 1) for 2008, and 2011-2014 (Sta. MB18, 8MGR1). MB18 is closest to the Maumee River mouth. The GLENDA data have been converted from biovolume to biomass using the conversions in both Reavie et al. [30] and Elbagoury [21]. Data sources: Table 2.
However, being spatially averaged, gradients in Chl-a from near the Maumee River mouth to offshore were observed, with higher than observed concentrations at MB18 and 8M, and lower concentrations at GR1 (Figure 8 [30] and Elbagoury [21]. Data sources: Table 2.  The GLENDA data have been converted from biovolume to biomass using the conversions in both Reavie et al. [30] and Elbagoury [21]. Data sources: Table 2. In addition to the GLENDA cyanobacteria observations, the modeled maximum rolling 30-day average of cyanobacteria concentration was compared to the NOAA bloom severity index (cyanobacteria index, CI; Figure 10). The CI was determined over 10-day intervals from remote sensing, taking the highest cyanobacterial chlorophyll-related index at each pixel available from any of the daily images within a 10-day period. The algal bloom severity was determined from the annual CI [53][54][55], which was the average of the 10-day intervals around the maximum severity of the bloom [56]. Intense algal blooms typically lasted 30-40 days [55]. When the index was in the range of 2-4, blooms were regarded as mild, with severe blooms above 4. Both the highest severity index and simulated maximum  [30] and Elbagoury [21]. Data sources: Table 2. Figure 10. Observed western Lake Erie bloom severity index [55] compared to simulated maximum 30-day average cyanobacteria concentration. Figure 10. Observed western Lake Erie bloom severity index [55] compared to simulated maximum 30-day average cyanobacteria concentration.

Discussion
The physical and biogeochemical state variables were simulated with similar RMSE to those from shorter-term simulations as reported in the literature. This provides confidence in model-output over the entire 1980 to 2015 simulation ( Figure 11) and in the ability of the model to evaluate nutrient reduction management scenarios. Model results show interannual variation in simulated water temperature, but a long-term trend in temperature was not evident. Both PO 4 and TP have spikes in concentration that are visually more evident during the early 1980s and after 2000 (Figure 11). Similarly, total Chl-a and cyanobacteria concentrations were elevated during these times; however, cyanobacteria were also simulated to increase during the late 1980s to early 1990s. The simulated longterm TP concentrations were typically lower than 0.02 mg L −1 , as expected for the target TP load of 11,000 MTA established by GLWQA in 1978 [57].

Discussion
The physical and biogeochemical state variables were simulated with similar RMSE to those from shorter-term simulations as reported in the literature. This provides confi-1981 1983 1985 1987 1989 1991 1993 1995 1997 1999 Figure 11. Long-term simulations of daily depth-averaged water temperature, PO4, TP, total Chl-a, and cyanobacteria concentrations in the western basin.
The Maumee River is the major phosphorus source for the western basin [1,50]. For example, during 2007 (1 October 2006-30 September 2007, the Maumee River delivered the largest phosphorus and nitrogen loads to the western basin in 33 years of monitoring [58]. Increases in TP and cyanobacteria were strongly related to the increased PO 4 loads from the Maumee River during 1996-2006 [9]. However, modelling results show that elevated PO 4 does not occur in the lake during bloom years (2008,2011,2013,2014, and 2015; Figure 11). We believe this is because all available PO 4 goes to phytoplankton uptake, causing low water-column PO 4 concentrations within the blooms. This inverse relationship has also been modeled in Lake Erie [19], where the lakewide highest Chl-a (~50 mg m −3 ) was associated with the lowest PO 4 (~5 mg m −3 ). Hence, to give insight into the correlation between algal blooms and nutrients in the western basin, researchers must consider influent PO 4 loads, as opposed to in situ concentrations.
Maumee River PO 4 loads were primarily delivered through spring runoff. Due to the discrete nature of storm events, both discharge and PO 4 load were integrated (Figure 12a,b), following Richards et al. [58]. The cumulative discharge and PO 4 load were largest in 2015 and 2011, which correspond to the bloom years with the highest simulated cyanobacteria concentration (Figure 11). High discharge leads to high loads, consistent with Baker et al. [59], and the 2015 discharge was nearly two times higher than the 35-year average. The hydrological water year begins in March, whereas, in Richards et al. [58], the water year began in October, resulting in the largest flow, TP, and PO 4 loads in 2007 [58], which was not a bloom year (Figure 10), thereby suggesting PO 4 loads in spring are more impactful on summer bloom severity, compared to those during the preceding fall and winter.

Nutrient-Reduction Scenarios
An abundance of literature has recommended that P mitigation can control eutrophication in western Lake Erie [14,17,63]. The USEPA [12] established the target of limiting March-July dissolved reactive phosphorus loading from the Maumee River to 186 metric tonnes and total phosphorus loading to 860 metric tonnes, an approximately 40% reduction from 2008 loads (closest to the original 1978 Annex 3 target of 11,000 MT), in order to reach 'mild bloom' threshold of 9600 MTA [42]. Nutrient loads may be controlled by BMPs (best management practices) in the Maumee River watershed to reduce the P loads entering into the western basin [14,41]. However, the direct effect of the loads recommended in the Annex 4 targets, as well as the direct effect of achievable loads from the implementation of BMPs, has not been simulated over decadal timescales.
The 17,030 km 2 Maumee River watershed receives 934 mm yr −1 precipitation and is In terms of algal biovolume, bloom area, and duration, the bloom area in 2011 [55] was more than two times higher than that in 2008, and about four times greater than that from 2002-2010 [11,60]. The blooms in 2011 and 2015 were similar, while the PO 4 loads in 2011 were much lower than in 2015, suggesting other factors, such as air temperature, also played an important role in bloom development [61,62].
We compared the monthly average air temperature in 2011 and 2015 to the 1981-2010 climate normal (Figure 12c). During summer (Jul.-Sep.), when cyanobacteria growth was rapid, the air temperature in 2011 was 1.2 • C higher than the 1981-2010 climate normal, and, in July and August 2011, the air was close to the optimal growth temperature for cyanobacteria (24 • C) leading to a strong bloom. The average air temperature in late summer (Jul.-Sep.) of 2015 was 21.8 • C, lower than that in 2011 (22.4 • C) but higher than 1981-2010 (21.2 • C). Since the PO 4 loads in 2015 were much higher than in 2011, the blooms in 2015 were still high, similar to 2011. The interplay between loads and temperature ( Figure 12d) shows modelled cyanobacteria to peak in years when both air temperature and nutrient loads (Figure 12b) are high (2011 and 2015). Warm air/water or high loads, in isolation, were insufficient (e.g., 2010 or 2012).

Nutrient-Reduction Scenarios
An abundance of literature has recommended that P mitigation can control eutrophication in western Lake Erie [14,17,63]. The USEPA [12] established the target of limiting March-July dissolved reactive phosphorus loading from the Maumee River to 186 metric tonnes and total phosphorus loading to 860 metric tonnes, an approximately 40% reduction from 2008 loads (closest to the original 1978 Annex 3 target of 11,000 MT), in order to reach 'mild bloom' threshold of 9600 MTA [42]. Nutrient loads may be controlled by BMPs (best management practices) in the Maumee River watershed to reduce the P loads entering into the western basin [14,41]. However, the direct effect of the loads recommended in the Annex 4 targets, as well as the direct effect of achievable loads from the implementation of BMPs, has not been simulated over decadal timescales.
The 17,030 km 2 Maumee River watershed receives 934 mm yr −1 precipitation and is comprised of 76% row-crop, 11% urban, 8% forested, and 5% hay [14]. We parameterized the effects of reducing tillage, planting cover crops, and the addition of edge-of-field filter strips (Table 5) by scaling our Maumee River loads to match the changes given in the Bosch et al. [14] reduction scenarios. This was preferable to directly applying the Bosch et al. modelled loads, as they differed from observed loads by 3-4%. When BMPs were considered, the simulated maximum 30-day average cyanobacteria concentrations, PO 4 , and TP loads were smaller than those under current conditions, causing the western Lake Erie cyanobacteria to be lowest under the source-combined scenario (Figure 13a,c,e). For the different BMP scenarios, the cyanobacteria biomass was reduced by 1.98-12.21%. No-till was the least effective among individual BMPs (Table 6), potentially facilitating fertilizer accumulation in the surface soil layer, which is then flushed into the lake by runoff [64]. The source-combined BMP scenario reduced TP, PO 4 , TN, NO 3 , Chl-a, and CYANO average annual concentrations by 3.6%, 1.7%, 2.0%, 2.0%, 0.68%, and 12.2%, respectively ( Table 6). The relatively small change in nutrients arose because the source-combined scenario only reduced PO 4 and TP loads by 7% and 8%, respectively. Implementing BMPs to 100% of all row-crop land in the Maumee watershed above the moderate BMP level considered feasible gave a reduction in P yield of 30% [14]. Therefore, a 10-20% BMP reduction in P yield seems reasonably achievable but is considerably less than the Annex 4 recommended P reduction of 40%.
To test the effect of 20%, 40%, and 60% P load reduction scenarios, the TP, PO 4 , DOP, and POP concentrations in the Maumee River were all reduced accordingly (Figure 13b,d,f). The resultant simulated maximum 30-day average cyanobacteria biomass indicated that significant P reductions in the Maumee River can lead to decreases in cyanobacteria biomass by limiting growth through P limitation [65]. When the TP and PO 4 spring loads were reduced by 40% (849 and 181 MTA in 2008, respectively), the maximum 30day average cyanobacteria biomass was below the 'mild bloom' threshold of 9600 MTA during most years (from 2005 to 2015) (Figure 13b,d,f). The 20% reduction scenario, which approximates the presently feasible maximum BMP impact, shows the maximum 30-day average cyanobacteria biomass was to be above 9600 MTA continuously from 2008 through the end of the simulation, having an insufficient effect on severity of the blooms during this time ( Figure 10).    Table 6. Average annual TP, PO4, TN, NO3, Chl-a, and CYANO concentrations in western Lake Erie for various BMP  These simulated results are consistent with the Great Lakes Water Quality Agreement Annex 4 recommendation that a mild bloom can be achieved by limiting DRP and TP loads from the Maumee River during March and July to 186 and 860 MTA, respectively, which was approximately a 40% reduction from 2008 values [66]. Zhang et al. [10] explored the impacts of reducing TP and PO 4 loads (by reducing concentrations) from the Maumee River by 20%, 40%, 60%, and 80% on the changes in P pools and fluxes, as well as algal biomass in western Lake Erie. The 80% P load reduction gave 65 and 88% decreases in total algal biomass in late June of 1997 and from June to October in 1998. Similar work by Verhamme et al. [17] found that lower Maumee River March-July TP loads of approximately 400 MTA were necessary to achieve a bloom with a maximum 30-day average biomass equivalent to the bloom in 2012.
Our 1D model indeed shows P to be the limiting nutrient from May through December, which is in agreement with observations of P limitation in the north of the western basin (Pigeon Bay, [67]). Our 1D model was, however, unable to capture spatial gradients in nutrient limitation. For example, near the mouth of the Maumee River, phytoplankton growth has been observed to be P-replete during wet years, but with low or no correspondence between nitrogen limitation and size of the cyanobacterial bloom [68].
Similar to our findings, application of the 3D model ELCOM-CAEDYM to Edmonton (Canada) stormwater ponds Nakhaei [22] simulated that a reduction of influent P and N fractions by at least 50% was required to improve the trophic state of each pond from mesotrophic/eutrophic to oligotrophic/mesotrophic. Simulation of Lake Raven (Denmark; [32]) with 1D DYRESM-CAEDYM also predicted that a substantial (40-50%) reduction in external TP loading would be required to meet phytoplankton biomass requirements. However, some model applications have suggested N control may also be necessary. Application of 0D WASP to Lake Winnipeg (Canada; [40]) showed that a 10% reduction in P loads decreased cyanobacteria and peak chlorophyll-a concentrations but promoted growth of non-N-fixing cyanobacteria. They modelled that increasing N:P loading ratio (P reduction > 12% and N reduction < 7%) would be effective for improving water quality in the lake. Modelling of Lake Okaro (New Zealand; [69]) suggests that N loading reduction may be more successful than reducing P loading alone, because N was modelled to be a major limiting factor for cyanobacteria growth. Given the observations of N limitation near the Maumee River [68], future modelling exercises should assess combined N and P reduction scenarios.

Conclusions
AED-GLM reproduced the long-term temperature and water quality in western Lake Erie without model drift, yielding RMSE that was comparable to single year simulation studies in the literature. The algal blooms in the early 1990s and recent years (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) were simulated, without a need to reinitialize or recalibrate in individual years. In agreement with Annex 4, a 40% reduction in P loads from the Maumee River was necessary to achieve a 'mild bloom' during most years. However, by scaling Maumee River concentrations to account for the documented effect of BMPs, the achievable reductions of 10-20% in PO 4 and TP loads would not be sufficient to achieve the bloom reduction goal in Annex 4 for western Lake Erie. We recommend future modelling efforts link long-term simulations directly to output from watershed models and consider combined N and P reduction scenarios.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/su13147516/s1, Figure S1: Calculated RMSE temperature at Sta. W1,2 for 1994 (a) and STN357 for 2008 (b), from time series shown in Figure 3, Table S1: Summary of Water Quality Data Sources for the Detroit River, Table S2: Summary of Water Quality Data Sources for the Maumee River, Table S3: Description, default and assigned values of the parameters in GLM, Table S4: Description, default and assigned values of the parameters in AED, Table S5: Description, default and assigned values of phytoplankton parameters in AED.