Bimodality and alternative equilibria do not help explain long-term patterns in shallow lake chlorophyll-a

Since its inception, the theory of alternative equilibria in shallow lakes has evolved and been applied to an ever wider range of ecological and socioecological systems. The theory posits the existence of two alternative stable states or equilibria, which in shallow lakes are characterised by either clear water with abundant plants or turbid water where phytoplankton dominate. Here, we used data simulations and real-world data sets from Denmark and north-eastern USA (902 lakes in total) to examine the relationship between shallow lake phytoplankton biomass (chlorophyll-a) and nutrient concentrations across a range of timescales. The data simulations demonstrated that three diagnostic tests could reliably identify the presence or absence of alternative equilibria. The real-world data accorded with data simulations where alternative equilibria were absent. Crucially, it was only as the temporal scale of observation increased (>3 years) that a predictable linear relationship between nutrient concentration and chlorophyll-a was evident. Thus, when a longer term perspective is taken, the notion of alternative equilibria is not required to explain the response of chlorophyll-a to nutrient enrichment which questions the utility of the theory for explaining shallow lake response to, and recovery from, eutrophication.

If we examine the relationship between both TP, TN and chlorophyll-a we can see that whilst some of the time series show some signs of bimodality in the Chl-a kernel density plots (Supplementary Figure. 1) the correlation between TP, TN and Chlorophyll-a are generally very linear, including the lakes "3000005" and "4500005". There are, however, some anomalous sites, such as lake "4200003" where there is very low chlorophyll-a for a given TP. This is, in fact, the same outlier dots visible in the TP -chlorophyll-a-a scatter plots of the main text (Fig. 2).
This different relationship between TP and chlorophyll-a could be argued to be a sign of alternative equilibria, however, a simpler explanation for this difference lies in the fact that the lake "4200003" is an extreme outlier in terms of the molar TN:TP ratio (Supplementary Figure. 2,bottom panel). Of all lakes with long-term data available, this lake is the only case with molar TN:TP ratios are well below Redfield ratio. This implies not TP but TN limitation of phytoplankton growth, which may explains its outlier status when only TP and Chlorophyll-a are considered (Supplementary Figure.  The simulation approach is described in the methods and is designed to assess the effects of the efficacy of the indicators used to identify the presence, absence or prevalence of the alternative stable states in simulated and real world data. The simulations were constructed to resemble real world data as much as possible, as such they are based on the true range of real nutrient data and the variability of the Chlorophyll-a response and for 99 lakes (five-year mean data consists of 99 lakes, see methods for details). We chose a constant number of years per lake (20 years per lake, which equals the maximum time series length of the real data) to better compare the different alternative-stable state scenarios, without adding the confounding factor of time-series length (variable time-series length is used in the scenarios in the main text to be as close to the real world data as possible). We used the following scenarios: 1. No alternative stable states (main text Fig.3 and Supplementary Figure.   further corroborate this, where R² increases with higher numbers of mean years. Furthermore, we find a more unimodal distribution of model residuals and a unimodal distribution of kernel density plots with increasing numbers of years in the mean calculation. In combination, these proxies are diagnostic of the absence of alternative stable states in the simulated data. The reason for the higher R² and the unimodality of the residuals and the response is the underlying linear relationship between the simulated (nutrient) and (Chl-a) over multi-year time scales, which is blurred by interannual short-term variability.
In contrast to a simple linear response of to (no alternative stable states), the unstable alternative simply increases the year-to-year variability of the response. This is reflected in a slightly lower R 2 compared with the scenario where alternative stable states are absence and a slightly broader spread in the kernel density plot.
However, due to the instability of the two states in this scenario, the response still follows a linear pattern.
Supplementary Figure 3: Scenarios 1 (no alternative stable states) and 2 (unstable alternative states). Example of a single time series from the simulated dataset of 99 lakes (grey dots are single-year values & black dots are five-year rolling means) and proxies for the entire simulated dataset. The bootstrap procedure for each of the simulation scenarios ran with 500 iterations. The scatter plots contain only the first 120 iterations of the bootstrap procedure for better overview. The distributions of R 2 , residuals and kernel density pots based on these multiple runs are shown as density plots with the error bars being the highest density interval with 95% credible interval.
Supplementary Figure 4: Scenarios 3 (alternative stable states within one time series) and 4 (alternative stable states in time series from different lakes). Example of a single or two time series from the simulated dataset of 99 lakes (grey dots are single-year values & black dots are five-year rolling means) and proxies for the entire simulated dataset. The bootstrap procedure for each of the simulation scenarios ran with 500 iterations. The scatter plots contain only the first 120 iterations of the bootstrap procedure for better overview. The distributions of R 2 , residuals and kernel density pots based on these multiple runs are shown as density plots with the error bars being the highest density interval with 95% credible interval.
For scenario 3 and 4, the example time series also show how the 5-year means smooth the year-to-year variation, making the underlying alternative stable states more visible within the time series of a single lake and in the comparison of two different lakes (Supplementary Figure. 4, uppermost panels). Here, the multimodality of the density distribution of the five-year means for (representing two different, stable Chl-a responses to the same nutrient concentration) already suggest alternative stable states Supplementary Figure. 4, uppermost panels).
In contrast to the scenarios without alternative stable states, the scatter plots reveal a clear two-level response of (simulated Chlorophyll-a) response to the same (simulated nutrient concentration) with a higher number of mean years. Furthermore, the R² of the GLMs do not increase as the number of years represented in the means increases, reflecting the two distinct linear responses of y (chlorophyll-a-a) to the x (nutrient). Finally, the residuals of the GLMs reveal clear multimodality and the high-bandwidth kernel density plots do not converge to a normal distribution with a higher number of mean years, further corroborating the R² proxy and scatter plots.
Thus, the use of this gradient of simulations containing different degrees of alternative equilibria show that the diagnostics of the R 2 , residual pattern and kernel density plot have patterns that are characteristic of both the presence, absence and even prevalence of alternative equilibria in the data. Thus, the combination of proxies can be used to test both the presence and absence of alternative stable states in real world data.

Combining data from different lakes with our bootstrap approach results in better detection of alternative stable states than traditional approaches
In combination with the proxies, the hierarchical bootstrap approach has another advantage. As shown in scenario 4, with our approach we find alternative stable states even if the alternative stable states appear in the time series of different lakes. Such a scenario is especially likely for short time series, where the entire time series may only be in one alternative state. We often have short time series in the data (in the dataset with five-year means, most of the time series (89 out of 99) are ten years long or even shorter, see also the analysis of long-term real data below).
Hence, our approach will detect alternative stable states in cases where an independent investigation of time series from single lakes might fail, which makes it a powerful means to detect both the presence and absence of stable states in multi lake data sets of variable length.
Supplementary note 3: Alternative stable state assessment for real data with limited data range.
The existence of alternative equilibria in shallow lakes is generally considered possible within a limited range of nutrient enrichment. Thus, it could be that the alternative stable states present in the data are obscured when the full dataset is analysed. In order to check this we repeated the analysis on a more limited TN and TP concentration range.
We filtered and re-analyzed the data, only keeping data points within the following two ranges: -TN concentration = 0.5 -2 mg / L -TP concentration = 0.05 -0.4 mg / L.
In the filtered data, 1329 out of the original 2876 single-year data points, 289 out of 1028 three-year mean data points and 212 out of the 864 five-mean year data points remained. The filtered data consisted of data points from 550, 48 and 27 lakes for the single-year data, three-year means and five-year means, respectively.
The smaller range resulted in lower R² of the models, but the patterns revealed by the multi-year means result in higher R² compared with single-year data was largely consistent. The exception being the five-year mean TN models for which both, the single-year and mean data resulted in very low R² (Supplementary Figure. 5). This can be explained by the fact that at the restricted nutrient range, TP exerts more influence on chlorophyll-a than TN and this becomes more evident at the multi-year scale. This difference in influence of nutrients over different portions of the enrichment gradient has been demonstrated by recent work using boosted regression trees on the Danish part of the data set 2 , where the predictive power of TP of chlorophyll a was highest in the low range of TP enrichment whereas TN becomes more important at higher levels of enrichment. The filtering of the data to have reduced range of enrichment cut out the portions of the gradient where the influence of TN is greatest. Due to the lower number of samples, the errors of all proxies were higher, reducing the clarity of the results compared with the full data set. Still, we see no indication of alternative stable states in the scatter plots (two groups of dots are not appearing, year and five year means) this data follows the patterns of the simulated data set with no alternative equilibria, as the Kernel density plots indicate no bimodality, the R 2 of the models increase in their explanatory power, though the multi-year mean data only reliably explains more variance than its single year equivalent for TN and there is no bimodality in the residual of the models. The distributions of R 2 , residuals and kernel density plots based on these multiple runs are shown as density plots with the error bars being the highest density interval with 95% credible interval.
Results of the analysis from only U.S. lakes (Supplementary Figure. 8) shows results of the diagnostic tests i) R 2 ii) patterns in models residuals and iii) shape of the kernel density plot the temporal scale increases (single year, three year and five year means). These data follow the patterns of the simulated data set with no alternative equilibria, as the R 2 of the models increase in their explanatory power and the multi-year means always explain more variance than their single year equivalent and there is no bimodality in the residual of the models or in the kernel density function plot of chlorophyll-a.
Supplementary Figure 8: Results of the analysis from only U.S. lakes showing results of the diagnostic tests i) Kernel density plots, ii) change in R 2 and iii) patterns in models residuals as the temporal scale increases (single year, three year and five year means) this data follows the patterns of the simulated data set with no alternative equilibria, as the Kernel density plots indicate no bimodality, the R 2 of the models increase in their explanatory power, though the multi-year mean data only reliably explains more variance than its single year equivalent for TN and there is no bimodality in the residual of the models. The distributions of R 2 , residuals and kernel density plots based on these multiple runs are shown as density plots with the error bars being the highest density interval with 95% credible interval.