Interannual variability is a pervasive feature of aquatic ecosystems. This variability results from short- and long-term dynamics of biotic and abiotic origin, inclusive of multiannual variability and long-term trends. Although understanding short-term variability and forecasting directional change are important research efforts, far less attention has been paid to oscillatory, or wave-like dynamics that play out over multiple years, in aquatic ecosystems. In this essay, we argue that understanding these modes of variability—in addition to directional trends and intraannual patterns—and their underlying causes are necessary for understanding aquatic ecosystem functioning over long time periods for effective conservation and management. Fortunately, given the growing availability of multidecadal data, development of statistical tools, and the urgent need to forecast change, the field can readily adopt multiannual dynamic thinking into our understanding of aquatic ecosystems. Environmental change occurs over a broad range of timescales. Aquatic ecosystems can change rapidly from disturbances that drastically affect structure and function. Other changes progress more slowly, due to processes such as climate change, eutrophication, changes in watershed land use and flow regime, biodiversity loss, and biological invasions. These long-term drivers tend to cause directional, slow change or the abrupt crossing of thresholds, leading to temporal trends or regime shifts. However, other processes operating at long timescales drive variability to ecosystem structure and function without necessarily resulting in directional change. Multiannual dynamics, or wavy, periodic, and quasiperiodic oscillations operating over timescales from 2 years to over a decade, are often a substantial source of variability that can be independent of long-term trends. Although multiannual dynamics are often treated as "operating in the background," drivers of oscillations and trends operate at all timescales, in some cases individually et al synergistically, to regulate the structure and function of aquatic ecosystems. These multiannual dynamics have been shown to be important in long-term studies such as the eutrophication and recovery of Lake Washington (Hampton et al. 2006), the effect of climate oscillations on calanoid copepods (Fromentin and Planque 1996), overexploitation as in the northwestern Atlantic cod collapse (Hutchings and Myers 1994), and species invasion as in the effects zebra mussels on the Hudson River (Strayer et al. 2014). In some ecosystems, there may be complete absence of trends in a variable of interest over multiannual timescales, but not an absence of pattern. For example, monthly mean nitrate concentrations for the past 40 years in the Des Moines River have no discernable trend despite being quite variable, ranging from below detection to greater than 18 mg L−1 (Fig. 1A). The lack of a trend is somewhat surprising given the history of land use change and agricultural intensification in the region during this time period (Yu and Lu 2018). However, while there are no long-term trends in nitrate concentrations in the river, there are strong oscillatory patterns in the time series. A wavelet analysis of the multiannual dynamics of the nitrate time series reveals that there are repeating oscillations at seasonal, annual, 3–5 year, and 10–14 year timescales (Fig. 1B) (see Supporting Information Appendix S1 for method details). This example illustrates the rich information that can be gleaned from multiannual pattern analysis of long-term data. Interest in characterizing multiannual dynamics is not new, per se. Oceanographers have long appreciated multiyear dynamics in currents, ocean–atmosphere connections, and impacts of this variation on distributions and populations of marine fauna (Di Lorenzo et al. 2013; Tommasi et al. 2017). Lotic ecologists have analyzed variations in stream and river discharge to assess directional, extreme, and periodic changes on long timescales (Palmer and Ruhi 2019). Similarly, lake researchers have conducted long-term studies of ecological variation in the context of the interactions of external drivers with internal processes (Hampton et al. 2006; Keitt and Fischer 2006). Yet, we argue that as research in limnology and oceanography is increasingly marked by long-term, intensive ecosystem monitoring, cross-disciplinary research, big data, and open science, the discipline is well positioned to begin routinely incorporating multiannual dynamics into our analysis of changes in aquatic ecosystems. Additionally, the increasing ubiquity of environmental change due to multiple anthropogenic stresses playing out over varying spatial and temporal scales necessitates disentangling multiannual trends and oscillations in order to effectively manage and conserve aquatic ecosystems long term. To that end, in this essay, we discuss the drivers that lead to multiannual variability, the consequences of multiannual variability on ecosystem functioning, and suggest strategies for characterizing and incorporating multiannual variability into aquatic ecosystem research, conservation, and management. Drivers of variability in ecological processes operate at time scales from days to decades. For many ecological processes, temporal patterns emerge at multiannual timescales, such as the population dynamics of long-lived organisms. While multiannual dynamics may interact with, and can even depend in part, on processes occurring at shorter time scales, they often cannot be observed, predicted, or understood from short timescale dynamics alone. If fact, multiannual oscillations can greatly affect the conclusions drawn about the trajectory of an ecosystem depending on the length and position of the observation window (e.g., colored panels in Fig. 2A). In the Des Moines River example, strong decadal oscillations in nitrate concentrations produced apparent shorter-term trends (e.g., decreasing concentrations from c.a. 1981–1990; red line in Fig. 1B) that could be erroneously attributed to directional changes brought on by improvements in watershed nutrient management. Multiannual dynamics manifest across multiple timescales, overlaying and interacting with each other. As an illustration of the consequences of superimposed dynamics in aquatic ecosystems, consider the long-term variation of primary productivity. Ultimately, primary productivity is driven by temperature, light, mixing, nutrient availability and grazing. The short-term diel and seasonal oscillations of light availability, temperature, and nutrient limitation are well understood in aquatic ecosystems (Fig. 2B). However, some of these same drivers also vary at longer timescales. Interannual variability in temperature, nutrient loading and deposition, precipitation, residence time, and grazer populations can all be driven by climate oscillations (Paerl et al. 2015; Carey et al. 2016; Nergui et al. 2016) (Fig. 2C). These multiannual oscillations vary in their period from a few years to decades. Finally, these oscillations can also interact with directional changes such as eutrophication, which stimulates primary productivity over time (Fig. 2D). All of these dynamics are superimposed resulting in a highly variable time series of measured primary productivity (Fig. 2A). The numerous phenomena that generate multiannual patterns manifesting in population dynamics, community composition and structure, and aquatic ecosystem function can be placed into two broad categories: intrinsic cycles and periodic oscillations driven by external forcing. Internal dynamics can either lead to the rise of chaos or intrinsic cycles, such as predator–prey dynamics, in the system, and will tend to have high amplitudes and consistent cycle lengths, even in the absence of a periodic external driver. For example, density-dependent reproduction and stage structure prominently contribute to population cycles (Myers 2018). Systems can also exhibit periodic dynamics that arise from external drivers that also fluctuate periodically. The El Niño-Southern Oscillation (ENSO), Pacific Decadal Oscillation, and other climate teleconnections are all examples of external quasiperiodic dynamics that operate at multiannual timescales. Climate oscillations are tied to periodic flood regimes (Palmer and Ruhi 2019), shifts in oceanic currents, upper layer mixing dynamics, patterns of extreme weather, and continental scale patterns of nutrient deposition (Nergui et al. 2016), which results in year-to-year variability in environmental parameters. In addition to linear trends, other forms of nonlinear directional change such as hysteresis and time lags also contribute to long-term variability. For example, legacy accumulation of nutrients in watersheds and aquatic sediments often results in a time lag between the reduction in nutrient inputs and an improvement in water quality (Van Meter and Basu 2017; Kusmer et al. 2019). Likewise, the trajectories of estuarine ecosystems recovering from eutrophication follow different paths than when nutrients were increasing, and typically do not return to the same starting point (Duarte et al. 2009). All of these drivers of variability—cycles, periodic dynamics, and trends—are interacting to shape the dynamics in aquatic ecosystems. What external forces or intrinsic dynamics may be driving the strong oscillations in nitrate concentration in the Des Moines River? Variability in river flow is a likely driver of the variability in nitrate concentrations as baseflow contributes substantially to nitrate flux. In this agriculturally rich river network, tile drainage contributes half of the annual water yield in many of the subbasins (Schilling et al. 2019). We can evaluate this hypothesis by comparing the decomposed dynamics, or wavelets, of driver (flow) and response variable (nitrate) time series. There was high wavelet power at similar timescales in both nitrate concentration and river flow in the Des Moines River (Fig. 1C,D, respectively), suggesting a relationship between these variables' dynamics. We can more specifically test this hypothesis by evaluating the wavelet coherence, which measures whether two time series have correlated magnitudes of oscillation and consistent phase differences (time lags) throughout a time series, as a function of timescale (see Supporting Information Appendix S1). The wavelet coherence analysis revealed that the dynamics were significantly related at annual (p < 0.001), 3–5 year timescales (p < 0.001), and 10–14 year (p = 0.036) timescales. We can further hypothesize what drivers control variability in flow at these timescales and could use wavelet coherence analyses to evaluate support for these predictions. For example, at particular timescales nitrate concentration and flow may be driven in part by multiannual climate oscillations such as ENSO (Jones et al. 2012). The incentive to examine multiannual dynamics is growing as the need to scale across time, and space is increasingly important for predicting and coping with large-scale environmental changes (Payne et al. 2017; Dietze et al. 2018). In addition to changing climate means and extremes, the autocorrelation structure of climate is changing (Lenton et al. 2017). Understanding the resulting ecological change and grounding long-term forecasting requires understanding patterns and drivers of multiannual oscillations in order to accurately capture variability that is independent of, or synergistic with, long-term directional changes. Fundamentally, we cannot contextualize or predict future dynamics if we do not understand the role of multiannual dynamics as a driver of variability in aquatic ecosystem structure and function. Long-term management and conservation of aquatic ecosystem services will be strengthened by gains in understanding long-term dynamics. These benefits are evident in the example of the Des Moines River nitrate dynamics. The river supplies drinking water for the City of Des Moines (Iowa). Due to human health concerns, the water utility must periodically remove nitrate from the source water prior to distribution, which is an expensive process. Analyzing multiannual dynamics revealed previously uncharacterized, but predictable patterns in nitrate concentration peaks and troughs at 3–5 and 10–14 year timescales (Fig. 1B) in addition to the well-known seasonal oscillations that are already incorporated into the water utility's management plan. In this instance, identifying the multiannual oscillations provides key information for long-term water source and budget management for the water utility. Many of the environmental threats to aquatic ecosystems and the services they provide are influenced by multiannual dynamics. The development of marine harmful algal blooms is influenced by the local response to climate oscillations, which can be modified by the co-occurrence of oscillations such as ENSO and Pacific Decadal Oscillation (Moore et al. 2008). The variability in ice cover duration of lakes and rivers in the northern hemisphere is also driven by high frequency climate oscillations (Schmidt et al. 2019) in combination with directional shifts to earlier ice-off dates due to climate change (Sharma et al. 2016). The changes in ice cover duration due to both oscillations and climate change influence ecosystem function in both the winter and following ice-free seasons (Hampton et al. 2017). In lotic ecosystems, species with life histories driven by changes in flow regimes necessitate the identification of the dominant frequencies of variation at both subannual and multiannual scales in order to effectively manage and conserve these organisms (Palmer and Ruhi 2019). These examples illustrate the important role that multiannual oscillations play in dictating the interannual variability in aquatic ecosystems. Detecting multiannual dynamics requires long-term data. The benefits of long-term studies within limnology and oceanography are well recognized and can serve as a foundation for further improving the understanding of multiannual ecosystem dynamics (Hampton et al. 2019). Robust investigation of wavy patterns requires time series containing multiples instances of the pattern of interest. Given the period lengths of many climate oscillations and other ecological cycles, this often means multidecadal time series. Even longer time series are needed to identify nonstationary patterns such as changes in the period or amplitude of oscillations. Despite the potential insights provided by analysis of long-term data, long-term data collections comprise only approximately 5% of limnological and oceanographic research published in the past four decades (Xenopoulos 2019). For the existing long-term data sets, there are often complications due to observation error generated from methodological imprecision, changes in sampling methods, detection limits, and shifts in personnel over time. Managing observation error requires statistical techniques that can account for both observation and process error and still allow patterns and trends to be detected (Hampton et al. 2019).Yet, given the ubiquity of multiannual dynamics in aquatic systems, where data are available we need statistical tools adept at decomposing and detecting these important patterns. This can be difficult given that multiannual dynamics often manifest as superimposed trends and multiple wave-like components having specific timescales of variability. Fortunately, a number of possible analysis tools exist, and we highlight a handful of these approaches that are notable due to their longstanding use in ecology or their growing popularity. Ultimately, the choice of statistical method for exploring multiannual dynamics should be dictated by the research question. Autocorrelation functions have been long available and a widely used tool to identify periodic oscillations in time series. However, it can be difficult to detect when many oscillations are present with differing periods (i.e., multiple periodicities), such as in the Des Moines River nitrate example (Fig. 1). Additionally, irregular oscillators such as ENSO may not produce clear and consistent peaks in the autocorrelation function. An alternative is autoregressive integrated moving average modeling which can fit complex patterns arising from both density-dependent processes and environmental forcing (Ives et al. 2010), but the meaning of higher-order coefficients is often opaque. Applying windowing or temporal filters to time series can amplify patterns at certain timescales while diminishing others, but is also prone to producing spurious, yet apparent cycles as an artifact. Dynamic factor analysis (Zuur et al. 2003), which is a state-space modeling approach, can identify the common patterns underlying the dynamics in a set of time series, such as shared trends and cycles, and can reveal multiannual wavy patterns. An additional benefit of state-space modeling is the treatment of missing data and observation error, which can be common issues in ecological time series. Finally, spectral techniques based on Fourier and wavelet analysis have great strength for analyzing periodic dynamics, even when multiple periodicities are present. In particular, wavelet analysis reveals timescale (the inverse of frequency) specific patterns and changes through time in system behavior, and there are robust tests for time- and timescale-dependent relationships among variables of interest and drivers. However, while spectral techniques are useful for analyzing periodic dynamics, they are not well suited to trend analysis and require evenly spaced, continuous time series for analysis. Fortunately, we are well positioned as a discipline to incorporate multiannual dynamic thinking into our understanding of aquatic ecosystem function and change. The continued growth and acceptance of open science practices in our field, coupled with a mounting collection of long-term data sets is creating the opportunity to ask multiannual questions in many systems. Additionally, the continued development of statistical methods that disentangle complex, nonstationary, and interacting signals is providing tools necessary for addressing these questions. As the environment continues to rapidly transform in response to numerous interacting stressors, incorporating multiannual dynamic thinking into our understanding of aquatic ecosystems will help us meet the challenges of a changing planet. Appendix S1: Method details Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article. This research was supported by the National Science Foundation Office of Advanced Infrastructure award #1839011 to G.M.W. and #1839024 to MLP. J.A.W. was also supported by a NatureNet Science Fellowship. We thank William Crumpton and Greg Stenback for providing the compiled Des Moines River nitrate and flow time series.
