Spatial identification of critical nutrient loads of large shallow lakes Implications for Lake Taihu (China)

of critical nutrient loads for different locations throughout the lake. Our analysis shows that the largest part of Lake Taihu follows a nonlinear load-response curve without hysteresis. The corresponding critical nutrient loads vary within the lake and depend on management goals, i.e. the maximum allowable chlorophyll concentration. According to our model, total nutrient loads need to be more than halved to reach chlorophyll-a concentrations of 30 e 40 m g L (cid:1) 1 in most sections of the lake. To prevent phytoplankton blooms with 20 m g L (cid:1) 1 chlorophyll-a throughout Lake Taihu, both phosphorus and nitrogen loads need a nearly 90% reduction. We conclude that our approach is of great value to determine critical nutrient loads of lake ecosystems such as Taihu and likely of spatially heterogeneous ecosystems in general. © 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license


Introduction
Extensive eutrophication threatens numerous lakes worldwide, leading to major challenges for drinking water supply, food security and public health (Brooks et al., 2016). To minimize negative consequences, water quality managers look for successful and efficient measures. In this process, they need to carefully walk on the thin line between the desired water quality and the available financial resources (Groffman et al., 2006). Importantly, success of management actions depends on the identification of ecological thresholds (Kelly et al., 2015). These thresholds set limits to sustainable use of resources within the safe operating space (Rockstr€ om et al., 2009) and exceedance leads to undesired ecological consequences such as toxic phytoplankton blooms (Groffman et al., 2006). Currently, there is a great interest in ecological thresholds (Kelly et al., 2015). In the first place, ecological thresholds support adaptive management by defining limits to maximum stress, and, secondly, ecological thresholds support recovery by setting a restoration goal (Kelly et al., 2015). A common threshold used in limnology is the concept of critical nutrient loads.
We define critical nutrient load as a maximum nutrient load an ecosystems can absorb, while remaining in a good ecological state. This definition of the critical nutrient load includes both gradual and sudden changes in water quality.
To estimate the critical nutrient load of a specific lake, it is important to know the type of load-response curve. The loadresponse curve of lakes is generally defined as the response of phytoplankton (i.e. chlorophyll-a) to a range of nutrient loads. Literature distinguishes three types of load-response curves which differ in linearity and the presence of hysteresis ( Fig. 1) (e.g. Scheffer et al., 2001). The first type is a linear load-response curve for which a reduction in nutrient load results in a proportional decrease in chlorophyll-a, irrespective of the present nutrient load (Fig. 1a). The critical nutrient load for lakes with a linear loadresponse curve is quite sensitive to the human-defined maximum allowable chlorophyll-a concentration. These human-defined maxima closely relate to the lake's societal functions (e.g. recreation or drinking water). Depending on the maximum allowable chlorophyll-a concentration, the critical nutrient load varies between a relatively high ( Fig. 1a blue line (1)) and low value (Fig. 1a red line (2)). The second type of load-response curve is nonlinear (Fig. 1b). The critical nutrient load for lakes with a nonlinear loadresponse curve depends on the positioning and slope of the fold. Therefore, the critical nutrient load only slightly depends on human-defined limits to the maximum allowable chlorophyll-a concentration. This can be seen from the relative small difference between the critical nutrient loads for the two hypothetical maximum allowable chlorophyll-a concentrations ( Fig. 1b blue (1) and red line (2)). Finally, the third load-response curve is nonlinear as well, but includes hysteresis that emerges from strong positive feedbacks leading to two alternative stable states (Fig. 1c). Lakes exhibiting a nonlinear load-response curve with hysteresis have two critical nutrient loads, one for eutrophication and one for oligotrophication ( Fig. 1c denoted by CL a and CL b ). In this case, the critical nutrient load is nearly independent of human-defined chlorophyll-a limits, due to the abrupt character of the shift (Fig. 1c blue and red line, with no differences between the critical nutrient loads).
Multiple methods exist to estimate load-response curves and the corresponding critical nutrient loads of lakes. First, field data can be used for hindcasting. For instance, hindcasting has been applied for Lake Veluwe, the Netherlands (Ibelings et al., 2007) and a couple of US Lakes (Baron et al., 2011). While this method gives good results for specific cases, it requires a large dataset to cover the full range of load-responses for both eutrophication and oligotrophication. Such load-response data are scarce in many cases (Baron et al., 2011;Capon et al., 2015). Additionally, hindcasting always lags behind as both eutrophication and oligotrophication need to have occurred to prove the presence of hysteresis. Hence, this method is unsuitable to estimate critical nutrient loads for lakes that are yet to be restored.
Another option is to determine load-response curves and critical nutrient loads through laboratory or field experiments as has been done for Hickling Broad (UK) using mesocosms (Barker et al., 2008). In this method, the critical nutrient load is determined in small containers or cattle tanks in which each batch of replicates receives a certain amount of nutrients. These experiments are relatively easy to perform and have the advantage of replication, however, they may not represent the field situation (Stewart et al., 2013). Indeed, essential feedbacks like the interaction between phytoplankton and macrophytes, or the effect of wind fetch could be easily missed.
Finally, a modelling approach can be used. Modelling has the advantage that most important feedbacks can be incorporated and that scenario simulations can be run nearly 'unlimited' times without affecting the real lake ecosystem. Studies using a modelling approach to determine load-response curves of lakes are numerous (e.g. Janse et al. (2010) and Kong et al. (2016)). However, in these cases only a single load-response curve is identified that should hold for entire lakes. While single load-response curves may be valid for small and homogeneous lakes, in case of large shallow lakes they critically ignore spatial variation in lake characteristics and connectivity .
Lake Taihu (Southeast China) is a good example of a large shallow lake with high spatial variation ( Fig. 2)  . In its pristine state, nutrient loads per unit area of lake surface were below 0.4 gP m À2 yr À1 and 8 gN m À2 yr À1 (Yan et al., 2011). Macrophytes were established at the shores and in the bays, whereas they were absent in the lake's centre due to strong wind forces . In recent decades, the nutrient load rose above 0.93 gP m À2 yr À1 and 19 gN m À2 yr À1 (2012) resulting in excessive phytoplankton blooms that threaten millions of people depending on Taihu (Xu et al., 2015b). These blooms concentrate primarily in the north and centre of the lake, while macrophytes still flourish in the east (Zhao et al., 2013). Clearly, assuming spatial homogeneity is invalid for Lake Taihu. Instead, high spatial heterogeneity in the abundance of macrophytes and phytoplankton suggest variation in the responses to eutrophication within the lake. This variation results in a spatial patterning of critical nutrient loads in Lake Taihu.
In this study, we aim to determine the spatial pattern of critical Fig. 1. Three types of load-response curves: a) linear; b) nonlinear without hysteresis; c) nonlinear with hysteresis (alternative stable states). The horizontal solid lines depict two human-defined restrictions for the maximum allowable chlorophyll-a that are strict (blue, number 1) and less strict (red, number 2). The blue and red dotted line show the corresponding critical nutrient load (CL) (e.g. Scheffer et al., 2001). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) nutrient loads of Lake Taihu, using a novel and comprehensive modelling approach that accounts for connectivity and spatial variation in lake characteristics. We consider determining critical nutrient loads in a spatial context as an important scientific innovation. This implies that the message presented here goes beyond the specific application of our method to Lake Taihu. Our proposed method includes two steps. First, we spatially determine loadresponse curves by applying bifurcation theory for different locations within the lake. Second, for each location the critical nutrient load is defined. Together, these steps result in a map of Lake Taihu, showing critical nutrient loads in a spatial context. This is an important result for lake restoration management of poorly mixed, spatially heterogeneous lakes such as Taihu.

