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

Abstract

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⁻¹), PO₄ (0.01 mg L⁻¹), NH₄ (0.03 mg L⁻¹), NO₃ (0.68 mg L⁻¹), total chlorophyll a (18.74 μg L⁻¹), chlorophytes (3.94 μg L⁻¹), cyanobacteria (12.44 μg L⁻¹), diatoms (3.17 μg L⁻¹), and cryptophytes (3.18 μg L⁻¹) 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.

1. 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 future trends in agricultural practices and climate-driven meteorological forcing, which favor bloom development with significant interannual variation in severity [11], suggesting algae blooms will continue and potentially worsen in the future.

The USEPA [12] report recommended a total P (TP) spring load of 860 metric tons and a dissolved reactive phosphorus (DRP) load of 186 metric tons from the Maumee River to achieve a western-basin bloom (90% of the time) no greater than that observed in 2004 or 2012. The 860 metric ton target is approximately a 40% reduction of the 2008 spring load, as determined from modeling studies. The DRP reduction addressed the shift from point-source to agriculturally derived non-point-source P loading. To assess the ability of this reduction to be realized in the Lake Erie watershed, Makarewicz et al. and Bosch et al. [13,14] modelled implementation of agricultural best management practices (BMPs; e.g., reduced tillage, cover crops, and filter strips) to control P loads within the Maumee River. However, the direct impacts of BMPs on water quality and nuisance bloom formation within western Lake Erie remain uninvestigated.

The more recent 2016 GLWQA Annex 4 applied a suite of models to relate P loads to overall phytoplankton biomass, occurrence of cyanobacteria blooms, degree of hypoxia, and presence of Cladophora. Focusing again on 2008, a 50% reduction in maximum western basin chlorophyll a (Chl-a) concentration was achieved by reducing TP loads from the Maumee River (in particular, P from March to July storm events). Particulate P was a significant fraction of the bioavailable P pool [15], and so reducing the Maumee DRP load alone would be insufficient to prevent bloom development. A sensitivity analysis confirmed that P from the Detroit River was not a driver of cyanobacteria blooms. Overall, a mean annual western basin P load of 2193 MTA was recommended using a combination of 2D (EcoLE, [10]), 3D (ELCOM–CAEDYM, [16], WLEEM, [17]), and semiempirical mass balance models [18]. WLEEM investigated the relationship between P loads and cyanobacteria response over March through November of 2008 and 2011 to 2014. They modelled that 890 metric tons of total P from the Maumee River over March to July produced a ‘mild’ cyanobacterial bloom in western Lake Erie [17].

However, these models had several limitations that make it difficult to delineate the impacts of long-term trends in climate-driven meteorological forcing from agricultural practices on the relationships between nutrients and nuisance blooms in western Lake Erie. The models were reinitialized annually each spring, making it (1) impossible to capture the cumulative effects of changes in nutrient loads over years or decades, and (2) difficult to account for changes in nutrient loads on model initialization. Typically, the initial nutrient concentrations are not changed within model specified initial conditions, causing the lake to respond slowly to load reductions over the residence time (~2 months for the western basin), which overlaps the spring bloom. Moreover, (3) nutrient loads are simply scaled (e.g., [10,16,17]), and so the effectiveness of BMPs has not been directly tested.

To address these issues related to seasonal simulations with static initial conditions, there is a need to run long-term models (much longer than the residence time and sufficient to resolve decadal changes in agricultural practices and climate forcing). However, 3D hydrodynamic–biogeochemical lake models are often under-calibrated due to the long run-times (e.g., ~3 weeks for a 6 month simulation of Lake Erie with ELCOM–CAEDYM; [19]). These uncertainties in model parameters and initial conditions are the primary model input errors [20], which can result in poor representation of processes, such as seasonal phytoplankton succession (e.g., [21,22]). In the present study, the one-dimensional vertical coupled hydrodynamic biogeochemical Aquatic Ecosystem Dynamics–General Lake Model (AED–GLM; [23]) has been applied to simulate western Lake Erie from 1979 through 2015, forced with both historical observed and modelled BMP tributary loads. GLM–AED runs significantly faster (~1 min for 36 years) compared to 3D models, allowing for long-term simulations, in-depth sensitivity analysis, and the utilization of an automated calibration procedure to optimize parameter settings (e.g., [24]).

The objectives of the present study were to: (1) calibrate AED–GLM using model-independent parameter estimation; (2) apply the model to simulate the interannual variation 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. List of acronyms.

