A modeling analysis of spatial statistical indicators of thresholds for algal blooms

Predicting algal blooms both within and among aquatic ecosystems is important yet difficult because multiple factors promote and suppress blooms. Statistical indicators (e.g., variance and autocorrelation) based on time series can provide warning of transitions in diverse complex systems, including shifts from clear water to algal blooms. Analogous spatial indicators have been demonstrated with models and empirical data from vegetated terrestrial ecosystems. Here, we test the applicability of spatial indicators to algal blooms using a nutrient‐phytoplankton spatial model. We found that standard deviation and autocorrelation successfully distinguished bloom state and proximity to transitions, while skewness and kurtosis were more ambiguous. Our findings suggest certain spatial indicators are applicable to aquatic ecosystems despite dynamic physical–biological interactions that could reduce detectable signals. The growing capacity to collect spatial data on algal biomass presents an exciting opportunity for application and testing of spatial indicators to the study and management of blooms.

Algal blooms have large impacts on aquatic ecosystems. In oceans, the spring bloom supports growth of zooplankton and eventual fish production (Mann 1993). In more nutrient-rich lakes and coastal oceans, blooms leading to very high concentrations of algae have adverse effects including accelerated nutrient cycling, depletion of hypolimnetic dissolved oxygen that causes fish kills (Swingle 1968), and toxin release that harms grazers and higher consumers (Christoffersen 1996;Ibelings and Chorus 2007). Blooms may also disrupt aquaculture, recreation, and drinking water supplies (Dodds et al. 2009) and blooms with negative impacts are designated harmful algal blooms, or HABs. High-profile blooms, for example, in Lake Erie, underscore the scientific consensus that HAB occurrence is increasing worldwide (Heisler et al. 2008;Ho and Michalek 2017).
Understanding algal bloom drivers and dynamics is critical to mitigating their negative impacts. Excess nutrient loading, especially phosphorus in inland waters (Schindler et al. 2016), has been widely studied as a main driver of algal blooms. Grazing, temperature, physical mixing, and other factors promote or suppress phytoplankton and thereby affect blooms (Paerl et al. 2001). The varied drivers of algal blooms complicate predictions of the timing and location of blooms. Reliably predicting blooms can improve prevention strategies (e.g., decreasing nutrient loading) and provide time to minimize bloom impacts (e.g., using algicides).
Algae blooms in eutrophic waters are examples of critical transitions (i.e., abrupt shifts in response to small changes in a forcing) in ecosystems (Carpenter et al. 1999;Batt et al. 2013;Cottingham et al. 2015). An emerging body of literature from fields as diverse as financial markets and human physiology has found that generic "early warning indicators" can provide information on system state and proximity to thresholds (reviewed by Scheffer et al. 2015). Early warning indicators are statistics that change in predictable ways as transitions are approached due to critical slowing downslowed recovery from perturbations as a system approaches a transition (Scheffer et al. 2009). Variance and standard deviation (SD) of ecosystem state variables have been widely used as early warning indicators (Carpenter and Brock 2006;Pace et al. 2013) and are expected to increase as a critical transition is approached (Biggs et al. 2009). Decreasing return rates near critical transitions due to critical slowing down also lead to increasing autocorrelation (Dakos et al. 2012). Changes in skewness, another early warning indicator, can capture the dependence of return rate on perturbation direction and response to large external fluctuations as thresholds are approached (Guttal and Jayaprakash 2008).
In ecology, early warning indicators have been studied in a variety of ecosystems using modeling, laboratory, and field experiments as well as historical observations and a few whole-ecosystem manipulations (Scheffer et al. 2015). Most studies, especially in aquatic ecosystems, have focused on temporal statistics as early warning indicators (Carpenter et al. 2011;Pace et al. 2017;Wilkinson et al. 2018). Analogous changes in spatial statistics are also expected from critical slowing down near transitions and can potentially provide information on ecosystem state and proximity to thresholds without extensive prior data (Dai et al. 2013). In terrestrial ecosystems including grasslands and shrubland, spatial statistics of vegetation coverage change near critical transitions driven by precipitation and fire frequency (Kéfi et al. 2014;Ratajczak et al. 2017). Spatial early warning indicators have not been studied as widely in aquatic ecosystems, but there are a few examples where changes in spatial statistics are associated with transitions. Donangelo et al. (2010) extended a simple eutrophication model to a spatial grid and found pattern formation near the transition to a degraded state, as well as an increase in spatial variance before temporal variance. Litzow et al. (2008) and Cline et al. (2014) observed changes in fish distribution and spatial variance during regime shifts in the fish community. Rindi et al. (2017) found that the spatial recovery length from perturbations, another indicator of critical slowing down, increased in the rocky intertidal near a transition from canopy to turf dominated benthic algae. A recent study used spatial mapping of phytoplankton pigments to test for indicators of transitions during a whole-lake nutrient manipulation (Butitta et al. 2017). The utility of spatial early warning indicators in these studies is promising, but additional modeling tests with precisely known system dynamics and transitions can provide understanding of when spatial indicators are expected to give reliable information on ecosystem state and threshold proximity and guide future research.
The spatially patchy distributions observed in algal blooms and captured by spatial models (Franks 1997) suggest that spatial early warnings may also be useful for analyzing and managing aquatic ecosystems. It is unknown if the highly dynamic (time-varying) spatial distribution of algal blooms could limit the possibility that spatial patterns and statistics provide early warning of critical transitions. Here, we ask (1) can spatial indicators be used to distinguish between states (e.g., bloom vs. nonbloom) in aquatic ecosystems that undergo critical transitions? and (2) do spatial indicators provide reliable warning of approaching algal bloom critical transitions?
To answer these questions, we applied a published spatial model of algal blooms with known critical transitions (Serizawa et al. 2008). We evaluated the ability of spatial early warning indicators to distinguish algal bloom states and transition proximity in aquatic ecosystems by simulating the model through time for a range of nutrient input levels and then calculating spatial early warning statistics.