Methods
To determine the different load-response curves and the corresponding critical nutrient loads for different locations throughout Lake Taihu we developed Spatial Ecosystem Bifurcation Analysis (SEBA) using model simulations. Like classical ecosystem bifurcation analysis, SEBA evaluates the equilibrium response of lake dynamics along the axis of external nutrient load (Scheffer et al., 2001). This equilibrium approach has the advantage of largely nullifying the dependence of the model results on the initial conditions on which the opposing transient approach heavily depends. Data to define these initial conditions is rarely available, especially on a spatial scale. Please note that in this study equilibrium refers to reoccurring seasonal patterns for subsequent years, instead of a constant value throughout time which is seen in studies using mini-models (e.g. Scheffer et al., 2007). SEBA differs from classical ecosystem bifurcation analysis by accounting for spatial heterogeneity in lake characteristics such as wind dynamics and depth. Technically, this comes down to substituting of ordinary differential equations for each substance concentration C used in classical ecosystem bifurcation analysis dC dt ¼ f ðtÞ by partial differential equations dCx;y dt ¼ f ðx; y; tÞ . Substitution leads to equations that not only account for bio-physicochemical processes, but also include diffusion and advection processes. Additionally, lake characteristics such as wind and depth will differ spatially. Subsequently, SEBA results in load-responses and critical nutrient loads for different locations throughout the lake.
Applying bifurcation analysis within a spatial context has implications for the definition of the external nutrient load. There are two options to define the external nutrient load: i) as the nutrient supply to a specific part of the lake or ii) as the total external load to the whole lake. The first definition might appear attractive because of its local character. However, this definition has several drawbacks which specially hamper the comparability of the different Map of Lake Taihu with the main inflow (black dots) and outflow (white dots) rivers. Grey dots show rivers that contribute to both in-and outflow, depending on the season. Numbers refer to 1. Zhihugang River, 2. Yugang, 3. Liangxi River, 4. Wangyu River, 5. Longtanggang River, 6. Xujiang River, 7. Guanjin Kou,8. Lianhu Bridge,9. Taipu River,10. Tuanjie Bridge,11. Tiaoxi River,12. Sainli Bridge,13. Hangchang Bridge,14. Yangjiapugang River,15. Changxinggang River,16. Chendonggang River,17. Caoqiao River,18. Yapu Bridge. THL01-THL08 are monitoring stations. load-response curves and critical nutrient loads within one lake (see Appendix A for details). The second definition is identical to the definition used in classical ecosystem bifurcation analysis with simplifies the comparison between both approaches. This definition is also more useful for lake managers because their data and targets relate to the total external nutrient load to the whole lake rather than the local load. Therefore, we consider it most informative to put the lake wide external nutrient load on the x-axis of the load-response curve while putting the local response on the yaxis.
To perform SEBA, we selected the ecosystem model PCLake because it is, as far as we know, the most extensively used food web model for aquatic ecosystems applied for bifurcation analysis. The advantage of PCLake over other eutrophication models is its focus on the competition between three types of phytoplankton (diatoms, green algae and cyanobacteria) and macrophytes . Additionally, zooplankton, fish and benthic organisms are included in PCLake as well as the bio-physicochemical interactions between water and sediment (Appendix B.1). By using such a comprehensive model, we account for a wide variety of essential processes that, besides nutrient load, indirectly influence the shape of load-response curves of shallow lakes. Since spatial heterogeneity is key to our research, the modelling setup requires a spatial approach. Therefore we couple PCLake with the interface Delft3D-DELWAQ using the Database Approach To Modelling (DATM) Van Gerven et al., 2015). The advantage of Delft3D-DELWAQ is the connection with the hydrodynamic modules Delft3D-FLOW and Delft3D-WAVE. We use the hydrodynamics as previously defined and simulated for Lake Taihu on a 3Dgrid (500 Â 500 m grid) (Liu et al. (2013) and Appendix B.2). In this hydrodynamic simulation, Lake Taihu is connected with the 18 most important rivers for in-and outflow ( Fig. 2 and Appendix B.3).
We coupled PCLake at an aggregated horizontal 2D-grid (~1500 Â 1500 m grid) using the default parameter settings as given by (Janse et al. (2010)). Lake specific settings such as sediment and wind characteristics are defined using information from literature (see Appendix B.4 for references and calculations). The full set of state variables, parameters and equations can be found in Appendix C as a DATM-file. We used nutrient load data and mass balances of Taihu (Yan et al., 2011;Janssen et al., 2014) in combination with watershed data to estimate the load contribution by each river (Yu et al. (2007) and Appendix B.5). Since the original resuspension equation of PCLake lost its meaning within a spatial setting, we substituted it by a dynamic equation based on shear stress and wind (Appendix B.6).
Next, the model is validated based on vegetation, chlorophyll-a and nutrient data. First, in an earlier study we have shown that PCLake predicts the spatial distribution of macrophytes in the oligotrophic situation of Lake Taihu well with only 7e24% discrepancy between data and model . Second, in the present study we validate the model based on a dataset of 8 variables for the period 1996-2006 from Sun et al. (2010) (see Fig. 2 for measurement stations). We use R 2 adj and the mean relative absolute error (in the literature referred to as both RE and MARE, here we use RE) for model validation (Equations (1) and (2), Appendix B.7): with n the sample size, p the total number of explanatory variables (including the intercept), d the observed data and m the model output. High R 2 adj -scores indicates good resemblance of the seasonal pattern in the data by the model. Low RE-scores show that the model output falls well within the range of the measurement data. To put our results in perspective, we compared our output with the results obtained by Trolle et al. (2014) using their statistics for R 2 and RE for chlorophyll-a (Fig. 3). The best fit of chlorophyll-a is found for the field stations located in Meiliang Bay (THL01, THL03, THL04, THL06). In the lake centre the model has a poorer fit to the field data (THL07 and THL08). There are two methodological explanations for this mismatch. First the floating phytoplankton could be missed in the sampling procedure resulting in an underestimation of the occurring phytoplankton blooms in the lake centre . A second explanation is a shortcoming in the model. While the model covers the effect of wind on the water movement, the effect on floating biomasses is not included. It is known that floating phytoplankton biomass are produced in the lake centre, but are driven by the wind against water flow in northwestern direction (Wu et al., 2015). In the bays this effect is much smaller due to lower fetch. This shortcoming creates the risk to underestimate the critical nutrient load in the lake centre and overestimate it in at the north-western shores.
Last, the RE-scores for the concentrations for total nitrogen (TN) and phosphorus (TP) are sufficiently low, indicating that the level of total nutrient concentrations is simulated well. However, R 2 adj -scores for the seasonal trends of the nutrient concentrations have their shortcomings. The low fit of the seasonal trend of nutrients can be ascribed to the rough estimations of the seasonal pattern for the nutrient loads since seasonal measurement data were not at our disposal. Additionally, the contribution of the nutrient compounds (e.g. nitrate or ammonium) to the total nitrogen and total phosphorus pools were unavailable to us and needed to be estimated as well. Since we are interested in average yearly chlorophyll-a concentrations that respond to a certain level of nutrients, the effect of the R 2 adj -scores for the seasonal trends of nutrient concentrations will be of less importance than the RE-scores.
With the validated model we performed SEBA. Similarly to classical ecosystem bifurcation analysis, both load-response curves for eutrophication and oligotrophication are determined in order to ascertain the presence of hysteresis. Hysteresis would mean a dissimilarity in shapes between the load-response curves for eutrophication and oligotrophication. First, we construct the eutrophication load-response curve starting with initial settings of an oligotrophic lake (see Appendix B.8 for settings). The initial settings contain a single value per state variable for the entire lake which might raise the question whether this is too limited for our spatial approach. However, since we adopt an equilibrium approach, initiation of the model with the assumption of initial spatially homogenous conditions is valid since each location in the lake ultimately develops to its own specific equilibrium. We impose a yearly constant, but seasonally and spatially variable nutrient load. We run the model to equilibrium, which we define here as a recurring identical yearly pattern while preserving the seasonal variability. We assume equilibrium been reached if 99% of the grid cells have a deviation between the chlorophyll-a output of the last five years of no more than 5%. We repeat the run 50 times with the same initial conditions, though each time with a higher nutrient load than the previous run. In all repeated runs the relative contribution of nutrient load by each river and the ratio between the various kinds of nutrients are kept invariant. For each run, we calculate the yearly average chlorophyll-a concentrations per individual grid cell. Loadresponse curves for eutrophication per grid cell follow from the total nutrient load of each of the 50 runs and the 50 corresponding calculated yearly averaged equilibrium chlorophyll-a concentrations. Next, in order to construct load-response curves for oligotrophication we repeat the same procedure with 50 runs, though starting with initial settings of a eutrophic lake and imposing decreasing nutrient loads. The imposed nutrient load per unit area of lake surface ranges between 0.12 and 3.00 gP m À2 yr À1 with a corresponding N-load of 2.40e60.0 gN m À2 yr À1 (N-load:Pload ¼ 20). The highest nutrient load in our simulations, 3.00 gP m À2 yr À1 and 60.0 gN m À2 yr À1 is nearly a threefold of the current nutrient load (Xu et al., 2015b). Nutrient loads for individual rivers are calculated using the proportional contribution as calculated in Table V in Appendix B.4. The nutrient load of each river is based on concentrations in order to account for seasonality in loads. To speed up runtime performance we run the model in parallel using the High Performance Computing cloud IAAS infrastructure offered by SURFsara (https://www.surf.nl).
The results of SEBA are highly data-rich and therefore a simple representation of the load-response curve for each grid cell as seen in Fig. 1 is inconvenient. Instead, we use two objective, though exclusive indicators that respectively assess the linearity and the presence of hysteresis. These objective indicators distinguish the three types of load-response curves. First, we determine the linearity of the load-response curve of each grid cell which enables us to distinguish the linear (Fig. 1a) from the nonlinear types of loadresponse curves (Fig. 1b and c). A test for linearity might seem superfluous because of the use of the nonlinear model PCLake, however, a linear output would still be possible in case the positive feedbacks are significantly weakened. To determine the linearity we calculate the predictive power of a linear model to resemble the load-response curve. We use the adjusted coefficient of determination (R 2 adj , equation (1)) as an objective indicator for linearity ( Fig. 4a). Essentially, the load-response curves calculated by PCLake depend on a set of mathematical equations that lack stochasticity and are low in computational errors. Therefore, the value of R 2 adj purely results from the deviation of the load-response curve from the linear model. The load-response curve can be regarded as linear in case R 2 adj approaches one, while lower scores indicate a nonlinear load-response curve.
Second, we assess the presence of hysteresis to distinguish the load-response curve of Fig.1c from the load-response curves in Fig.1a and b. Therefore, we look for major differences between the surface area (integral) below the load-response curves for eutrophication and oligotrophication (Fig. 4b). In case of hysteresis, the surface area below the load-response curve of oligotrophication should be substantially larger than the surface area below the load-response curve of eutrophication. From earlier studies we know that such differences are at least 10 mg L À1 (gP m À2 yr À1 ) À1 (Janse et al., 2008).
In case load-response curves are linear to moderately nonlinear, maximum allowable chlorophyll-a concentration needs to be defined according to the societal functions of the lake (Fig. 1). Lake Taihu serves different societal functions, including drinking water, fisheries, irrigation, shipping and flood control. It is inevitable that each societal function has different requirements with respect to chlorophyll-a concentrations. Drinking water requires low chlorophyll-a concentrations whereas higher chlorophyll-a concentrations are allowed for fisheries or irrigation. From earlier studies we know that at 20 mg L À1 chlorophyll-a, bloom formation starts in Lake Taihu which becomes especially harmful if concentrations exceed 40 mg L À1 chlorophyll-a (Xu et al., 2015b). We focus on a range of maximum allowable chlorophyll-a concentrations between 20 and 50 mg L À1 chlorophylla, to cover the full width of bloom formation and beyond. The critical nutrient load of a specific grid cell is calculated as the maximum nutrient supply to the entire lake, for which chlorophyll-a concentrations in this specific grid cell remains below the maximum allowable chlorophyll-a concentration. Fig. 3. Validation of model outcomes for chlorophyll-a concentrations at 7 different stations in Lake Taihu. Ideally, a model has R 2 -scores of 1 and an RE of 0: "optimal model". To put our results in perspective, we compared them with the results obtained by Trolle et al. (2014) for Lake Engelsholm with the three lake ecosystem models PCLake, DYCD and PROTECH.

Results
We found that the indicator for linearity, R 2 adj , keeps largely below 0.4, suggesting that largest parts of Lake Taihu react nonlinearly to eutrophication and oligotrophication (Fig. 5, histogram and load-response curves for location 3 and 4). Yet, there are exceptions to this pattern as for example in Dong Taihu Bay where the indicator for linearity is high (Fig. 5, location 1). Dong Taihu Bay is an isolated area that shows no response at all. Additionally, close to the water inlet at the south-western shore the R 2 adj is higher, corresponding with a slightly nonlinear load-response curve (Fig. 5, location 2).
In our search for hysteresis we compared the surface areas below the two load-response curves for eutrophication and oligotrophication where a large difference would indicate hysteresis. We found a maximum difference between the surface areas of less than 1.3 mg L À1 (gP m À2 yr À1 ) À1 . In most cases we found even less than 1% deviation between the surface areas below the two load-response curves.
These values are considerably lower than 10 mg L À1 (gP m À2 yr À1 ) À1 (Janse et al., 2008). Therefore, it is more likely that this deviation results from small computational errors rather than from hysteresis. In summary, these analyses indicate a nonlinear load-response curve without hysteresis for major parts of Lake Taihu (Fig. 1b) with a few regions that show a linear load-response curve (Fig. 1a). For both types of load-response curves, critical nutrient loads will partly depend on human-defined limits to the maximum allowable chlorophyll-a concentration. We calculated the critical nutrient load for Lake Taihu using a range of possible human-defined maxima for yearly average chlorophyll-a, varying from 20 mg L À1 chlorophyll-a to 50 mg L À1 chlorophyll-a. Results show large spatial variation in critical nutrient loads (Fig. 6). Intriguingly, the pattern of spatial variation depends on the chosen maximum allowable chlorophyll-a concentration. In case of 20 mg L À1 chlorophyll-a as maximum allowable chlorophyll-a concentration, 84% (the area under a normal distribution minus a one-side tail of one standard deviation) of the critical nutrient loads of Lake Taihu are below 0.19 gP m À2 yr À1 (3.8 gN m À2 yr À1 ) (Fig. 6a). The lake centre shows lowest critical nutrient loads indicating that this area is most likely to reach the maximum allowable chlorophyll-a concentration of 20 mg L À1 chlorophyll-a first. Contrarily, Dong Taihu Bay has a higher critical nutrient load indicating that this isolated region is least vulnerable to increased nutrient load. The shores also show higher critical nutrient load especially in the west of the lake. Even at extremely high nutrient loads, a 20 mg L À1 chlorophyll-a threshold is unreachable for about 10% of the regions within Lake Taihu. These regions are indicated in grey in Fig. 6.
Next, we look at the critical nutrient loads for the maximum allowable chlorophyll-a concentration of 30 mg L À1 (Fig. 6b). Model simulations show that 84% of the critical nutrient loads are below 0.27 gP m À2 yr À1 (5.4 gN m À2 yr À1 ). The spatial pattern is similar to the 20 mg L À1 chlorophyll-a with most vulnerable areas in the lake centre. At extremely high nutrient loads, about 15% of the lake will not reach the maximum allowable chlorophyll-a concentration of 30 mg L À1 chlorophyll-a.
With a maximum allowable chlorophyll-a concentration of 40 mg L À1 chlorophyll-a the most vulnerable areas moves towards Gonghu bay and to a lesser extent to Zhushan Bay and Meiliang Bay (Fig. 6c). Less than 80% of the lake area reaches this maximum allowable chlorophyll-a concentration and of these regions the critical nutrient loads are between 0.25 and 0.47 gP. .m À2 yr À1 (5 gN m À2 yr À1 e9.4 gN m À2 yr À1 ). In case of a maximum allowable chlorophyll-a concentration of 50 mg L À1 chlorophyll-a an interesting pattern emerges (Fig. 6d). In our simulations, the lake centre is unable to reach this maximum allowable chlorophyll-a concentration. Instead, critical nutrient loads are reached first at the outer edges of the northern bays including Zhushan Bay, Mailiang Bay and Gonghu Bay. Additionally, the southeast shore reaches the maximum allowable chlorophyll-a concentration of 50 mg L À1 . The critical nutrient loads within the regions that reach the maximum allowable chlorophyll-a concentration of 50 mg L À1 is between 0.32 and 0.78 gP m À2 yr À1 (6.4 gN m À2 yr À1 e15.6 gN m À2 yr À1 ). Only 25% of the lake area reaches the maximum allowable chlorophyll-a concentration of 50 mg L À1 chlorophyll-a.

Discussion
SEBA provided insight in the critical nutrient loads of Lake Taihu by first defining the load-curves and then identifying the critical nutrient load given a certain management target. We now will discuss each of these two steps and compare our findings with historical data, lab experiments and other models. Next, we will reflect on the predictability and applicability of SEBA. Finally we Fig. 6. Critical nutrient load per grid cell for TN and TP (g m À2 yr À1 ) for four maximum allowable chlorophyll-a concentration between 20 and 50 mg L À1 . The grey colors indicate areas that never reached the maximum allowable chlorophyll-a concentration within our simulations. The cumulative graphs below the maps show the area of the lake that reached the maximum allowable chlorophyll-a concentration at a specific critical P-load (TN:TP ¼ 20). will discuss the opportunities of SEBA with respect to methodology, scientific insights and societal relevance.

Step 1: load-response curves
In the first step of our analysis we identified variation in the shapes of load-response curves. Hysteresis as assumed in, for instance, large shallow Lake Apopka (USA) (Bachmann et al., 1999), is not detected within Lake Taihu. Apparently, the lake's surface area is too large, allowing wind forces to dominate and preventing vegetation growth in most part of the lake (Zhao et al., 2013). Additionally, Taihu has heterogeneous lake characteristics which are believed to weaken large-scale hysteresis as well (Van Nes and Scheffer, 2005;Van de Leemput et al., 2015). Hysteresis in isolated areas within Lake Taihu is also absent, possibly due to sufficient interaction between isolated regions and the lake centre . Accordingly, Taihu joins the long list of large lakes without hysteresis Capon et al., 2015;Bunting et al., 2016). Instead, we found nonlinear load-response curves in largest part of the lake (Fig. 1b).
Nonlinearity within load-response curves results from nonlinear interactions (e.g. Scheffer et al., 2001). A well-known non-linear interaction occurs between nutrients and phytoplankton growth. The growth of phytoplankton is limited at extreme low nutrient load due to nutrient limitation and limited at high nutrient load due to light limitation. In contrast, phytoplankton growth is highest at medium loads . Apparently, these nonlinear interactions are strong enough for nonlinear load-response curves to occur, but too weak for hysteresis. At some locations in Lake Taihu, however, the load-response curves are nearly linear. These are regions where positive feedbacks are suppressed e.g. due to strong hydrological influences by rivers that dominate over biological processes (e.g. southwest Taihu), or due to the great distance from the actual nutrient source leading to limited impact of external nutrient loads (e.g. Dong Taihu Bay).
Beside the differences between load-response curve types (Fig. 1), also differences within types can be seen. An important difference is expressed in the maximum chlorophyll-a concentrations reached by the nonlinear load-response curves. Indeed, the maximum allowable yearly averaged chlorophyll-a concentrations of 50 mg L À1 is only feasible in a limited area of the lake. These are regions where nutrient concentrations rise to unprecedented levels while light availability is sufficient for phytoplankton growth as for example in the northern bays. In the middle of the lake it seems that these high averages for chlorophyll-a concentrations cannot be reached due to high suspended matter concentrations (Zhang et al., 2014).

Step 2: critical nutrient loads
In the second step of our analysis we determined the critical nutrient loads for different sections in the lake. The differences between and within load-response types, result in a range of critical nutrient loads for Lake Taihu. Besides, critical nutrient loads depend on the maximum allowable chlorophyll-a concentrations because of the gradual slope in Taihu's nonlinear load-response curves. Therefore, for multiple maximum allowable chlorophyll-a concentrations, a range of critical nutrient loads have been defined.
The low critical nutrient loads in the lake centre are noteworthy in this respect. As mentioned in the validation process, the critical nutrient load has the risk to be underestimated in the lake centre due the absence of wind-driven floating phytoplankton movement within the model. The current version of PCLake only accounts for phytoplankton within the water column; however, phytoplankton can have sufficient buoyancy to float on the water surface making their horizontal movement highly susceptible to wind forces. Adding floating phytoplankton to the model is far from trivial, thus we restrict ourselves to a thought experiment to elucidate the effect of this shortcoming. Intuitively, addition of wind-driven floating phytoplankton movement could improve the critical nutrient load prediction. Likely, this addition will lead to lower predicted phytoplankton biomass in the lake centre where they are skimmed and to higher predicted phytoplankton biomass at the northwestern shores where they accumulate. Hence, the extended model will show higher critical nutrient loads at the lake centre and will show lower critical nutrient loads at the north-western shores. This difference between the potential extended and the present model, points us at interesting aspects in the theory of critical nutrient loading within a spatial context. Indeed, the extended model probably approximates better where the problems become apparent. However, the present model might better indicate where the problems arise. While good localizations of problems are important for symptom treatment, they will not help to localize the problem's origin. It is like mopping the floor, without checking whether the tap still runs.
To check the validity of the critical nutrient loads found in our study, we compare our results with other studies (Table 1). First we made a comparison with historical measurement data. From these data we know that blooms have occurred at nutrient loads of around 0.4 gP m À2 yr À1 and 8 gN m À2 yr À1 (Duan et al., 2009;Yan et al., 2011, with blooms defined as a detectable reflectance signal at the wavelength of 900 nm). Furthermore, previous regime shift analysis with data of Lake Taihu suggest a regime shift between late 1980s and early 1990s when nutrient loads were between 0.38 and 0.58 gP m À2 yr À1 and 7.6e15 gN m À2 yr À1 (Xu et al., 2015a;Yan et al., 2011). These load values exceed our predicted critical nutrient loads for the minimum allowable chlorophyll-a concentration of 20 mg L À1 and fall in the same range as the critical nutrient loads found for the minimum allowable chlorophyll-a concentration of 40 mg L À1 . Our results are thus well in line with historical data.
Next we make a comparison with lab experiments. In the nutrient addition experiment by Xu et al. (2015b) nutrient thresholds were identified. In this bioassay, either nitrogen or phosphorus was added daily to 100 L floating mesocosms in Lake Taihu, containing water with the natural phytoplankton assemblage. After 10 days, Xu et al. (2015b) defined the critical nutrient load as the point where the nitrogen or phosphorus load was sufficiently high that it caused the microcosm chlorophyll-a concentration to exceed the maximum allowable chlorophyll-a concentration of 20 mg L À1 . It appears that our model simulations predict lower critical nutrient loads than found by Xu et al. (2015b). This is especially true for the phosphorus threshold. An explanation for this discrepancy could be neglected lake-size-related processes in the bioassay such as wind forces that act on the sediment. Bioassays could therefore be regarded as "small lakes" of which earlier studies found that critical values are lower (Van Geest et al., 2003;Janse et al., 2008). Estimating critical nutrient load using bioassays have thus the risk of overestimation. The deviation is unlikely caused by the lacking floating phytoplankton transport in the model because higher chlorophyll-a at one location is probably compensated by lower concentrations elsewhere. Additionally the chlorophyll-a estimations for Meiliang Bay in the validation step were reliable.
Finally, we compare our results with the empirical model of Vollenweider, 1975. The phosphorus thresholds found by Vollenweider fall well in the range of what is found in our analysis (20 mg L À1 e30 mg L À1 ). Vollenweider's thresholds for the nitrogen are slightly lower than what have been found in our study (20 mg L À1 e30 mg L À1 ).

Reflections on the predictability and applicability of SEBA
Now that we have compared our findings with field data, lab experiments and other models, we reflect on the predictability of SEBA. SEBA is based on a seasonal equilibrium approach. As mentioned in the methods, this equilibrium refers to reoccurring seasonal patterns for subsequent years. By adopting an equilibrium approach all internal processes including internal loads are in equilibrium. This means that, for instance, there is a seasonal internal nutrient load from sediments, which, at a yearly basis, equals the nutrient uptake by sediments. The net internal nutrient load thus equals zero when adopting a seasonal equilibrium approach. From a transient point of view this seems an imperfection of SEBA. Internal loading is, however, in fact a delayed external load. SEBA thus indirectly accounts for the internal nutrient load.
A major factor that affects the predictability of SEBA is that in reality, seasonal equilibria are unlikely to occur because local and global external dynamics alter the equilibrium point of lake ecosystems. As a result, ecosystems continuously aim for a moving equilibrium which is never reached. Moreover, all kind of stochastic processes keep ecosystems away from their equilibrium state, including weather variations. Despite this fact, the equilibrium approach is commonly used to determine critical nutrient loads of, for instance, lakes (e.g. Scheffer, 1990;Janse, 1997;Janssen et al., 2014). A major reason to adopt the equilibrium approach can be found in the difficulty to define the exact initial conditions of ecosystems on which the opposing transient approach heavily depends. In order to identify the loading threshold of the regime shift using a transient approach, the initial conditions as well as the course of the environmental development need to be exactly defined. However, this information is rarely available. In the equilibrium approach, initial conditions only matter in case lake dynamics show hysteresis. In case of hysteresis and giving a yearly constant nutrient load, there are two equilibrium points to which lakes may develop, either clear or turbid. The equilibrium approach, therefore, requires only a rough estimation of the initial conditions for a highly oligotrophic and a highly eutrophic condition of the lake ecosystem. In this study, the dependence of the critical nutrient load on the initial values of nutrients in sediment was tested by running the model for each nutrient load using both highly oligotrophic initial conditions and highly eutrophic initial conditions. We did not find significant differences between these equilibrium outputs and therefore conclude that the initial conditions have no effect on our results.
Using an equilibrium approach instead of a transient approach limits the predictability of the outcomes. The found critical nutrient loads are indications of the nutrient load at which a shift may be expected to happen. Due to transient dynamics, however, the response of lakes may be different than would be expected based on the identified critical nutrient load. For example, external nutrient loads may require a more drastic decrease to compensate internal loads that emerge at the onset of recovery. On the other hand an equilibrium approach has the advantage that it widens the application of SEBA, since it nullifies the dependence on predefined initial and boundary conditions that delineate a possible future scenario. The outcome of the equilibrium approach used in SEBA thus provides broader applicable guidelines of the critical nutrient loads than a transient approach, however at the expenses of being less precise in the timing of the regime shift.

Methodological opportunities for SEBA
In this study, we introduced the novel technique of Spatial Ecosystem Bifurcation Analysis (SEBA). In contrast to conventional methods such as Vollenweider's models and bioassays, SEBA indicates spatial variation in the critical nutrient loads. Spatial variation is, together with temporal variation and methodological anomalies, a major factor in the uncertainty of the critical nutrient loads (Janse et al., 2010). SEBA elucidates an important part of critical nutrient load uncertainty. The spatial patterns in critical nutrient loads give important insights in the most vulnerable regions of Lake Taihu. We found that vulnerability of regions to nutrient load depends on the maximum allowable chlorophyll-a concentration. With increasing maximum allowable chlorophyll-a concentrations, the most vulnerable region moves from the centre of Lake Taihu towards the borders of the northern bays. Indeed, the 2007's water crisis started at a drinking water station in one of these northern bays Zhang et al., 2010). This is important information for drinking water companies.
Compared to 0D-models used in classical ecosystem bifurcation analysis, SEBA requires longer computational time due to the spatial grid, as well as due to the threefold increase to reach equilibrium due to advection and diffusion. With the present increase in computational power, calculation time is less of an issue.

