Eutrophication management in a Great Lakes wetland: examination of the existence of alternative ecological states

. The degradation and loss of ecologically important wetlands has been a topical issue in the Great Lakes region, where 60 – 80% of the coastal wetlands have been lost since the 1800s. The present modeling study aims to guide the restoration efforts in Cootes Paradise marsh, one of the most degraded shallow wetlands in Southern Ontario. We use a process-based eutrophication model designed to reproduce the biotic competition among multiple phytoplankton and macrophyte functional groups. Our primary focus is to offer guidelines for wetland restoration by characterizing the ecophysiological processes of the autotrophic assemblage, such as the nutrient uptake from the water column and / or the sediment pore waters, the relative ability to harvest light and fuel photosynthesis, and temperature control of the algal / macrophyte growth and basal metabolism. We predict that the additional reduction of external phosphorus loading in Cootes Paradise could induce an abrupt, non-linear shift from the current turbid phyto-plankton-dominated state to a desirable clear macrophyte-dominated state. The emergence of this critical (or tipping) point, where the shift to another ecological state may occur, can be accelerated by the presence of a thriving macrophyte community with an enhanced ability to sequester phosphorus. However, it may also be delayed by the presence of a suite of biogeochemical mechanisms (often referred to as “ feedback loops ” ), such as the remobilization of legacy P due to sediment diagenetic processes, wind resuspension, bioturbation, hydraulic loading from local tributaries, water-level ﬂ uctuations, and the leachable P pool of dead plant material that can be returned into the water column through senescence. Our study identi ﬁ es the restoration actions required to minimize the likelihood of prolonged hysteresis and to facilitate a shift to a desirable ecological state in the foreseeable future. The areal expansion of aquatic vegetation will not only lead to the establishment of a thriving meadow and emergent vegetation community, but may also pave the way for submerged macrophytes through a suite of synergistic mechanisms. Additional point-source loading reductions will facilitate the transition to an alternative stable clear macrophyte-dominated state, but could also consolidate the future resilience of the marsh.