Methods
We use a model created by Serizawa et al. (2008) (see Supporting Information Appendix A for complete description).
The model is defined on a two dimensional spatial grid (180 × 180 cells) and represents dynamics of phosphorus and phytoplankton in pelagic systems, which interact nonlinearly via Holling type-II responses. Neighboring grid cells exchange components physically via diffusion and advection, the latter generated by randomly seeded eddies. Like Serizawa et al. (2008), we use the nondimensional version of the model. This transformation reduces the number of parameters and does not affect the relative magnitude of the resilience indicators (see eqs. 3 and 4 in the Supporting Information for transformation between nondimenional and dimensional values). Because this model is dimensionless, the spatial components (extent, resolution, diffusion, and advection) of the grid are defined relative to each other. We have chosen this framework as a basis for exploring spatial indicators conceptually with a model that is reasonably complex but also tractable. The model framework creates patterns that are highly dynamic both spatially and temporally and includes advection and diffusion. Advection is an important spatial flux in pelagic ecosystems and has not been considered in the previous studies of spatial resilience indicators in terrestrial ecosystems. The original model was modified to include stochasticity in the phytoplankton dynamics to represent processes not explicitly included in the model and environmental noise, for example, local variations in nutrient inputs or weather (eq, 1, Supporting Information Appendix A). Critical slowing down and resilience indicators are based on a system's response to such perturbations changing as a transition is approached. Simulations and calculations were carried out in R 3.4.1 (R Core Team 2017); code is available on GitHub at https://github.com/cbuelo/ SpatialBloomIndicators.
The model exhibits critical transitions (Hopf bifurcations, where stable fixed points transition to limit cycles or vice versa) in phosphorus-phytoplankton dynamics as the phosphorus input rate, i, is varied. These critical transitions are functions of the model parameter values and define the bloom states, which Serizawa et al. (2008) identify using stability analysis. All analyses use the parameter values from Serizawa et al. (2008, Table 2, set I), and included in the Supporting Information Table S1. At low i (0.3 ≤ i < 0.5), the system is in a constant, low-phytoplankton, nonbloom stable state. As i increases, the system enters a stable limit cycle at the low-input transition (i ≈ 0.5); this cycling bloom state is characterized by repeating cycles in phosphorus and phytoplankton concentrations. At i above the high-input transition (i ≈ 1.25), the system returns to a constant-bloom stable state with high phytoplankton concentrations (Serizawa et al. 2008, fig. 4b).
Integration was carried out using the Euler-Maruyama method with a small time step (0.025) to approximate continuous development of the model (Higham 2001). The phosphorus input rate was constant for each simulation. Separate simulations, each with a fixed and distinct phosphorus input rate, were compared to assess the ability of spatial indicators to distinguish different states or indicate proximity to thresholds. We define "reliable warning" as an indicator for which neighboring bloom states have distinct indicator values and, within a given state, the indicator changes unambiguously as a transition is approached. For each phosphorus input rate simulation, the model was first run deterministically as in Serizawa et al. (2008; and see Supporting Information). The longrun (t = 1000) deterministic state was used as the starting state for stochastic simulations. Stochastic simulations were run for 500 time units, with system state and spatial statistics retained every 1 time unit after a spin-up period of 100 time units. For each of the 400 "snapshots" at a given phosphorus input rate, the following spatial statistics were calculated from each grid's phytoplankton concentrations: mean, SD, skewness, kurtosis, Moran's I, and autocorrelation range (AC range). We tested kurtosis as a potential spatial indicator, expecting that there may be a changing proportion of extreme values as the transitions approached. Moran's I measures the degree of correlation between neighboring grid cells (analogous to lag-1 temporal autocorrelation), while AC range is the maximal distance over which grid concentrations are correlated; these autocorrelation measures explicitly depend on the spatial distribution of phytoplankton concentrations. The other measures are sample statistics of the 180 × 180 values in each "snapshot" and do not depend on the specific spatial distribution of the values on the grid. Skewness and kurtosis were calculated using the moments package in R. AC range was computed from the semivariogram using the gstat package function fit.variogram() and exponential fit. We evaluated the robustness of our findings by varying the physical forcing. In extreme cases, high diffusion and low advection would eliminate spatial patterns and simulate 180 × 180 synchronous cells. Relative to the base case from Serizawa et al. (2008), we repeated the simulations for two additional cases: a high diffusion case (double the base case diffusivity) and low advection case (half of the base case advection velocity). Increased diffusivity decreased local concentration gradients and decreased advection slowed patch formation.