Scientific opportunities for SEBA
SEBA allows for unique insights in local responses of large lake ecosystems such as Taihu and should be applicable to spatially heterogeneous ecosystems suffering from eutrophication in Table 1 Overview of critical loads per unit area of lake surface for Lake Taihu, calculated by our study and others.

Method
Critical P-load Critical N-load b At a threshold of 20 mg L À1 chlorophyll-a, the concentration at which bloom formation is expected in Taihu according to Xu et al. (2015b). general. First, the wide range of load-response curves and critical nutrient loads shows that simple theory (e.g. Scheffer et al., 2001) breaks down in heterogeneous lakes. Critical nutrient load as estimated by conventional methods such as Vollenweider and bioassays will result in one single value which is not always representative for the whole lake . Indeed one may question the meaning of this value, if the variation in critical nutrient loads is huge. SEBA leads to a range of critical nutrient loads that apply to different sections of the lake. The highest critical nutrient load within the range marks the load at which the first improvements appear. The lowest critical nutrient load marks the load at which the risks on blooms in the entire lake are diminished. Second, dynamics of Lake Taihu can only be understood when connectivity is taken into account. Within Taihu we found relatively isolated regions such as Dong Taihu Bay. In these regions, prevailing processes such as eutrophication and wind forces that govern major part of the lake are locally unimportant. This can explain why macrophytes still prevail in Dong Taihu Bay, notwithstanding the high eutrophication (Zhao et al., 2013). A simple 0D-model simulation based on average values of Lake Taihu likely leads to wrong conclusions for Dong Taihu Bay. Critical nutrient loads can thus only be understood if connectivity is accounted for.
An important next step would be to identify the dependence of critical nutrient loads on the TN:TP ratio of the nutrient load. For simplicity, here we assumed a fixed TN:TP ratio. Xu et al. (2015b), however, assumed P-limitation while defining the critical N-load and vice versa assuming N-limitation to establish the critical Pload. The best solution probably lies in the middle since in nature the TN:TP ratio changes due to natural variability and human activities (Beusen et al., 2016). Information on the dependency of critical nutrient loads on the TN:TP ratio will contribute to the ongoing debate on nutrient control for lakes. In a previous study on Lake Taihu the need for dual nutrient management is expressed (Paerl et al., 2011). In a dual management approach, the TN:TP load ratio is kept constant. Others suggest that the ratio in loads could be used to promote specific desired species as each species has its optimal TN:TP ratio (Anderson et al., 2002).
Finally, SEBA can be applied to other stressors as well. For instance, heavy metals pesticides or pharmaceuticals are major threats to aquatic ecosystems. Like we have done here for chlorophyll-a, critical thresholds for these stressors can be determined spatially as well. Another interesting future step is to look after the effect of climate change to critical nutrient loads as this might be a major factor changing the hydrodynamics of lakes and thereby the spatial patterns.