HFCB	harmful cyanobacteria
AED	Aquatic Ecosystem Dynamics
GLM	General Lake Model
P	phosphorus
N	nitrogen
PEST	Model-Independent Parameter Estimation
GLWQA	Great Lakes Water Quality Agreement
TP	total phosphorus
DRP	dissolved reactive phosphorus
USEPA	United States Environmental Protection Agency
BMP	agricultural best management practice
Chl-a	total chlorophyll a
MTA	metric tons per annum
ELCOM	Estuary and Lake Computer Model
CAEDYM	Computational Aquatic Ecosystem Model
WLEEM	Western Lake Erie Ecosystem Model
EcoLE	Ecological Model of Lake Erie
GLERL	Great Lakes Environmental Research Laboratory
NOAA	National Oceanic and Atmospheric Administration
ECCC	Environment and Climate Change Canada
EMRB	Environmental Monitoring and Reporting Branch
OCWA	Ontario Clean Water Agency
GLENDA	Great Lakes Environmental Database
NDBC	National Data Buoy Centre
GREEN	green algae
CYANO	cyanobacteria
DIAT	diatoms
CRYPT	cryptophytes
DO	dissolved oxygen
PO4	dissolved reactive phosphorus
NO3	nitrate
NH4	ammonium
RSi	reactive silica
DON	dissolved organic nitrogen
DOP	dissolved organic phosphorus
DOC	dissolved organic carbon
PON	particulate organic nitrogen
POP	particulate organic phosphorus
POC	particulate organic carbon
RS	relative sensitivity
ARS	absolute relative sensitivity
RMSE	root-mean-square error
DYRESM	Dynamic Reservoir Simulation Model
SWAT	Soil and Water Assessment Tool
CI	cyanobacteria index
WASP	Water Quality Analysis Simulation Program

.

2. Methods

2.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², a volume of ~20 km³, 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.

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).

To compare the modeled phytoplankton concentrations to observations, the observed GLENDA biovolumes (µm³ L⁻¹) were converted to Chl-a biomass (μg L⁻¹) using two conversions: (1) we applied the formula ${\log }_{10}\left(\mathit{biovolume}\right)=a+b{\log }_{10}\left(\mathrm{chlorophyll} a\right)$ with species-specific values of a and b [30], and (2) the biovolumes were multiplied by density (1g cm⁻³) and 7.5% of the phytoplankton biomass used to estimate Phyto C, followed by a Phyto C: Chl-a = 50:1 (μg L⁻¹) ratio being applied [21].

2.2. 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₄), nitrate (NO₃), ammonium (NH₄), 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⁻¹) and mineralization of dissolved organics were modeled as migration with photoinhibition (Rdoc_minerl = 0.001 day⁻¹; Rdon_minerl = 0.005 day⁻¹; Rdop_minerl = 0.001 day⁻¹). 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.

2.3. 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. Initial biogeochemical conditions specified in the model for 1 May 1979. For data sources see Table S2.

DO (mmol O2/m3)	SiO2 (mmol Si/m3)	NH4 (mmol N/m3)	NO3 (mmol N/m3)	PO4 (mmol P/m3)	PON (mmol N/m3)	DON (mmol N/m3)	POP (mmol P/m3)
376.76	114.29	9.66	47.42	2.97	30.18	80.88	3.48
DOP (mmol P/m3)	POC (mmol C/m3)	DOC (mmol C/m3)	GREEN (mmol C/m3)	DIAT (mmol C/m3)	CYANO (mmol C/m3)	CRYPT (mmol C/m3)
2.32	41.67	250	1.59	1.81	0	0.03

), 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).

2.4. 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²³ 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₄, total Chl-a, and the four phytoplankton groups to changes in parameter values were evaluated according to relative sensitivity (RS; [23,40]):

$$RS= \frac{\raisebox{1ex}{$\Delta C_i$}\!\left/ \!\raisebox{-1ex}{$C_{is}$}\right.}{\raisebox{1ex}{$\Delta {\beta }_j$}\!\left/ \!\raisebox{-1ex}{${\beta }_{js}$}\right.}$$

(1)