Results
Model simulations generated the expected spatial distributions: spatially uniform low and high phytoplankton concentrations in the low and high phosphorus input stable states, respectively, and patchy patterns of high and low concentrations in the intermediate-input cycling state (Fig. 1). The spatial patterns in the cycling state changed and repeated through time, as in the study by Serizawa et al. (2008; and see animated Fig. S1 in Supporting Information).
For the stable states at both low and high phosphorus input rates, spatially uniform phytoplankton concentrations were maintained through time and had relatively constant mean grid phytoplankton concentration and spatial indicators at a given input rate (Fig. 2). In the intermediateinput rate cycling bloom state, temporally cycling spatial patterns in phytoplankton concentration resulted in cycling spatial statistics. Both the grid mean and spatial indicators exhibited repeating cycles in time for a given phosphorus input level.
All spatial indicators, except for Moran's I, were highly variable in the intermediate-input cycling bloom state and more constrained in both the nonbloom and constant-bloom stable states (Fig. 3). Mean phytoplankton concentration increased across the range of phosphorus input rates; however within the cycling bloom state, there was significant overlap in the distributions of grid mean phytoplankton concentration between adjacent input rates (Fig. 3A). Both (SD) and skewness had "humped" patterns within the cycling state where indicator values peaked in the middle of the phosphorus input range and declined near the transitions (Fig. 3B,E). In the cycling state, SD values had less overlap between adjacent input rates near the transitions and no overlap with the stable states, whereas skewness had a large degree of overlap especially at the high input transition. Autocorrelation range (AC range) had the opposite pattern within the cycling state; values were highest and variable near the transitions and decreased at phosphorus input rates in the middle of the cycling state (Fig. 3C). Moran's I was the only indicator that was highly constrained (near 1) in the cycling bloom state and was also the only indicator to decrease steadily within increasing input rates in the constant-bloom stable state (Fig. 3D). Kurtosis, while constrained in the stable states, was highly variable in the cycling bloom state with a high degree of overlap and no trends near the transitions (Fig.3F).
Increasing diffusion or decreasing advection did not have a strong effect on the statistical moment spatial indicators (mean, standard deviation, skewness, and kurtosis) in any of the bloom states (Fig. 3A, B, E, F, respectively). Increasing diffusion increased both autocorrelation indicators relative to the base case. AC range was most strongly effected in the cycling bloom state (Fig. 3C) while the increase in Moran's I was stronger in the nonbloom and constant-bloom stable states (Fig. 3D). Decreased advection had minimal impact on all statistics.
To determine whether spatial indicators can reliably distinguish proximity to thresholds, we selected the spatial indicators from Fig. 3 that were best at distinguishing between bloom states and studied them at a higher resolution of phosphorus input rate near the transitions. Bloom state could be unambiguously inferred from SD, AC range, and Moran's I at phosphorus input rates near the critical transitions (Fig. 4). SD increased steadily with phosphorus input rate through the transition from the low-input stable state to intermediateinput cycling state, with little or no overlap between distributions of SD at adjacent input rates (Fig. 4A). At the transition from the cycling bloom state to constant-bloom stable state, SD decreased and was relatively constant at phosphorus input rates above the transition (Fig. 4B). Autocorrelation range was low and overlapped at input rates from 0.3 to 0.4, increased sharply at i = 0.45 before the low-input transition occurred, then gradually declined from i = 0.5 to i = 0.7 (Fig. 4C). AC range increased gradually with overlapping distributions as the high-input transition was approached at i = 1.25, then fell sharply and overlapped from i = 1.3 to 1.45 (Fig. 4D). Moran's I also increased sharply at i = 0.45 below the low-input transition and was high (near 1) in the cycling bloom state (Fig. 4E). Moran's I remained near 1 as the high-input transition was approached in the cycling bloom state and declined slightly at i = 1.25 before falling sharply at i = 1.3 (Fig. 4F). From i = 1.3 to i = 1.45 in the high-input stable bloom state, Moran's I declined steadily with relatively little overlap among adjacent input values.