INTRODUCTION
In the Great Lakes region, a significant number of ecologically and economically important wetlands have been subjected to a variety of natural and anthropogenic disturbances, such as cultural eutrophication, land-use changes (e.g., urbanization), increased water levels, and bioturbation by the invasive common carp, Cyprinus carpio (Lougheed et al. 1998, Croft andChow-Fraser 2007). As a characteristic example, the Cootes Paradise marsh, located at the western end of Lake Ontario, has a long history of ecological degradation owing to excessive nutrient loading, algal blooms, and declines in the macrophyte community (Fig. 1). Cultural eutrophication, associated with point and non-point nutrient loading sources in the surrounding watershed, has been identified as the primary culprit for the impairment of this shallow ecosystem (Thomasen and Chow-Fraser 2012, Kim et al. 2016). The predominance of the invasive common carp also exacerbated the water turbidity through their nesting and feeding behavior that uproots submerged plants and perturbs bottom sediments (Thomasen and Chow-Fraser 2012). Consequently, the areal macrophyte coverage has declined from over 90% to less than 15% over the past ninety years, while the native macrophyte community has largely been extirpated by alien species (Chow-Fraser 2005). Since the mid-1990s, total phosphorus (TP) concentrations have varied anywhere from 202 AE 119 µg TP/L, at the innermost western side of the marsh, to 160 AE 81 µg TP/L, at its central area, and therefore the remediation efforts have a long way to go until the restoration target of 50 µg TP/L is realized (Theÿsmeÿer 2011). The same trends hold true for the prevailing chlorophyll α (42 AE 33 µg Chl a/L) and Secchi disk depth (<0.5 m) levels relative to the targeted values of 20 µg Chl a/L and > 1.5 m, respectively (Theÿsmeÿer 2011).
Restoration actions in Cootes Paradise marsh have been historically influenced by the popular paradigm of the existence of two alternative stable states: a clear, low-nutrient, macrophytedominated state and a turbid, high-nutrient, phytoplankton-dominated state (Scheffer et al. 2001). The notion that an abrupt, non-linear shift to a more desirable ecological state may be feasible was substantiated by the prevailing water quality conditions in 1997, when abnormally low spring temperatures caused a delay in fish migration into the marsh and allowed zooplankton to escape the predation pressure of planktivorous fish (Lougheed et al. 2004). This unexpected anomaly subsequently cascaded onto the system as a zooplankton-mediated control of phytoplankton biomass, water-transparency improvement, and proliferation of submerged vegetation in unvegetated shallow locations of the marsh (Lougheed et al. 2004). Bio-manipulation practices, such as the physical exclusion of large carp (>40 cm) from the system via the construction of a barrier (or fishway) at the outlet of Cootes Paradise, were intended to implement shock therapy measures that would facilitate the rapid restoration of aquatic vegetation (Lougheed et al. 2004). Consequently, the reduction in carp biomass from an estimated average marsh-wide biomass of 80-5 tons km −2 brought about a distinct waterclarity improvement, but the light availability near the sediment surface of the deep, open-water sites has not been adequate to consistently support the growth of submerged aquatic vegetation (Lougheed et al. 2004). Thus, notwithstanding the recent progress , the fact that the prevalence of favorable water quality conditions is typically short lived suggests that the remedial efforts to reduce external nutrient loading have not yet reached the critical levels that will allow the establishment of a resilient clear state in the marsh (Chow-Fraser 2005, Kim et al. 2016. Recent research has identified two major factors that may determine the success of the on-going restoration efforts in Cootes Paradise marsh and its likelihood to experience hysteresis patterns, whereby the degraded state will not be restored by an equal-sized reversal of the external nutrient loading that originally led to the system collapse (Scheffer et al. 2001, Schallenberg andSorrell 2009). First, the excessive legacy phosphorus in the sediments can exert control on the status of the wetland (Kelton and Chow-Fraser 2005). This is a critical factor that typically hinders restoration efforts in anthropogenically impacted freshwater systems (Orihel et al. 2017). In particular, Kelton and Chow-Fraser (2005) showed that internal loading due to resuspension, mineralization, and reflux may contribute up to 34% of TP loading to Cootes Paradise. This value was corroborated by Kim et al.'s (2016) net internal loading estimate of 36% for the growing season, determined via a mass-balance approach. Likewise, Mayer et al. (2005) estimated gross internal P loading rates ranging from 0.27 to 5.25 mgÁm −2 Ád −1 based on calculated diffusive fluxes at various sediment sites in Cootes Paradise marsh, which are approximately on par with Kim et al.'s (2016) estimated net contribution of 7.39 mgÁm −2 Ád −1 , considering the fact the former study did not correct for the temperature-dependence of sediment diagenetic processes. Some disparity between net and gross internal phosphorus loading estimates should also be expected due to differences in the time scales and methodology (Orihel et al. 2017). Water-level fluctuations are another critical factor that can modulate the response of Cootes Paradise to remedial measures, given that lower water levels and thus smaller water volumes imply lower dilution of point source and internal loads and correspondingly higher nutrient concentrations (Chow-Fraser 2005). Lower water levels may also ease the transmission of wind energy to the bottom sediments which in turn can magnify the release of phosphorus due to stirring and mixing (Prescott andTsanis 1997, Chow-Fraser 2005). Moreover, the disappearance of submerged macrophyte taxa in Cootes Paradise in the late 1990s was attributed to low water levels, which induced desiccation and increased reburial of propagules due to wind resuspension (Chow-Fraser 2005).
Considering the aforementioned factors that may impede the local restoration efforts, the goal of the present modeling study is to offer prescriptive guidelines for the management of Cootes Paradise marsh through three main objectives, aiming: (1) to characterize the ecophysiological processes responsible for resource (nutrient, light, temperature) procurement by the autotrophic (phytoplankton and macrophyte) assemblage in the water column and/or the sediment pore water; (2) to identify critical P loading levels below which  Kim et al. [2018] with the publisher's permission).
an abrupt, non-linear shift from the current turbid phytoplankton-dominated state to a clear macrophyte-dominated state can occur; (3) to examine the relative importance of a suite of biogeochemical mechanisms that may determine the resilience and future sustainability of the desirable ecological state. Our intent with these objectives is to inform the local research agenda and refine the ongoing monitoring efforts in Cootes Paradise marsh accordingly.

METHODS
The present study is based on the recently developed wetland eutrophication model or WEM (Kim et al. 2018) that characterizes the interplay among phytoplankton processes, fate and transport of phosphorus in the sediments, and macrophyte dynamics over the course of a 17-year period (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012). WEM's state variables are the dissolved-(DP w ) and particulate-(PP w ) phase phosphorus in the water column, three phytoplankton-(PHYT i , i = group A, B, and C) and macrophyte-(MAC j , j = meadow, emergent, and submerged) functional groups, sediment-particulate phosphorus (PP sd ) in labile (G1), refractory (G2), and inert (G3) forms, as well as dissolved phosphorus in sediment pore waters (DP sd ). Several biogeochemical processes are considered to determine the P fate and transport in the water column, including the uptake by the autotrophic assemblage, respiration/mortality of the simulated biotic components (phytoplankton and macrophytes), bacterialmediated mineralization, sediment resuspension, organic matter decomposition, particle settling, diffusive transport from the sediments, and advective mass exchanges with the Hamilton Harbour (Kim et al. 2018). DP sd is replenished through temperature-dependent sediment decomposition and is lost via macrophyte uptake and sediment diffusive reflux. PP sd is removed through sediment resuspension and permanent burial into deeper sediment layers. The model also contains a massbalance component of suspended solids (TSS w ) that accounts for the inflowing and outflowing particulate material, sedimentation, and the particles returned into the water column through wind resuspension.
WEM's phytoplankton submodel is based on Arhonditsis and Brett's (2005), and Ramin et al.'s (2011) models and simulates three phytoplankton functional groups, each with distinct ecophysiological characteristics. A diatom-like group (PHYT i = A ) represents an r-strategist with fast growth/metabolic rates, fast settling rates, and superior phosphorus kinetics; a cyanobacterialike group (PHYT i = C ) resembles to a K-strategist with slow settling rates, and inferior phosphorus kinetics; and a third group (PHYT i = B ) specified as an intermediate between diatoms and cyanobacteria. Similar to the phytoplankton governing equations, WEM simulates macrophyte growth as a function of nutrients, temperature, and light availability, but there are four key differences between macrophyte and phytoplankton processes with respect to resource (light and nutrient) procurement: (1) photosynthesis rates by phytoplankton and submerged macrophytes depend on the vertical light attenuation within the water column, whereas meadow and emergent plants do not experience light limitation (Kim et al. 2018); (2) meadows and emergent macrophytes exert light-shading effects on submerged macrophytes and phytoplankton (Kim et al. 2018); (3) macrophytes uptake phosphorus from both water column and interstitial waters with an uptake ratio that varies according to the relative abundance of the two pools, whereas phytoplankton rely exclusively on water-column nutrient availability (Granéli andSolander 1988, Christiansen et al. 2016); and (4) phosphorus dynamics within the phytoplankton cells account for luxury uptake (that is, phytoplankton uptake rates are modulated by both intracellular and extracellular nutrient concentrations, for example, Schladow andHamilton 1997, Arhonditsis et al. 2002), whereas phosphorus retention within the macrophyte tissues assumes that half of the uptaken phosphorus is rapidly released through respiration, and the other half is recycled through the slow decomposition of dead plant tissues (Asaeda et al. 2000). Details of the differential equations are presented in Kim et al. (2018), while the final model parameter values are presented in Appendix S1: Table S1. Our recent work presented a detailed sensitivity analysis (multiple-linear regression analysis, self-organizing maps) of WEM through which we identified critical parameters that shape the simulation of the interplay among phytoplankton, macrophytes, and nutrient loading (Kim et al. 2018). The present study goes one step further by examining the ability of WEM v www.esajournals.org to recreate the observed patterns in Cootes Paradise marsh, and subsequently to delineate the critical phosphorus loads and degrees of success of other remedial measures (e.g., restoration of native macrophyte community) that could trigger an abrupt shift to desirable water quality conditions.
The performance of WEM was evaluated at both the monthly and seasonal (May-October) scales. Model fit was examined based on four statistical measures: (1) coefficient of determination (r 2 ), (2) root mean squared error (RMSE), (3) relative error (RE), and (4) modeling efficiency (MEF).  (Mataya et al. 2017). Some of the predicted sediment reflux rates were validated against measured values from both the marsh itself and a neighboring eutrophic system, the Bay of Quinte, in northeastern Lake Ontario (Dittrich et al. 2015). We examined the competition mechanisms between submerged macrophytes and phytoplankton biomass under an incremental reduction of P loading until 50% of the current levels (27.12 kg/d). The latter exercise was combined with a local sensitivity analysis (AE 30%) to selected model parameters: diffusivity, decomposition rate of particulate phosphorus in the sediments (PP Dcp(k) : conversion of PP sd to DP sd , see Appendix S1: Table S1), submerged macrophyte phosphorus P uptake fraction between sediment and water column, and light attenuation coefficients for both submerged macrophytes and phytoplankton. To account for the sensitivity of our projections to macrophyte P sequestration (Kim et al. 2018), we also assessed how submerged macrophytes would influence the prevailing water quality conditions to varying levels of exogenous loading reduction (up to 70%) and submerged macrophyte physiological attributes (respiration and mortality rates, and phosphorus fraction in biomass).

Model performance
Based on the detailed literature review presented by Kim et al. (2018; see their Table 2) and plausible ranges derived from previous modeling work in Hamilton Harbour and Cootes Paradise marsh (Ramin et al. 2012, Kim et al. 2013, we created 50,000 parameter vectors and nutrient loading levels (model runs) with Latin Hypercube Sampling that were then used to calculate the average values for each model end point, as well as the 95% uncertainty bounds. Prior to this exercise, we determined the behavioral space of the model that delineated the lower and upper acceptable values for each state variable. This acceptable predictive envelope was used to scrutinize our model runs and exclude the ones that represented ambient levels/conditions never registered in Cootes Paradise marsh during our study period. WEM displayed variant degree of fidelity in reproducing the average monthly TP, DP w , Chl a, and TSS w concentrations v www.esajournals.org in Cootes Paradise ( Fig. 2 and Appendix S1: Fig. S1). The vector of calibration parameter values that provided the best fit is presented in Appendix S1: Table S1.
Among the goodness-of-fit statistics calculated for the monthly scale (Table 1), we found satisfactory RE values for TP, Chl a, and TSS w concentrations, ranging from 55% to 65%, but an exceptionally high RE for DP w (136%). MEF values were also suggestive of model performance that was (at best) as good as the observed average of the corresponding water quality variables. Likewise, the high RMSE values (TP 83.2 µg/L, DP w 31.6 µg/L, Chl a 28.9 µg/L, and TSS w 42.6 mg/L) largely reflect the model's inability to capture ambient levels that deviated significantly from the typically prevailing water quality conditions. Namely, the model was unable to capture the unusually low TP concentrations at Cootes Paradise marsh in 1997, resulting in an overestimation Fig. 2. Predicted versus observed monthly TP, DP w , and Chl a concentrations in Cootes Paradise Marsh. Solid lines indicate the predicted average values, while dotted lines delineate the 95% uncertainty bounds, as determined by the exogenous nutrient loading variability as well as the model sensitivity to key ecophysiological parameters used during our calibration exercise. by approximately 50 µg TP/L on a monthly basis during the growing season. In contrast, the model significantly underestimated the substantially higher than average TP and Chl a levels on a monthly basis observed in 2001 and 2007 ( Fig. 2; see also following discussion). When the model's predictive capacity was further examined at a coarser (seasonal) resolution, RE (TP 23.7%, DP w 62.8%, Chl a 30.6%, TSS w 60.1%) and RMSE (TP 54.8 µg/L, DP w 12.0 µg/L, Chl a 15.9 µg/L, and TSS w 33.8 mg/L) improved significantly, but MEF values (TP − 0.51, DP w − 4.32, Chl a 0.01, TSS w − 3.35) were still mostly negative ( Table 1). The model reproduced year-to-year variability of the seasonal TP average values, including the last five years of the study period (2008)(2009)(2010)(2011)(2012), when TP was lower than the long-term seasonal average (Appendix S1: Fig. S2). However, similar to the finer-monthly-resolution, the predicted seasonal TP differed significantly from the observed values in the marsh in 1997 and 2007. Interestingly, if we consider the 95% uncertainty bounds, as modulated by the exogenous nutrient loading variability as well as the model sensitivity to key ecophysiological parameters used during our calibration exercise, the model satisfactorily captured the yearto-year variability of the seasonal DP w and Chl a average levels in the marsh (Appendix S1: Fig. S2). Peaks in TSS w concentrations (Appendix S1: Fig. S1) largely coincided with elevated TP concentrations in the water column (Fig. 2), while the highest total phytoplankton (and functional group) biomass levels coincided with DP w minima during the summer period ( Fig. 2 and Appendix S1: Fig. S3). Sediment dissolved phosphorus (DP sd ) and macrophyte biomass (MAC j = EMG,MDW ) levels exhibited a sinusoidal pattern and generally peaked in the summer (Appendix S1: Figs. S1 and S4). Notably, our simulations of sediment phosphorus concentrations (DP sd and PP sd ) were indicative of a gradual decline over the last 15 years in the marsh (Appendix S1: Fig. S1), while the submerged (MAC SUB ) macrophyte community appears to have reached a point of collapse after 1999-2000 (Appendix S1: Fig. S4).

Characterization of the P cycle in Cootes Paradise marsh
Daily phosphorus fluxes associated with each ecophysiological process throughout the growing (May-October) season averaged across the entire study period (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) are illustrated in Fig. 3. External TP and DP w loads contributed approximately 14.91 and 12.21 kg/d and were distinctly greater than internal loads (i.e., TP fluxes from sediment and macrophytes). In particular, internal loads stemming from sediment fluxes by diffusion and resuspension collectively contributed 12.0 kg/ d to the water column. Our model suggested that approximately 6.23 kg/d is lost from the water column through sedimentation and 22.30 kg/d through outflows to Hamilton Harbor. Moreover, an additional 6.90 kg/d is lost from the marsh as P that is sequestered in phytoplankton cells, and flows out into Hamilton Harbor. The permanent loss of phosphorus from the system through burial to the deeper sediments accounts for 2.25 kg/d, resulting in a net phosphorus settling rate related to the sediments (i.e., sedimentation + macrophyte mortality − sediment reflux/resuspension − macrophyte uptake − burial) of −9.18 kg/d. Based on their low abundance levels in the system during the study period, macrophytes appear to play a small role in the phosphorus cycle, removing a net of 1.35 kg/d from the sediments and contributing 0.32 kg/d into the water column via respiration. The bacteria-mediated mineralization of particulate material further supplies 0.75 kg/d, while phytoplankton mortality/basal metabolism is responsible for an additional 1.82 kg/d. By the same token, our model also predicts that the phytoplankton assemblage uptakes 13.22 kg/d from the water column during the growing season.

TP-Chl a relationship versus nutrient/hydraulic loading variability
One of the critical features of a eutrophication model is its capacity to recreate the fundamental relationships between key water quality variables and external drivers of the system, such as nutrient and hydraulic loading variability (Arhonditsis and Brett 2004). In this regard, we compared the TP-Chl a relationship between observed and predicted seasonal (May-October) averages during the growing period (May-October) from 1996 and 2012 ( Fig. 4 and Appendix S1: Fig. S5). WEM does capture the linear TP-Chl a relationship in the marsh with similar slope coefficients, that is, 0.35 versus 0.39 μg Chl aÁL −1 (μg TPÁL −1 ) −1 . The difference between the intercepts of the two equations stems from the Chl a-TP pair associated with the aforementioned outlying year 1997, when unusually v www.esajournals.org favorable water quality conditions prevailed in the system. Despite their strong causal linkage, TP and Chl a concentrations are also modulated by the riverine nutrient and hydraulic loading. Considering that the three local creeks are the largest TP sources in Cootes Paradise Marsh, including the urban-to-stream runoff (see Fig. 1 in Kim et al. 2016), we regressed the observed and predicted TP and Chl a concentrations against the riverine flow-weighted average TP levels and the hydrological flushing rates (i.e., number of water-renewal cycles per season; Fig. 5). Both observed and predicted seasonal TP concentrations displayed a positive relationship with tributary TP concentrations and a negative one with the flushing rates (Figs. 5a-d). Notably, the empirical evidence is suggestive of a stronger signature of riverine nutrients on ambient TP levels relative to the one depicted by our model (slope coefficients 2.22 vs. 0.87 μg TP/L in the water column per μg TP/L inflowing in the marsh). In contrast, the flushing rates appear to influence more strongly the predicted than the observed TP concentrations (−1.32 vs. −1.14 μg TP/L in the water column per renewal cycle during the growing season). Observed Chl a response to nutrient/hydraulic loading were drastically different from the predicted responses (Fig. 5e-h). In particular, the model plausibly predicts a positive linkage with flow-weighted TP concentrations (Fig. 5f) and a negative one with the flushing rates (Fig. 5h). Considering the wide array of confounding factors that can shape these relationships (some of which were inevitably unaccounted for by the model), observed phytoplankton biomass bore almost no relationship with the inflowing P masses into the marsh (Fig. 5e) and a counter-intuitive positive association with the flushing rates (Fig. 5g).

Phytoplankton-macrophyte competition patterns under different nutrient loading regimes
We examined the competition mechanisms and likelihood of coexistence patterns between submerged macrophytes and phytoplankton biomass under different nutrient loading reduction scenarios. The same exercise was also repeated by considering the role of additional ecophysiological/biogeochemical processes that could potentially modulate the nature and strength of the phytoplankton-macrophyte relationship, such as the diffusive nutrient exchange between the water column and interstitial waters, decomposition rate of organic matter in the sediments, the fraction of P uptaken by submerged macrophytes from the water column pool relative to that in the sediment porewater, and the ability of submerged macrophytes and phytoplankton to shape the vertical light attenuation and consequently the broader illumination of the water column. Chl a concentrations and submerged-macrophyte density displayed an inverse linear relationship in response to nutrient-loading reduction, shifting from a eutrophic state, where a seasonal phytoplankton biomass of nearly 45 μg Chl a/L coexisted with submerged-macrophyte density of 20 g DW/m 2 , to a more desirable water quality state, where the standing phytoplankton biomass is reduced to 20 μg Chl a/L and submerged macrophytes proliferate in the system, that is, 120 g DW/m 2 (Fig. 6a). The strength of this inverse relationship is influenced by the suite of confounding factors (associated with respiration and mortality rates, and phosphorus fraction in biomass) examined to varying degrees. Interestingly, the strength of diffusive exchange and the sediment organic matter decomposition rates did not seem to exert a major control on the resource competition within the autotrophic assemblage and the associated coexistence patterns (Fig. 6b, c). In stark contrast, our model predictions appeared to be particularly sensitive to the specification of factors that determine the relative ability of macrophytes to compete for light and nutrients relative to phytoplankton. Specifically, Chl a concentrations are moderately decreased and macrophytes remained at very low density levels with decreasing nutrient loading, when macrophytes were specified to have greater reliance on bioavailable phosphorus in the water column than on the interstitial water pool (Fig. 6d). Light attenuation was another critical factor to shape the submerged macrophyte-phytoplankton relationship (Fig. 6e, f), whereby the shading effect of phytoplanktonand less so that of submerged macrophytes-on the water column not only determine the composition of the autotrophic assemblage, but also the steepness of the trajectory followed during the nutrient loading reduction. Most notably, assigning a low light-attenuation coefficient to phytoplankton resulted in a moderate decrease of Chl a concentrations and very low macrophyte density until a 30-35 μg Chl a/L level was reached (Fig. 6f). Below that threshold level, phytoplankton followed a steeper declining trajectory that was accompanied by a commensurate increase in macrophyte abundance for each incremental nutrient loading reduction. v www.esajournals.org A final modeling experiment examined the response of phytoplankton to TP loading reduction in the presence of a submerged macrophyte assemblage with lower metabolic loss (respiration and mortality) rates and a greater amount of P sequestered per unit of biomass produced (Fig. 7).
We found that a nutrient loading reduction from anywhere between 10-20% (≈23.1 kg/d) and 60-70% (≈9.5 kg/d) could lead to an approximate decrease from 42 to 25 µg Chl a/L for the default macrophyte specification (ωR SUB = 0.015 d −1 , ω DSUB = 0.0480 d −1 , and μ P/DW = 0.0025 g Pg/DW; v www.esajournals.org 10 February 2021 v Volume 12(2) v Article e03339 Fig. 6. Analysis of the submerged macrophyte-phytoplankton relationship in response to (a) TP loading reduction by 5% increments from the present conditions (dots from left to right). The rest of the panels present the same exercise coupled with variations of the (b) diffusion rate (ω Dif ), (c) sediment decomposition rate (ω Dcp(k) ), (d) submerged-macrophyte P-uptake ratio from the water column (1-Ψ Pup(j) ), (e) submerged-macrophyte lightextinction coefficient (kext SUB ), and (f) phytoplankton light-extinction coefficient (kext Chla(i) ). See Appendix S1: Table S1 for the complete list of model parameters and to Kim et al. (2018) for the WEM differential equations. v www.esajournals.org 11 February 2021 v Volume 12(2) v Article e03339 see also Appendix S1: Table S1) that recreates a scenario of low macrophyte density with weak P sequestration and high metabolic losses. In contrast, the same nutrient loading reduction strategies, under a scenario of a thriving macrophyte community with strong P sequestration and low metabolic losses, led to Chl a concentrations consistently lower than 20 µg/L. Thus, the two distinct macrophyte specifications appear to render a fairly resilient system with respect to the ability of nutrient loading reductions to induce dramatic changes in the prevailing water quality conditions. Consistent with the findings by Kim et al. (2018); see their Fig. 4, we also identified a critical submerged-macrophyte characterization between the two extreme states, that is, a threshold biomass density and associated impact on P cycle, where the benefits from the management actions of loading reductions in the watershed can be capitalized for greater water quality improvements (Fig. 7).

DISCUSSION
Degraded freshwater ecosystems often experience a time lag before showing any response to restoration actions, with shallow wetland ecosystems, like Cootes Paradise marsh, requiring more time and effort to recover than deeper lake and riverine ecosystems (Zedler 2000, Thomasen andChow-Fraser 2012). The degree of success for wetland management implementation depends on a suite of biogeochemical mechanisms that can prolong the prevalence of undesirable conditions, such as remobilization of legacy P due to sediment diagenetic processes, wind resuspension, bioturbation, hydraulic loading from the local tributaries and/or water-level fluctuations, the amount of dead plant material comprising a leachable P pool that can be returned into the water column through senescence or bacterial decomposition (Moss et al. 1996, Ibelings et al. 2007). The present modeling study predicts that Fig. 7. Predictions of seasonal average Chl a concentrations in response to TP loading reductions at various degrees of submerged aquatic vegetation (SAV) P sequestration and metabolic loss rates. The latter processes were emulated by a 20% reduction of macrophyte respiration (ωR SUB = 0.015 → 0.012 d −1 ) and mortality (ω DSUB = 0.0480 → 0.0384 d −1 ), as well as a fivefold increase of P content in macrophyte tissues (μ P/DW 0.0025 → 0.0125 g P g DW −1 ). Circles and diamonds represent TP loading reduction scenarios from 10% to 70%. Each dot represents the average for each combination of P loading level and P sequestration/metabolic loss rates, calculated over all the growing seasons during the 1996-2013 period. a drastic reduction of the external nutrient loading in Cootes Paradise could induce an abrupt, nonlinear shift from the current turbid-phytoplankton dominated state to its former clear-macrophyte dominated state. However, the timing of the emergence of this tipping point largely depends on the presence of a thriving macrophyte community with an enhanced ability to sequester phosphorus. Thus, the resurgence and composition of the macrophyte assemblage, the characterization of their ecophysiological processes pertaining to natural-resource procurement, and their competition patterns with phytoplankton represent focal points to maximize the benefits from the local remedial measures.

Recovery and establishment of the optimal macrophyte community structure
In alignment with the objective to recover the macrophyte areal coverage and community structure, one of the overarching themes of the on-going management actions in Cootes Paradise marsh revolves around the plantings of emergent plants, mainly cattails (Typha sp.), and less so meadow and submerged macrophytes . Owing to the sustained highwater levels in Lake Ontario over the past 30 yr, a significant portion (≈25%) of the shoreline in Cootes Paradise marsh still remains without emergent vegetation and is thus subjected to erosion, which is further exacerbated by the wind fetch and waves. For emergent seedling germination and subsequent shoreline stabilization to occur, a maximum summer water level of less than 74.75 mean sea level is required (Chow-Fraser 2005), but this condition has been rarely met over the past two decades. and therefore, all emergent plant resurgence has been achieved through active plantings and subsequent expansion by rhizomatic growth. Meadow marsh is a priority habitat that is often used as an environmental indicator for Lake Ontario water-level regulation. The potential meadow zone in Cootes Paradise marsh is already vegetated, but the plant community is dominated by two highly aggressive non-native plant species; common reed (Phragmites australis) and Eurasian manna grass (Glyceria maxima; Theÿsmeÿer et al. 2016). Even though the presence of the former species has been effectively managed, one of the on-going challenges is to protect the intact native meadow marsh habitats in favor of the keystone plant, lake-bank sedge (Carex lacustris). Importantly, both empirical evidence and theoretical evidence suggest that if the restoration targets pertaining to emergent and meadow macrophytes are realized in Cootes Paradise, it will not only be beneficial in supporting native insects and wildlife, but it could also pave the way for submerged macrophytes to sustain themselves and effectively proliferate in the system through a series of synergistic mechanisms, such as stabilization of the sediments and increased competition with the algal assemblage , Kim et al. 2018.
In stark contrast, the submerged vegetation appears to have a competitive handicap against phytoplankton, when the areal coverage and density of emergent and meadow macrophytes decline (Kim et al. 2018;see also Fig. 6). Notwithstanding the capacity of submerged vegetation to uptake nutrients by/through two pathways (interstitial waters and water column), the competitive advantage of phytoplankton stems from the ability of the current algal assemblage to adapt in poorly illuminated environments that still prevail in Cootes Paradise marsh . In particular, the continuous presence of cryptomonads (Cryptomonas, Rhodomonas) and euglenophytes (Euglena, Lepocinclis, Phacus) in the summer phytoplankton community, known to display elevated acclimation to turbid and low-light conditions, lends support to the latter conjecture about the importance of light in shaping the phytoplankton-submerged vegetation relationship and the need to capitalize upon the potential synergistic effects among all the major functional groups of the macrophyte (emergent, meadow, submerged) assemblage in order to consolidate the shift toward an improved water quality state (Kim et al. 2018). In the same context, our analysis showed that when favorable competition traits for light are assigned to phytoplankton (e.g., low self-shading effects combined with elevated shading on submerged macrophytes), the submerged vegetation does not respond to the reduced nutrient loading regime until a 30-35 μg Chl a/L level is reached (Fig. 6f). Above that threshold level, phytoplankton varied with the incremental nutrient loading reduction, but macrophyte density remained at very low levels (≈20 g/m 2 ). This predicted control of v www.esajournals.org 13 February 2021 v Volume 12(2) v Article e03339 macrophyte abundance and light penetration driven by phytoplankton could indirectly reflect the shading effect of epiphytic algae, as existing evidence suggests that they may be the actual causal factor and the elevated phytoplankton growth is subsequent (Granéli andSolander 1988, Sand-Jensen andBorum 1991). It is also worth mentioning that WEM's characterization of the ecophysiological processes draws upon literaturebased estimates for submerged macrophyte species typically present in Cootes Paradise marsh, including the dominant ones Potamogeton pectinatus, P. foliosus, P. crispus, and several other less common species, such as Myriophyllum spicatum, Zannichellia palustris, Elodea canadensis, and Ceratophyllum demersum Stefan 2006, Lougheed et al. 2004, see also Table 2 in Kim et al. 2018).

Linking P dynamics with the presence of macrophytes in the marsh
Our analysis identified macrophyte P sequestration (amount of phosphorus uptake per unit of biomass produced) as an important regulatory factor that can influence the TP budget within the wetland. Consistent with empirical evidence of a predominantly acropetal translocation movement (Granéli and Solander 1988), our model postulates that the phosphorus uptake of submerged macrophytes takes place primarily from the sediment pore waters via the roots and less so from the water via the leaves. Importantly, the latter assumption was relaxed to accommodate the idea that the two dominant Potamogeton species (P. pectinatus, P. crispus) in Cootes Paradise have different vascular systems. Consequently, they may be able to display a higher degree of uptake through their shoot system and basipetal movement (P. pectinatus; Welsh and Denny 1979). Under this scenario, it was shown that greater reliance of submerged macrophytes on nutrients from the water column may not necessarily be advantageous (Fig. 6d). Namely, in the context of WEM parameterization, our simulated submerged plants can more easily reach critical biomass density levels, when their P requirements are mainly met from the substrate instead of the water. In doing so, macrophytes can be viable competitors against phytoplankton, which closely rely upon the time-varying P in the water column. It should again be noted that our study does not explicitly consider the fact that macrophytes are covered by an active epiphyte community with a higher nutrient uptake per unit area rate, which is another critical factor that can shape the competition mechanisms within the autotrophic assemblage (Howard-Williams and Allanson 1981).
Another critical facet of the macrophyte physiology with important implications for the ecosystem functioning is related to P losses from respiration and mortality (Kim et al. 2018, see also Fig. 7). Generally, phosphorus retention within plant tissues is modulated by different physiological strategies for nutrient storage and senescence. P translocated to below-ground biomass (roots and rhizomes) represents a longterm storage pool, whereas P contained in aboveground plant tissues (e.g., shoots) can be released into the water column within a shorter timeframe. Notwithstanding the conflicting empirical evidence on foliar release of phosphorus by actively growing and healthy macrophytes (Carignan andKalff 1982, Gabrielson et al. 1984), the general consensus in the literature is that their contribution is relatively minor compared with the quantities of nutrients released from the leachable P pool of dead plant material and returned into the water column through senescence or bacterial decomposition (Granéli andSolander 1988, Herb andStefan 2003). With regard to the contribution of decaying plant tissues, there appears to be a rapid initial loss of phosphorus due to leaching, which is independent of temperature and fragmentation of the plant material (Carpenter 1981). On the other hand, the slower release of refractory organic phosphorus is mediated by microbial mineralization, and as such is controlled by temperature, nutrient supply, availability of energetically favorable electron acceptors, and nature of the plant tissues (Granéli and Solander 1988). Emergent macrophytes typically possess large, perennial storage organs for carbohydrates and supporting tissues which are resistant to microbial decomposition, while submerged plants have mainly fine roots and have relatively low cellulose content and can thus be more easily subjected to mineralization upon death (Twilley et al. 1986, Leisti et al. 2016. The latter pattern is highly relevant to the recent (2014-2016) conditions in Cootes Paradise, when dominant submerged macrophyte species of the Pondweed family (e.g., P. foliosus) rapidly proliferated in the shallow, slow-moving waters of the marsh at the beginning of the growing season, when there were remarkably favorable water quality conditions (TP = 57 AE 38 μg/L, Secchi disk depth = 63 AE 18 cm), but soon thereafter they collapsed and the decomposition of their dead tissues brought about a significant water quality deterioration toward the end of the growing season, for example, TP = 143 AE 48 μg/L and Secchi disk depth = 40 AE 15 cm (Yang et al. 2020).

Role of internal nutrient loading, water-level fluctuations, and other uncertainties
Our internal loading estimate of 12.0 kg/d, emanating from sediment reflux (4.9 kg/d) and wind resuspension rates (7.1 kg/d) is somewhat lower than the value (16.2 kg/d) reported by Kim et al. (2016). Likewise, Kelton and Chow-Fraser (2005) projected an internal loading rate of 13.0 kg/d, which puts greater weight on net mineralization and/or release rate (8.8 kg P/d) and less so on wind resuspension (4.2 kg P/d) relative to our predictions. In the same context, Mayer et al. (2005) estimated phosphorus diffusive exchange ranging from 0.27 to 5.25 mgÁm −2 Ád −1 at various sediment sites in Cootes Paradise marsh, which is comparable with our temperature-corrected average estimate of 2.47 mgÁm −2 Ád −1 during the growing season. Considering the significant spatial variability of P reflux in the marsh (Kelton andChow-Fraser 2005, Mayer et al. 2005), our internal seasonal diffusive fluxes as well as the pore-water soluble reactive phosphorus concentration gradients (Fig. 8) are remarkably similar to the values registered in three sites in the Bay of Quite, Belleville 3.6 AE 2.5 mgÁm −2 Ád −1 , Napanee 2.7 AE 0.8 mgÁm −2 Ád −1 , and Hay Bay 2.3 AE 0.5 mgÁm −2 Ád −1 (see also Fig. 4 in Markovic et al. 2019); a neighboring embayment in northeastern Lake Ontario with long history of eutrophication problems. In the shallow sites of the latter system, Markovic et al. (2019) predicted long-term diagenetic P release that was half of the magnitude of the aforementioned short-term diffusive release. The long-term release is the combined effect of multiple processes of sediment P mobilization or diagenetic transfers from labile to more stable solid forms and thus represents the base potential for sediment P mobilization or (simply put) the long-term trends of P accumulation and burial. The significant discrepancy between short-and long-term P release is characteristic of shallow aquatic environments, like Cootes Paradise, where short-term summer pulses of increased mobilization of labile P forms (freshly deposited algal cells, redox-sensitive P extracted with bicarbonate dithionite, or other organic P forms, e.g., polyphosphate) could be triggered by intermittent changes in environmental conditions and phytoplankton dynamics (Roelofs et al. 2002, Tammeorg et al. 2015, Parsons et al. 2017, Markovic et al. 2019. Interestingly, Parsons et al. (2017) showed that short-term redox fluctuations resulted mostly in internal redistribution of redox labile P to alternative binding sites in the solid phase, rather than release to the aqueous phase for sediments collected from West Pond. Notwithstanding the uncertainty regarding the short-and long-term contributions of the nutrient-rich, anoxic sediments in Cootes Paradise marsh, our present analysis reinforces earlier predictions that internal loading is responsible for approximately 30-35% of the TP loading.
Along with their ability to shape the resurgence and structure of the macrophyte assemblage, the water-level fluctuations can also directly influence the water quality conditions in Cootes Paradise. Lower water levels or smaller water volumes lead to less dilution and higher nutrient concentrations. Wind energy is also transmitted more rapidly to the bottom when water levels are low, which in turn can magnify the stirring of the sediments and subsequently the replenishment of the water column with phosphorus from the sediment pool (Prescott andTsanis 1997, Chow-Fraser 2005). Likewise, warmer conditions can accelerate the loss of water through evapotranspiration and tighten the water column-sediment coupling in shallow waters (Jeppesen et al. 1997, Søndergaard et al. 2003, Drexler et al. 2004, Panin et al. 2006. We thus examined how ambient water quality would respond to distinct changes on the hydrological regime in the area, as approximated by the year-to-year variability of the water volume, flushing rates, and volume of outflowing water masses (Fig. 9). Our model predicts that the prevalence of warmer and drier conditions can be responsible for up to twofold increase in Chl a and TP concentrations and can undermine the on-going restoration efforts; especially if we consider that the water level can drop by nearly a meter between the beginning and end of the growing season. Similar negative relationships were recreated between water levels and TSS w concentrations or emergent macrophyte density, which were qualitatively on par with the empirical equations presented by Chow-Fraser (2005); see her Figs. 2a, 4. Consistent with these predictions, Smith et al. (2001) asserted that the recovery of emergent vegetation in 1999 could have been attributed to the lower water levels registered during that year, while Harwell and Havens (2003) proposed two threshold water levels for evaluating the resilience of submerged macrophytes; a maximum threshold, above which light could become limiting; and a minimum threshold, below which the prevailing conditions will be excessively dry. In fact, Chow-Fraser (2005) hypothesized that the elimination of submerged macrophytes from Cootes Paradise in 1999 was caused by the effects of desiccation due to low water levels, as well as the re-burial of propagules due to increased wind resuspension. The notion that the boundary between perennial emergent and submergent vegetation is a function of the water cycle is not new in Cootes Paradise marsh . What the present study underscores though is that the direct effects of the water cycle on water quality could pose another challenge in our restoration efforts of the system, and thus the active planting efforts for the establishment of a balanced macrophyte community should continue to be closely connected with the on-going planning efforts among the local stakeholders to further mitigate the impact of point-and nonpoint pollution sources.
Founded upon the characterization of the interplay among phytoplankton processes, fate, and transport of phosphorus in the sediments, and macrophyte dynamics, WEM generally reproduced the observed water quality conditions in Cootes Paradise marsh in a satisfactory manner. Nevertheless, there were instances in which the food-web configuration of the model was inadequate to capture the observed dynamics. The first example was the anomalous year 1997, when carp exclusion and low spring temperatures delayed the migration of spawning fish, released zooplankton from fish predation and ultimately led to unusually low ambient TP levels (Lougheed et al. 2004). A second major model mismatch was found in 2007 and could again be attributed to external factors, unaccounted for by the model. In particular, the significant underestimation of the high end of summer-early fall TP levels corresponds to a major pesticide spill that resulted in an excessively high diazinon concentration (≈94 µg/L) and led to massive killings of aquatic invertebrates/small fish shortly after (Theÿsmeÿer and Galbraith 2007). As 5-15 million fish died on the bottom of the marsh, the decomposition of their v www.esajournals.org tissues profoundly elevated nutrients, which the current model structure cannot reproduce without a major re-configuration (Theÿsmeÿer and Galbraith 2007). In designing the next iteration of WEM, we believe that its learning capacity will improve by introducing the dynamic interactions of the autotrophic assemblage with zooplankton (Thomasen and Chow-Fraser 2012). Another mechanistic augmentation will be to increase the granularity of the study of macrophyte growth and decay processes by focusing on different plant parts (shoots, roots, and tubers), but only if the empirical knowledge from the system is commensurate to the added mathematical complexity (Asaeda et al. 2000). However, above all else, we stress the need to improve the existing estimates of critical forcing functions, such as the non-point nutrient loading (e.g., estimates of daily flow discharge for Chedoke and Borer's creeks are still based on areal weighting factors and TP concentrations mainly from baseflow conditions) and water-level fluctuations (Kim et al. Fig. 9. TP and Chl a concentrations according to different hydraulic regimes between wet (closed circle) and dry conditions (open circle). The regimes are modulated by (a-b) flushing rate, expressed as number of times that the water in the marsh is renewed during the growing season; (c-d) outflow; and (e-f) water volume. The dry scenario was defined by an 1°C increase in water temperature, 50% decrease in water volume, and 70% decrease in outflows relative to the wet scenario.
v www.esajournals.org 2016). Regarding the latter input, the data uncertainty pertaining to individual fluxes in the water cycle (i.e., evapotranspiration, precipitation, groundwater, runoff, outflows) posed challenges in balancing the water budget. In particular, we note that the existing outflow estimates do not accurately capture the diurnal and semi-diurnal flow reversals between Cootes Paradise and Hamilton Harbor and therefore had to be adjusted by year-specific fudge factors (Skafel 2000).

Lessons learned, future perspectives
To recap, consistent with an extensive body of the limnological literature (Scheffer et al. 2001, Scheffer and Carpenter 2003, Folke et al. 2004), our modeling study predicts that the transition from the turbid phytoplankton-dominated state to clear macrophyte-dominated state is possible and can be accelerated by the combined effects of exogenous loading reductions and restoration of the macrophyte community. Importantly, the tipping point of the regime shift varies with the degree of proliferation of submerged macrophytes in the marsh through their ability to sequester P from the water column and recycle nutrients through respiration and/or senescence of their plant tissues. Our study renders support to the rationale and priorities of current restoration efforts to maintain system connectivity and restore the underlying conditions for biodiversity recovery and sustainability. The ongoing areal expansion of aquatic vegetation (60 ha in the 1990s, 131 ha in 2015, and 270 ha set as the latest objective) will not only lead to the establishment of a thriving meadow and emergent vegetation community in the system, but could also pave the way through a suite of synergistic mechanisms for submerged macrophytes to effectively proliferate and have a positive impact on the ecological status of Cootes Paradise marsh ). This is one of the key lessons learned from our analysis that could be used for wetland restoration practices elsewhere.
Our modeling analysis also contributes to the on-going discussion regarding the implementation of additional loading reduction from Dundas WWTP, which are predicted to facilitate both the transition to an alternative stable clear macrophyte-dominated state and its future resilience. While the uncertainties pertaining to water-level fluctuations, climate change, and episodic pollution events that stem from intense anthropogenic activities in the surrounding watershed, could pose challenges for the integrity of the ecosystem in the future, we believe that the positive trajectory of the restoration efforts over the past 20 yr offers plenty of leeway to continue the efforts to reduce point-and non-point sources loading. After all, the considerable socioeconomic direct and indirect benefits gained from ecosystem services of Cootes Paradise marsh make its restoration an absolutely worthwhile investment (Yang et al. 2020). Except from the uncertainties stemming from the data availability, there are at least two facets of the broader role of macrophytes that could be potentially incorporated in the next iteration of our model; namely, the allelopathy of submerged macrophytes to phytoplankton growth (Kim et al. 2018), and the refuge effect of submerged macrophytes to zooplankton. However, given that the potential role of allelopathy has not been studied in Cootes Paradise marsh and the insufficient available data for zooplankton, both processes have not been considered. Considering the central role of macrophyte restoration for the water quality of the marsh, another critical piece of missing information will be the characterization of certain critical ecophysiological processes pertaining to P macrophyte uptake from the sediments and water column with measurements from the system. After all, one of the major benefits of any modeling exercise is to pinpoint knowledge gaps (understudied processes, missing data for critical model end points) and shape future monitoring and research accordingly.