Extreme events in lake ecosystem time series

Climate change has generated growing interest in extreme events, and extremes are known to have important consequences for ecosystems. Theoretical mechanisms generating frequent extremes apply to both environmental and biological processes, yet past studies of ecological extremes have focused primarily on the abiotic environment. The rarity or commonness of extremes in biological time series is unknown. We evaluated the statistical tendency to produce extreme events in 595 biological, chemical, physical, and meteorological time series taken from 11 lakes. We found that extreme events were much more likely to occur for biological variables than for other categories. Additional analysis revealed that the tendency to produce extremes was driven primarily by within‐year dynamics, suggesting that processes occurring at short time scales underlie the high frequency of extremes in biological variables. These results should lead us to expect surprises in long‐term observations of biological populations.

time series (Palutikof et al. 1999;Katz et al. 2002Katz et al. , 2005. These large events can be minima or maxima, but henceforth we restrict our discussion to maxima because ecosystem variables are usually defined over the interval [0, 1). The "tailedness" of a distribution refers to the shape of the tail of a probability distribution, affecting its ability to produce frequent extremes and to generate record-breaking events. We measured tailedness by fitting the generalized extreme value distribution (GEV) to observations. The GEV has three parameters, location (l), scale (r), and shape (n). n measures tailedness: negative n indicates a bounded upper tail, n near 0 indicates a thin (exponential) tail such as demonstrated by a normal distribution, and positive n indicates that the maxima are fat-tailed.
Fat tails have been found in many types of data, although there is little empirical description of the tailedness of time series of most ecosystem variables, particularly biological variables. Anecdotally, ecologists often recall large events having been observed in their study systems, but rigorous quantification of the size and frequency of these events is lacking. Tailedness has been evaluated for environmental time series like flood stage (Villarini and Smith 2010), stream discharge and precipitation (Katz et al. 2002), sediment deposition (Katz et al. 2005), water temperature (Gaines and Denny 1993), wind speed (Gaines and Denny 1993;Palutikof et al. 1999), and wave force (Gaines and Denny 1993;Denny et al. 2009). Biological time series, by contrast, are rarely analyzed for tailedness. In the case of Segura et al. (2013), tailedness was calculated for the changes in population size of microbes, but not for population size itself. Schmitt et al. (2008) calculated the tailedness of the size of a single population of marine zooplankton, finding the population to be fat-tailed. We could not find more examples quantifying tailedness in biological time series. However, extreme events produced by fat-tailed time series are known to have significant impacts in other settings, making it important to identify their prevalence in biological variables.
Despite the paucity of studies of tailedness in biological variables, it is plausible that biological time series could be fat-tailed. One way for a biological time series to be fat-tailed is for a population to amplify the variability of the surrounding environment. As a heuristic example for how such enhancement may occur, consider a population governed by birth and survival processes that are in turn influenced by environmental (stochastic) drivers (Supporting Information: Extremes in a model population). Birth and death processes often act multiplicatively (Pielou 1977;Gurney and Nisbet 1998), which causes the population to have thicker tails than the environment (Fig. 1). The ubiquity of multiplicative processes in biological systems makes this mechanism for amplifying tailedness potentially widespread, although other mechanisms can generate fat tails such as some types of nonlinear dynamics (Sornette 2006;e.g., chapter 14).
Here, we document the fat-tailedness and time series characteristics of biological and environmental variables to better understand the frequency and origin of extreme events in an ecosystem context. Our goal is to fill the knowledge gap about how often biological variables may be fattailed, and how the tailedness of biological variables compares to the tailedness of their environment. We surveyed 595 long-term lake time series to quantify the relative frequency of extremes in biological, chemical, physical, and meteorological variables. These variables were measured within the same ecosystems, providing a direct comparison between biological tailedness and the tailedness of environmental variables, many of which are expected to affect variation in the biological variables. We then analyzed these time (D-F) Empirical density plots for the parent distributions in A-C (gray lines) and for annual maxima (blue, red, and black lines). (G) The blue, red, and black lines show density plots of n estimated from 500 simulations of the process exhibited in A-C. From theory, as the number of years becomes large, the estimates of n X will always exceed n u and n c , although this is not always the case for the simulated 30 yr of data. series to determine if they showed evidence for generating fat tails at time scales shorter or longer than the sampling interval. This comparison constitutes a direct empirical assessment of the tailedness of a large number of biological time series, and may provide a clue as to if and how biology amplifies variability in exogenous drivers to produce extreme events.

Data collection
Raw data were downloaded from the online databases of the North Temperate Lakes Long-Term Ecological Research (NTL-LTER) program (lter.limnology.wisc.edu/; Supporting Information). Data were downloaded on 18 April 2013, except wind data, which were downloaded on 30 March 2014. Most time series from the northern Wisconsin lakes begin in 1981. Most time series from the southern lakes begin in 1996, except fish time series that begin in 1995, zooplankton that begin in 1976, and three ice cover time series that begin in 1851, 1852, and 1877. Most time series end between 2010 and 2012. Measurement frequencies varied among time series; meteorological time series were typically daily, fish time series were annual, and other lake time series were quarterly or every 2-6 weeks, depending on the study site and time of year. Additionally, some lake time series were only sampled once during the winter, depending on ice conditions. Empirically, observation frequency did not have an influence on tailedness (Supporting Information Fig. S1). Northern Wisconsin lakes include Allequash Lake, Big Muskellunge Lake, Crystal Bog, Crystal Lake, Little Rock Lake, Sparkling Lake, Trout Bog, and Trout Lake. The only time series from Little Rock Lake was ice cover. Southern Wisconsin lakes include Fish Lake, Lake Mendota, Lake Monona, and Lake Wingra.
Meteorological variables included air temperature (daily minimum, maximum, range), sum daily precipitation, snow depth, and daily average wind speed. Physical variables included dissolved oxygen concentration, dissolved oxygen percent saturation, water temperature, the fraction of light at depth relative to light at surface, light extinction coefficient, Secchi depth, lake level, duration of the ice-free season, and the depth of the epilimnion. Chemical variables included alkalinity, bicarbonate-reactive silica (unfiltered), calcium, chloride, conductivity, dissolved inorganic carbon, dissolved organic carbon, dissolved reactive silica (filtered), iron, potassium, magnesium, manganese, sodium, ammonium, nitrate plus nitrite, pH, sulfate, total inorganic carbon, total organic carbon, total nitrogen (unfiltered), total phosphorus (unfiltered), and total particulate matter. Biological variables included the abundance of fishes and zooplankton, and epilimnetic chlorophyll concentration. Each genus of fish or zooplankton comprised its own time series.

Distribution fitting and waiting times
We fit a generalized extreme value distribution (GEV) to the time series of annual maxima for each variable; this is the theoretical limiting distribution of maxima that is approached for many statistical processes with large enough sample sizes, much like the normal distribution is the limiting distribution of sums of random variables (Coles 2001). The cumulative distribution function of the GEV is where the case of n 6 ¼ 0 is subject to the constraint that 11nðx2lÞ=r > 0. Variables that follow very different distributions can nonetheless have their tails compared by the GEV because it is the limiting distribution for maxima. The extremes to which we fit the GEV were defined via the block-maximum method (Palutikof et al. 1999) with blocks corresponding to years to yield annual maxima. Each time series contained between 15 and 158 annual maxima (median 5 29) and was categorized as biological (n 5 283 lake time series, including chlorophyll, 28 fish genera, and 34 zooplankton genera), chemical (n 5 220, e.g., ion concentration), physical (n 5 80, e.g., water temperature), or meteorological (12 regional time series, e.g., mean daily wind speed). We used the GEV to assess tailedness and estimate waiting time between record-breaking extremes. The GEV parameters were fit using maximum likelihood after removing a linear temporal trend from time series of maxima (Katz et al. 2005). We calculated the unconditional expected waiting time as the inverse of the probability p (calculated from the cumulative distribution function) of observing a value 10% over the current record, multiplied by the observation frequency (F, yr 21 ) of maxima in the time series (F was usually 1); thus, waiting time 5 1/p*F. This calculation does not account for possible autocorrelation in the annual maxima in the time series; if there is autocorrelation, then the expected conditional waiting time will depend on the most recent value of the time series. Statistical analyses were performed in the statistical programming language R (R Core Team 2014). In some cases (79 of the 595 time series; 0 biological, 44 chemical, 32 physical, 3 meteorological), waiting times were undefined because the value 10% over the record (x) was outside the support of the generalized extreme value distribution (GEV): x 2 ½l2r=n; 11Þ when n > 0; x 2 ð21; 11Þ when n50; x 2 ð21; l2r=n when n < 0 For example, the return level could be outside the support if the estimate of n is less than zero, and the return level is greater than the limit of the upper tail (l-r/n). Undefined waiting times were excluded from comparisons of waiting time between variable types. Note that removing the 79 time series with undefined (infinite) waiting times makes Fig. 2 conservative with respect to showing that biological time series have shorter waiting times than the other categories.

ARMA models
We used autoregressive moving average (ARMA) models to assess the relative contributions of processes occurring at time scales shorter (within-sample) or longer (between-sample) than a sampling interval (Supporting Information: Using ARMA to assess the time scale of sources of fat tails). ARMA models are a generic class of time series model capable of approximating many processes by predicting the current value of a variable from its past values and with past values of error terms. The parameter values of an ARMA model can be used to make general inferences about a variable's temporal dynamics, such as if it has long-term "memory," which is the autocorrelation of current values with values from many time steps ago. Furthermore, the residuals from the fitted ARMA models give information about the dynamics after memory has been removed. We fit ARMA models to parent time series, with the number of AR (p) and MA (q) parameters being chosen by AIC, subject to the constraint that 1 p 3 and q < p (Supporting Information: ARMA models). If the set of block maxima for a time series is fat-tailed and that for the ARMA residuals is not, this indicates that between-sample dynamics are responsible for generating the fat tail. Conversely, if both the time series and the ARMA residuals are fat-tailed, then within-sample dynamics are responsible at least in part for the generating the fat tail (Ives et al. 2003. ARMA is only an approximation of these between-sample processes, and could fail in some cases. Nonetheless, this approach can be informative and provides a common analysis for comparisons across our large set of diverse time series.

Frequency of extremes
We found a wide range of tailedness (n) in each variable category, but on average the biological variables had higher values of n than other categories of variables ( Fig. 2A; Supporting Information Table S1), even after controlling for length of time series and the identity of the lake (pairwise comparisons with biological time series: chemical, z 5 210.1, p 0.0001; physical, z 5 211.7, p 0.0001; meteorological, z 5 22.1, p 5 0.0234 (Supporting Information Table S2). For waiting times to break a current record, the pattern across categories was similar to that for n: biological variables were likely to break their current records sooner on average than other categories of variables (Fig. 2B). These results indicate that the biological variables tended to produce extremes more frequently than other categories of ecosystem variables.
To further compare the different categories while accounting for uncertainty in the estimates of n, we tested whether the time series were fat-tailed by, for each time series, comparing the null hypothesis that n 5 0 against the alternative n > 0 for positive estimates or n < 0 for negative estimates (1-tailed test with a 5 0.05, not correcting for multiple comparisons). We present results that are not corrected for multiple comparisons, and subsequently those with p-values that are corrected to maintain the false discovery rate (Benjamini and Hochberg 1995). Of the 283 biological time series, 38% Bio (283) Chem (220) Phys (80) Met (12) Tailedness (ξ) had an estimate of n statistically significantly greater than 0 (corrected 5 12%), by far the highest frequency of statistically significant fat-tailedness among categories (13% for chemical, 3.8% for physical, 0% for meteorological; corrected 5 1.8%, 0%, and 0%, respectively). Using the same procedure to investigate distributions with bounded upper limits (n < 0), only 0.7% of biological time series had statistically significant estimates of n < 0 (corrected 5 0.4%), which was the lowest frequency among variable categories (22% for chemical, 51% for physical, 42% for meteorological; corrected 5 8.2%, 25%, 25%, respectively). Among the biological variables, the fish, zooplankton and chlorophyll time series with the largest estimates of n were, respectively, Lepomis spp. in Lake Mendota (n 5 1.2; 90% CI: 0.7 < n < 1.8), Kellicottia spp. in Crystal Bog (n 5 1.2; 90% CI: 0.7 < n < 1.7), and chlorophyll in Crystal Lake (n 5 0.5; 90% CI: 0.2 < n < 0.8).

ARMA dynamics
If extreme events were generated largely by long-term dynamics created by correlations across multiple time steps in the parent time series, then estimates of n for the original time series should be uncorrelated with estimates of n for the ARMA residuals. Additionally, the magnitude of correlation across time steps, measured by the eigenvalue of the ARMA process, should be positively correlated with the tailedness. Both of these expectations were false. Estimates of n for the time series were positively related with the estimates of n calculated for the ARMA residuals (slope 5 0.40, p 0.0001; Supporting Information Fig. S2, Fig. 3). Moreover, estimates of n for the times series were negatively related with the eigenvalues (slope 5 20.23, p 0.0001; Fig. 4). These correlations are consistent with the inference that extreme events were generated primarily by processes that occur within a time step.

Discussion
Biological time series had distributions with fatter tails than time series of other ecosystem variables. As tailedness (n) increases, the tail probabilities become large or "fat" relative to distributions like the normal, implying that extreme events will be more common. A fatter tail also decreases the time that passes until an event of arbitrarily large size is observed, increases the largest value observed after an arbitrary amount of time, and makes the chances of observing "massive" and "big" events more similar. Thus, in our analyses biological variables have larger and more frequent extremes. As values of n increase for a distribution, the statistical moments of the distribution 1/n become undefined (Coles 2001); for example, at n 1 the mean (first moment) is infinite, and at n 0.5 the variance (second moment) is infinite. By contrast, biological variables are often characterized by distributions that have defined means and variances (e.g., lognormal). However, the estimates of n for some of the fish and zooplankton time series were greater than 1, implying that their annual maxima will not converge to a long-term mean as more data are collected. Therefore, even distributions like the lognormal would underestimate the frequency and intensity of extremes in these fat-tailed biological variables. This result underscores the challenge that fat-tailed distributions pose for management: prevalent extreme events can make it impossible statistically to characterize even the mean abundance of a species.
Compared to meteorological, physical, and chemical time series, the biological time series we analyzed had larger values for tailedness (n) and shorter expected waiting times until record-breaking events, even though the variables were measured in the same ecosystems. This result provides evidence to suggest that biological variables should be expected to produce larger and more frequent surprises than the other categories of variables. Overall, these analyses provide evidence that biological variables are more likely to be fat-tailed than other types of variables in the same ecosystems.

Mechanisms
Mechanisms for promoting or diminishing fat tails are common in biological processes (Sornette 2006). Negative density dependence can lead to thinner tails and is common in population dynamics (May 1976;Ziebarth et al. 2010), giving reason to expect less fat-tailedness in biological variables relative to the environmental variables that drive them. However, nonlinear models fit to time series of disease (Sugihara et al. 1990), flour beetles (Dennis et al. 2001), spruce budworms (Ludwig et al. 1978), and phytoplankton (Beninc a et al. 2008) can produce fat tails under certain parameter combinations, even deterministically. Nonlinear processes are common in biology (Hsieh et al. 2005), suggesting that they may produce fat tails in biology just as they are known to do in other settings (Barab asi 2005;Newman 2005).
Our ARMA analysis suggests that mechanisms within a time step ( 1 year) often generate observed fat tails, and that increasingly long memory cannot explain the observed fat tails. A simple process for generating such fat tails is for a population to amplify the variability of its environment. Multiplicative birth and survival processes influenced by stochastic environmental drivers are one way that fat tails could be generated within time steps ( Fig. 1 and Supporting Information: Extremes in a model population). Comparing alternative mechanisms using models tailored for specific biological systems is needed if we are to better understand how fat tails are generated.

Concluding remarks
Our analysis included all of the variables typically collected by limnologists (Magnuson et al. 2006), and because our results might have differed if we had investigated different subsets of the variables, we analyzed all of the data that were available in a comprehensive analysis. Estimating a statistic based on rare events requires long time series (Palutikof et al. 1999;Katz et al. 2005), and estimates of n will improve with continued data collection. In fact, this need for long time series may serve as an explanation for the relative paucity of empirical assessments of biological tailedness-time series like those used in this study have only recently matured to a length suitable for extreme value analysis, whereas long environmental time series are more common. Therefore, continuing long-term monitoring of ecological variables and characterizing the tailedness of additional time series are critical for understanding ecological extremes.
Because fat tails dictate the frequency and intensity of massive events produced by an ecosystem, biological fat tails stand to greatly impact the management of populations that are critical to human well-being, such as fish stocks, invasive species, pests and diseases, or phytoplankton that form harmful blooms. Although there are several classes of general causal mechanisms for fat tails (literature cited above), more research is needed to identify causes of fat tails for specific populations. Our study focused on 11 lakes, but similar analyses of additional ecosystems are needed to test the generality of the patterns found here. Tests in additional ecosystems are also crucial for determining whether certain taxa consistently have higher tailedness, or if certain ecosystem characteristics can affect the tailedness of its denizens. Such patterns could help make predictions of extremes for ecosystems without long-term records, and would provide further insight into the causes of biological tailedness, its impacts on ecosystem functioning, and its implications for human use and management of ecosystems.
It is dangerous to consider the future as a set of norms from the past (Oppenheimer et al. 2007;Burgman et al. 2012), and forecasting future events, especially extreme events, is difficult. Our empirical results indicate that biological time series present frequent surprises. If this pattern was caused by directional environmental change in climate, nutrient inputs, habitat loss, or other factors, then extremes may become more common in ecosystem time series.