Discussion
Spatial indicators can discern ecosystem state and proximity to thresholds, even in a model of a highly dynamic pelagic system with interacting physical and biological components. SD, AC range, and Moran's I differed between the cycling and stable states, allowing for classification of bloom state (nonbloom stable, cycling bloom, or constant-bloom stable) from a single spatial snapshot. At phosphorus input rates near both transitions in bloom state, trends in these indicators provided a reliable warning of approaching transition, indicating that changes in spatial indicators calculated from repeated sampling in time would occur prior to crossing thresholds. Skewness and kurtosis had somewhat distinct distributions between states but significant overlap within states, suggesting that these indicators would be unreliable for determining threshold proximity.
Spatial SD was one of the most robust indicators of both ecosystem state and threshold proximity. An approximately order of magnitude difference in SD on either side of the transitions clearly indicated the state of the system. In addition, the monotonically increasing and nonoverlapping distributions prior to the low-input transition provide a relative measure of proximity to the approaching threshold. At the high-input transition from the cycling bloom state to the constant-bloom stable state, there was a slightly less than order of magnitude difference in SD. The decline in SD as this transition was approached from lower phosphorus input levels may provide an indication of the approaching threshold, although there was some overlap between adjacent phosphorus input levels. The observed increases in SD as transitions were approached from the stable states match with findings of temporal early warning indicator studies (Carpenter and Brock 2006;Carpenter et al. 2011;Pace et al. 2013). The maximal values of SD occurred at intermediate input levels within the cycling state, corresponding to the maximum cycle amplitude in phytoplankton concentration ( Fig. 3 this study; fig. 4b from Serizawa et al. 2008).
Both measures of autocorrelation were indicative of bloom state and threshold proximity. AC range and Moran's I provided warning of the low input transition from the stable non-bloom state to the cycling bloom state, increasing sharply from i = 0.4 to i = 0.45 prior to the transition at i = 0.5. The steady decrease in Moran's I at input rates above the high-input transition combined with a fairly small degree of overlap in this statistic between adjacent input rates could provide early warning of transition from the stable, constantbloom state bloom state to the cycling bloom state (i.e., if input rates were high but decreasing, Fig. 4F). Moran's I provides clearer warning of proximity to transition to the cycling state (from either low or high input rates) than AC range. In contrast, AC range provides more information than Moran's I on proximity to transitions from within the cycling state but is more ambiguous in the high-input stable state. The differences between AC range and Moran's I are perhaps not surprising as the interactions between physical processes, biological processes, and return rate from perturbations have different effects on correlation between neighboring cells (Moran's I) and those further apart (AC range). The minimum in AC range within the cycling state likely corresponds to a maximum in return rate at the phytoplankton cycle amplitude maximum. The increase in AC range in the high diffusion case may be the result of decreased local gradients in phytoplankton concentration near patch edges.
While SD and Moran's I changed unambiguously as at least one transition was approached from at least one direction, they also displayed statistic-dependent differences in sensitivity to the low-and high-input transitions, depending on the direction from which the transition is approached (SD, Fig. 4B; Moran's I, Fig. 4E,F). While critical transitions are most frequently studied in the context of ecosystem degradation (e.g., increased nutrient inputs), we also observed early warnings in the opposite scenario which may be useful in systems undergoing remediation (e.g., Lake Washington, Hampton et al. 2006). Asymmetric warnings prior to transitions (i.e., whether approached from lower or higher input rates) have also been observed in studies of temporal resilience indicators using models with Hopf bifurcations (Batt et al. 2013). In our study, warning asymmetry may result from an interaction between the deterministic dynamics of the model and the structure of the stochasticity (Horsthemke and Lefever 1984;Benincà et al. 2011). These results reinforce previous suggestions that a resilience indicator approach should only be used when there is reason to suspect a critical transition and data on appropriate state variables are collected (Wilkinson et al. 2018, Gsell et al. 2016. In this study, both skewness and kurtosis were highly constrained in the stable states and much more variable in the cycling states. Despite these differences between bloom states, the lack of trends near transitions and overlap of the indicator distributions in different states limit the applicability of skewness and kurtosis as spatial indicators for this model. Contrasting with our results, Guttal and Jayaprakash (2009) found spatial skewness an unambiguous indicator of an impending regime shift using a two dimensional model of terrestrial vegetation collapse. We suspect these differences arise from the different model types used in each study; their model includes a single fold-bifurcation separating two spatially homogenous states while the model used in this study has a spatially patterned intermediate state separated from stable states by Hopf bifurcations. Overall, studies of spatial resilience are less numerous than studies of temporal resilience and more work is needed to understand differences among types of ecosystems.
This work demonstrates the potential applicability of spatial resilience indicators to pelagic aquatic systems. Spatial statistics have been successfully tested as resilience indicators in arid terrestrial ecosystem models containing critical transitions in vegetation state (Reitkerk et al. 2004;Kéfi et al. 2007;Guttal and Jayaprakash 2009;Dakos et al. 2011). Recent empirical tests in terrestrial ecosystems including experimental grassland manipulations (Ratajczak et al. 2017) and applications of remote sensing (Eby et al. 2017) have supported model findings. Many of these terrestrial vegetation models include diffusive exchange between neighboring grid cells. Compared to rooted and relatively stationary vegetation on land, the physical forcing of pelagic environments (e.g., both diffusion and advective currents) plays a more immediate role in determining the movement of freefloating phytoplankton and nutrients. These forces could "wash out" signals of aquatic ecosystem state and transitions. However, a few empirical studies have shown promising results in the application of spatial resilience indicators to aquatic ecosystems with organisms that have more control over their spatial pattern than phytoplankton, including fish (Litzow et al. 2008;Cline et al. 2014) and intertidal benthic algae (Rindi et al. 2017).
It is also possible that the rapid time scales on which algal blooms develop and spatial patterns form and change could limit the use of spatial resilience indicators. However, a recent empirical study applied several of the indicators used here to spatial data from an experimentally fertilized lake (Butitta et al. 2017). Butitta et al. (2017) found that SD was highest prior to and during an induced algal bloom and declined after fertilization ended and the bloom abated. AC range was highest prior to and after the bloom peak but declined for a short period coincident with peak bloom conditions. Both the SD and AC range observations of Butitta et al. (2017) are consistent with the findings in this study of a shift from the nonbloom stable state past the low-input transition into the cycling bloom state, and then a return to the stable nonbloom state when fertilization ceased. Contrasting with our findings, Butitta et al. found that skewness was elevated prior to and during the bloom and maximal just after the bloom peak. It is not possible to compare spatial statistics in the constant-bloom stable state as Butitta et al. (2017) stopped nutrient additions after temporal early warning signs were observed.
While our findings suggest that spatial indicators can differentiate bloom state and threshold proximity using a model that creates dynamic patches, evaluation and application of spatial indicators requires further research. The model could be expanded to include the vertical dimension of the water column and more ecological and physical detail (e.g., additional nutrients, light limitation, spatially and temporally varying zooplankton density). Most importantly, the model could be used in tandem with field experiments to improve understanding of spatial dynamics of blooms under diverse conditions. Empirical studies could also address questions related to the spatial and temporal scales at which spatial indicators occur, for example, do spatial indicators change prior to blooms and early enough to be useful for management? What spatial resolution and temporal sampling frequency is required? How do temporal and spatial indicators compare?
The increasing ease and cost-efficiency of collecting spatial data on phytoplankton biomass provides the opportunity to test and improve the predictions in this study. Technological advances have made it possible to spatially map photosynthetic pigments both in situ (Crawford et al. 2015) and via remote sensing (Tyler et al. 2016). We advocate additional work using modeling, empirical, and combined studies to further develop and test the applicability of these methods. Successful use of spatial indicators could allow the classification of the current state of individual aquatic systems and their proximity to thresholds. These goals are particularly important as continued nutrient inputs to inland waters and coastal marine systems increase the potential for eutrophication and harmful algal blooms (Sinha et al. 2017).