where i is the modelled value (output), j is the calibrated parameter value (input), $\Delta C_i$ is the change in the modelled value, $C_{is}$ is the initial modelled value, $\Delta {\beta }_j$ is the change in the parameter, and ${\beta }_{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. Summary of absolute relative sensitivity (ARS) values from the sensitivity analysis. The model parameters with the 6 largest ARS values (shaded) were calibrated with PEST [40]. The ‘--’ indicates no sensitivity.

Model Parameter	ARS of Modelled State Variable
	Water Temperature	PO4	TP	Chl-a	GREEN	CYANO	DIAT	CRYPT
$K_w$	0.02	0.37	0.09	0.86	4.89	3.19	0.58	4.6
$Wind\_factor$	0.06	0.28	0.07	0.35	0.66	0.93	0.66	0.22
$Lw\_factor$	0.66	0.32	0.24	0.78	2.62	5.50	2.58	12.88
$C_e$	0.01	--	--	--	--	--	--	--
$C_h$	0.01	--	--	--	--	--	--	--
$C_d$	0.01	--	--	--	--	--	--	--
$Theta\_sed\_frp$	--	0.65	0.07	--	--	--	--	--
$Fsed\_frp$	--	0.29	0.15	0.31	1.25	2.20	0.02	1.06
$Ksed\_frp$	--	0.25	0.13	0.26	1.06	1.87	0.03	0.91
$P_{max, GREEN}$	--	0.24	0.06	0.75	20.70	7.69	2.49	8.29
$V_{T, GREEN}$	--	0.03	0.03	0.30	9.80	4.98	0.06	8.44
$T_{std, GREEN}$	--	0.08	0.05	0.60	14.44	7.94	0.29	7.31
$T_{opt, GREEN}$	--	0.05	0.04	0.45	11.89	6.58	0.21	6.40
$T_{max, GREEN}$	--	0.01	0.01	0.16	9.67	2.58	0.04	2.90
$K_{r, GREEN}$	--	0.02	0.03	0.30	9.50	4.80	0.03	8.19
$P_{max, CYANO}$	--	0.23	0.04	0.72	9.39	12.94	2.48	8.60
$V_{T, CYANO}$	--	0.03	0.04	0.52	8.31	9.53	0.05	15.98
$T_{std, CYANO}$	--	0.004	0.03	0.35	8.18	8.49	0.02	6.89
$T_{opt, CYANO}$	--	0.004	0.03	0.30	6.53	8.89	0.03	5.39
$T_{max, CYANO}$	--	0.001	0.01	0.07	1.36	8.26	0.02	1.36
$K_{r, CYANO}$	--	0.002	0.04	0.45	7.87	8.96	0.13	15.10
$P_{max, DIAT}$	--	0.13	0.02	0.65	0.27	0.88	3.50	0.24
$V_{T, DIAT}$	--	0.001	0.06	0.05	1.76	2.69	9.93	0.96
$T_{std, DIAT}$	--	0.01	0.002	0.05	0.15	0.01	3.32	0.12
$T_{opt, DIAT}$	--	0.02	0.003	0.12	0.47	0.11	3.82	0.40
$T_{max, DIAT}$	--	0.10	0.01	0.38	1.93	0.57	2.94	1.84
$K_{r, DIAT}$	--	0.11	0.03	0.60	0.77	0.55	3.48	0.56
$P_{max, CRYPT}$	----	0.24	0.07	0.73	9.24	7.63	2.48	83.71
$V_{T, CRYPT}$	--	0.01	0.01	0.06	1.44	1.09	0.03	16.65
$T_{std, CRYPT}$	--	0.07	0.05	0.58	8.67	7.72	0.28	64.43
$T_{opt, CRYPT}$	--	0.05	0.04	0.47	8.21	6.72	0.21	57.72
$T_{max, CRYPT}$	--	0.008	0.01	0.12	2.92	1.91	0.03	17.90
$K_{r, CRYPT}$	--	0.004	0.01	0.07	1.38	1.01	0.001	16.98

Notes: $K_r$—phytoplankton respiration/metabolic loss rate of 20 °C (d⁻¹); $P_{max}$—maximum phytoplankton growth rate of 20 °C (d⁻¹); $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⁻¹); $Wind\_factor$—wind factor; $Lw\_factor$—longwave factor; $C_e$—bulk aerodynamic coefficient for latent heat transfer; $C_h$—bulk aerodynamic coefficient for sensible heat transfer; $C_d$—bulk aerodynamic coefficient for transfer of momentum; $Theta\_sed\_frp$—temperature multiplier for temperature dependence of sediment phosphate flux; $Fsed\_frp$—maximum flux of oxygen across the sediment water interface into the sediment (mmol P m⁻² d⁻¹); $Ksed\_frp$—half saturation constant for oxygen dependence of sediment phosphate flux (mmol O₂ m⁻³); $\mathrm{Water} \mathrm{temp}—$water temperature (° C); $\mathrm{TP}—$total phosphorus concentration (mmol m⁻³); $P{O}_4$—phosphate concentration (mmol m⁻³); Chl-a—total chlorophyll a (mmol m⁻³); GREEN—green algae (mmol m⁻³); CYANO—cyanobacteria (mmol m⁻³); DIAT—diatoms (mmol m⁻³); CRYPT—cryptophytes (mmol m⁻³).

) 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. Sources of measured data used in calibration (see also Figure 1). Shaded datasets had the most comprehensive spatial and temporal coverage and so were used for calibration with PEST.