Societal opportunities for SEBA
The information on the spatial pattern of critical nutrient loads in Taihu will help managers to develop restoration strategies more effectively. Managers have two options. First, in order to prevent any exceedance of the maximum allowable chlorophyll-a concentrations, managers could aim to reduce the present nutrient load below the lowest critical nutrient load found within the entire lake (law of the minimum). For Taihu this will mean a reduction of at least 15% of the present phosphorus load and over 20% of the nitrogen load to meet the standards of the 50 mg L À1 chlorophyll-a. To reduce structural emergence of blooms of 20 mg L À1 chlorophyll-a, nutrient loads need a more drastic reduction of nearly 90% of both nitrogen and phosphorus. However, if little exceedance of chlorophyll-a concentrations is allowed in certain areas of the lake, managers could select for the second option for which the reduction may be less strong. In that case managers could assign specific societal functions to certain parts of a lake. Societal functions with strict requirements to chlorophyll-a concentration can then be assigned to the least vulnerable regions.

Conclusion
SEBA (Spatial Ecosystem Bifurcation Analysis) provides detailed insights in the critical nutrient loads of large shallow lakes. These insights help lake managers in realizing spatially explicit management goals. Indeed it gives a better understanding of the critical nutrient loads throughout the lake, the local vulnerability to harmful phytoplankton blooms and the local recovery trajectories. We conclude that nutrient loads of Taihu have to be more than halved to reach chlorophyll-a concentrations of 30e40 mg L À1 in most sections of the lake. Preventing structural emergence of phytoplankton blooms containing 20 mg L À1 chlorophyll-a throughout Lake Taihu requires a more drastic nutrient load reduction of nearly 90% of both phosphorus and nitrogen. We feel that our approach contributes to the growing societal interest in critical nutrient loads of lakes such as Taihu and of spatially heterogeneous ecosystems in general.