Character	Location	Source	Identifier in Figure 1	Sample Date
Water temperature	Western basin surface and Sta. 45005	ECCC and NOAA & NDBC	Red and magenta dots	1979–2015 and 2005–2015
	Sta. W1, W2 and Sta. 357	Ackerman et al. (2001)	Blue circles and star	1994 and 2008
Nutrients (phosphorus and nitrogen)	Western basin	Ludsin et al. (2001)	Orange diamond	1980–1992
	Western basin	Thomas et al. (2014)	Black circles	1999–2015 (May–Sep.)
	1–11 m	ECCC	Green dots	1994–2015
	0–18.8 m	OCWA	Cyan dots	2001–2014
	1 m above bottom	EMRB	Black asterisks	1986–2015
	ER58 (integrated sample)	GLENDA	Black dots	1986–2015
	ER59 (integrated sample)
	ER60 (integrated sample)
	ER61 (integrated sample)
	ER91 (integrated sample)
	ER92 (integrated sample)
Chl-a	1–11m depth	ECCC	Green dots	1994–2015
	MB18 (integrated sample)	Verhamme et al. (2016)	Yellow dots	2008–2014
	8M (integrated sample)
	GR1 (integrated sample)
	ER58 (integrated sample)	GLENDA	Black dots	2001–2015
	ER59 (integrated sample)
	ER60 (integrated sample)
	ER61 (integrated sample)
	ER91 (integrated sample)
	ER92 (integrated sample)
Phytoplankton groups	ER58 (integrated sample)	GLENDA	Black dots	2001–2015
	ER59 (integrated sample)
	ER60 (integrated sample)
	ER61 (integrated sample)
	ER91 (integrated sample)
	ER92 (integrated sample)

, 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₄). 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].

2.5. 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₄ 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].

3. Results

3.1. 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.

3.2. 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⁻¹, comparable to 0.01–0.074 mg L⁻¹ for a 1D model of Lake Ravn [32] and 0.03 mg L⁻¹ for a 3D model of Lake Erie [49].

Depth-averaged PO₄ concentrations were both simulated and observed to be ~0.02 mg L⁻¹ 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₄, 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₄ was the only phosphorus boundary condition that was not directly measured. The parameters regulating sediment mineralization ($Theta\_sed\_frp$, $Fsed\_frp$, and $Ksed\_frp$; Table 4) were automatically calibrated; therefore, simulation of less PO₄ from the sediments suggests that PEST may be compensating for excessive PO₄ loads, either from the tributaries or internal cycling.

During the late 1980s through 2010, the PO₄ concentration increased, with increases in soluble P from the Maumee River watershed [50], which is in agreement with the simulations, showing PO₄ to be maximal during 2010 through 2014. The average RMSE of PO₄ was 0.01 mg L⁻¹, comparable to 0.007–0.061 mg L⁻¹ for simulations of Lake Raven [32] using 1D DYRESM–CAEDYM.

The annual and seasonal variations in simulated and sampled NO₃ and NH₄ were small (Figure 6 and Figure 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₃ and NH₄ was 0.68 and 0.03 mg L⁻¹, respectively, comparable to 0.036 mg L⁻¹ for NO₃ from 3D simulations of Lake Erie [49].

3.3. 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 $\mathsf{\mu }$g L⁻¹) was smaller than that between simulations and GLENDA observations (19.21 $\mathsf{\mu }$g L⁻¹ based on Reavie et al. [30], and 6.25 $\mathsf{\mu }$g L⁻¹ based on Elbagoury [21]) (Figure 8). The error was comparable to 10.4 and 12.76 $\mathsf{\mu }$g L⁻¹ 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].

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, transect running from Sta. MB18 to 8M and GR1) were not reproduced due to horizontal averaging in the 1D AED–GLM model framework. The model reproduced the observed variation at GR1 (RMSE = 17.24 $\mathsf{\mu }$g L⁻¹) near the mid-basin, while it underestimated the Chl-a concentration at MB18 (RMSE = 31.30 $\mathsf{\mu }$g L⁻¹) and 8M (RMSE = 20.57 $\mathsf{\mu }$g L⁻¹).

Peaks in both simulated and observed total Chl-a occurred in 2008 and 2011, with high cyanobacteria (Figure 9), and the lowest Chl-a appeared in 2009 and 2012, with low cyanobacteria. This is consistent with higher TP concentrations in 2008 and 2011, and lower values in 2009 and 2012 (Figure 4b,c). Both simulated and observed diatoms in this 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-day average cyanobacteria concentration occurred in 2011 and 2015, indicating these two years were with the most severe blooms. Overall, the simulations reproduced the increasing trend of blooms from 2002 to 2015 (Figure 10), as well as the maximum peaks in 2011 and 2015.

4. 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₄ 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 long-term TP concentrations were typically lower than 0.02 mg L⁻¹, as expected for the target TP load of 11,000 MTA established by GLWQA in 1978 [57].

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₄ loads from the Maumee River during 1996–2006 [9]. However, modelling results show that elevated PO₄ 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₄ goes to phytoplankton uptake, causing low water-column PO₄ concentrations within the blooms. This inverse relationship has also been modeled in Lake Erie [19], where the lakewide highest Chl-a (~50 mg m⁻³) was associated with the lowest PO₄ (~5 mg m⁻³). Hence, to give insight into the correlation between algal blooms and nutrients in the western basin, researchers must consider influent PO₄ loads, as opposed to in situ concentrations.

Maumee River PO₄ loads were primarily delivered through spring runoff. Due to the discrete nature of storm events, both discharge and PO₄ load were integrated (Figure 12a,b), following Richards et al. [58]. The cumulative discharge and PO₄ 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₄ loads in 2007 [58], which was not a bloom year (Figure 10), thereby suggesting PO₄ loads in spring are more impactful on summer bloom severity, compared to those during the preceding fall and winter.

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₄ 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₄ 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² Maumee River watershed receives 934 mm yr⁻¹ 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. Scenario descriptions for the BMPs applied to the Maumee River watershed. Our load reductions were obtained by scaling Maumee River P concentrations to match those achieved in Bosch et al. [14].

Scenario Name	Scenario Description	Change in TP Load (%)	Change PO4 Load (%)
No-till	No-till corn and soybean implemented in random 25% of row-crop agricultural land	+2.46	−1.54
Cover	Rye cover crop planted between soybean and corn crop in random 25% of row-crop agricultural land	−2.37	−1.66
Filter	Filter strip (10 m with 25% trapping efficiency) in random 20% of row-crop agricultural land	−2.97	−3.08
Random combined	Combination of three BMPs on same 25% of Maumee row-crop agricultural land; randomly distributed among sub-watersheds	−0.88	−4.41
Source combined	Combination of three BMPs on same 25% of Maumee row-crop agricultural land; distributed in high source sub-watersheds	−6.92	−8.04
Mouth combined	Combination of three BMPs on same 25% of Maumee row-crop agricultural land; distributed in sub-watersheds near river mouth	−0.86	−0.18

) 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₄, 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. Average annual TP, PO4, TN, NO3, Chl-a, and CYANO concentrations in western Lake Erie for various BMP implementation conditions from May 1979 to December 2015.

BMP Scenarios	TP (mmol m−3)	PO4 (mmol m−3)	TN (mmol m−3)	NO3 (mmol m−3)	Chl-a (μg L−1)	CYANO (μg L−1)
Baseline	0.56	0.07	63.54	32.84	17.58	10.17
No-till	0.57	0.07	63.70	32.98	17.56	9.96
Cover	0.55	0.07	62.84	32.43	17.55	9.76
Filter	0.55	0.07	63.15	32.65	17.54	9.80
Random combined	0.56	0.07	62.74	32.38	17.51	9.42
Source combined	0.54	0.07	62.28	32.17	17.46	8.93
Mouth combined	0.55	0.07	62.90	32.54	17.56	9.87

), 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₄, TN, NO₃, 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₄ 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₄, 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₄ spring loads were reduced by 40% (849 and 181 MTA in 2008, respectively), the maximum 30-day 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).

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₄ 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.

5. 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–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₄ 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.