Slowing down of North Pacific climate variability and its implications for abrupt ecosystem change

• regime shifts
• North Pacific
• marine ecosystems
• abrupt change
• climate change

Sea surface temperature (SST) fluctuations are well known to exhibit “red” spectra––with increased power at lower frequencies––even when forced by purely “white” noise from the atmosphere (1, 2). Individual realizations of a stationary red noise process will typically drift one side of the mean and later switch over, resembling an irregular oscillation. One such “oscillation,” the Pacific Decadal Oscillation (PDO) index (3, 4), is described by the variation of the dominant spatial pattern of SSTs >20°N in the Pacific over time, having removed the seasonal cycle and an overall warming trend. Shifts in the PDO between its positive and negative phases––especially around 1977 and 1989––have been linked to abrupt changes in salmon productivity (3, 5), drought regimes in the United States (6), changes in Indian summer monsoon rainfall (7), and a host of other ecological and climatic time series (4). Whereas some early work described the PDO itself as exhibiting regime shifts (4, 8), implying the existence of alternative climate attractors, subsequent work has found no evidence for nonlinearity in the PDO or other North Pacific climate indices (9). Instead, North Pacific climate variability can be characterized as a linear, stochastic, red noise process (10). North Pacific marine ecosystems have also been widely described as exhibiting regime shifts between alternative attractors in response to stochastic fluctuations in the physical climate system (9). Whereas there is evidence for nonlinearity in some North Pacific ecosystem time series (9, 11), other populations appear to linearly track the climate forcing, especially if their generation time matches the damping timescale of the forcing (11). Thus, one “double-integration” model for the overall behavior is that SSTs (and indices derived from them such as the PDO) integrate white noise forcing from the atmosphere to produce red noise (1), and ecosystems further integrate this red noise to create time series that are even redder (12). The redder the ocean climate variability, the more closely a marine ecosystem will track it (13), making the ecosystem more prone to large changes driven by the climate variability (13).

Here we consider whether North Pacific climate variability has changed over the observational record. A simple generic model for SST fluctuations (10) is the first-order linear autoregressive process:xt=αxt-1+σηt[1]where the subscript t denotes time, x is the time series, α is the lag-1 autocorrelation [AR(1)] coefficient (0 < α < 1, for red noise), and η_t is Gaussian white noise of amplitude σ. Existing work has used a constant value of α = 0.95 (corresponding to a decorrelation timescale of 20 mo) to mimic the PDO index (10). If, however, there has been a change in the spectrum of SST variability we expect this to manifest as a change in α, with an accompanying change in variance (14, 15). Hence, for each of a series of North Pacific datasets [the PDO index (3), HadISST (16), ERSST v3 (17), and HadSST3 (18) detailed below], we estimated the AR(1) coefficient (α) and calculated the variance, in a sliding window moved through the dataset, to look for any trends (Methods). In each case, the results were tested against a null model for North Pacific SST fluctuations, with fixed α = 0.95 (10), which was used to generate 10,000 realizations of a series of the same length as the time series being tested (Methods)––the idea being that a relatively short realization of such a stationary red noise process can display a wide range in trends of lag-1 autocorrelation and variance.

Slowing down of fluctuations is visible in the original PDO index (Fig. 1A), with a strongly increasing trend in the AR(1) coefficient (Kendall τ = 0.941, Methods), from an estimated α = 0.57 to α = 0.95 over the 112.5-y record (Fig. 1B, Methods). Variance also generally increases (Fig. 1C), but not as strongly across the whole record (τ = 0.599), due to a downturn in recent decades. The observed trend in AR(1) in the PDO index lies outside the distribution resulting from 10,000 realizations of the null model (Fig. 1D). Hence we can reject with high confidence (P < 0.0001) the null hypothesis that the observed trend in autocorrelation in the PDO index is the result of a red noise process with fixed α. The increasing trend in variance in the PDO index (τ = 0.599) is less significant and could occur by chance with P = 0.1411 (Fig. S1A). However, other aspects of North Pacific climate are showing increases in variance (19), and if we just consider the data up to 1989 the PDO index shows a strong increase in variance (τ = 0.909), noting that the decrease in variance appears to begin in 1960 in Fig. 1C because the indicators are plotted at the midpoint of the window used to calculate them and at this point data from 1989 enter the sliding window.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Slowing down observed in the PDO index, 1900–2012. (A) The original PDO index. (B) Estimated increase in AR(1) coefficient, using a window length of 675 points (half the series), without detrending (Methods), results plotted in the middle of the sliding window. (C) Estimated increase in variance. Trends in the indicators are expressed as Kendall τ-values. (D) Range in τ-values expected from ensembles of 10,000 realizations of a null model with fixed α = 0.95. Red vertical line denotes the τ-value found in the PDO index A (P < 0.0001). (E and F) Sensitivity analyses of the AR(1) coefficient estimate and variance respectively by testing the value of τ for a variety of window lengths and bandwidths.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S1.

Assessing the significance of trends observed in (A, C, and E) variance and (B, D, and F) variance at high frequencies (having used a low bandwidth for detrending the original data), for (A and B) Mantua PDO index, (C and D) HadISST average North Pacific index, and (E and F) HadSST3 average North Pacific index. The histogram in each case represents 10,000 runs of a null model (described in the main text), the red line is the result from analysis of the corresponding index, and P values for each hypothesis test are shown above each plot.

Nonstationary behavior can cause autocorrelation to increase; hence we examined the effect of further detrending the data before analysis. As the bandwidth is decreased, the lowest frequencies, including the (multi-) decadal “oscillation” itself, are the first to be removed and the results are limited to ever shorter timescale fluctuations. As a further sensitivity analysis, we also varied the sliding window length in which the indicators are calculated. The positive trend in AR(1) is robust to varying sliding window length and filtering bandwidth (Kendall τ = 0.77–0.98, Fig. 1E). The positive trend in variance is also fairly robust (Fig. 1F), except when using a very short filtering bandwidth for detrending. This leaves only the highest frequency variability in the index and its variance decreases (e.g., τ = −0.903 with the shortest filtering bandwidth). This is consistent with a shift in power from high to low frequencies, which can also be seen in the changing power spectrum of the data (Fig. S2A). Comparing to the null model (also filtered with the shortest bandwidth before analysis), the decline in variance at high frequencies has P = 0.0006 and we can reject the null model at 5% significance (Fig. S1B).

• Download figure
• Open in new tab
• Download powerpoint

Fig. S2.

Power spectra of SST variability (log–log plots) for (A) Mantua PDO index, (B) HadISST average North Pacific index, and (C) HadSST3 average North Pacific index. In each case, power spectra for the first half of the data are shown in gray and the second half in orange. A regression line has been fitted through each of these (black and red, respectively), showing a shift in power from high to low frequencies in all three time series. There is no overlap between the 95% confidence intervals on the gradient of the regression lines when using the PDO index (A), an overlap of ∼4% the width of the 95% confidence intervals for HadISST (B), and an overlap of ∼75% of the width of the confidence intervals for HadSST3 (C).

To establish whether slowing down has anything to do with the spatial pattern of the PDO we examined the HadISST (16) dataset. Removing seasonal and overall warming trends and then simply averaging HadISST over the North Pacific domain 20–60°N (Methods), slowing down is again visible by eye (Fig. S3A) and confirmed by a strongly increasing trend in the AR(1) coefficient (τ = 0.938, Fig. S3B) with estimated change in α from 0.68 to 0.92 (the 95% confidence interval on the change in α hardly alters these results, given to two significant figures). There is a weaker increasing trend in variance (τ = 0.275, Fig. S3C). For the trend in AR(1), a null model of fixed α = 0.95 can be rejected at P < 0.0001 (Fig. S3D), whereas the positive trend in variance is not significant (P = 0.3518, Fig. S1C), again due to a decline in variance in recent decades. The positive trend in the AR(1) coefficient is robust to varying detrending bandwidth and sliding window length (τ = 0.70–0.95, Fig. S3E). Variance again decreases at the highest frequencies (τ = −0.892 for the shortest filtering bandwidth, Fig. S3F), with P = 0.0007 against the null model (Fig. S1D), consistent with a shift in power from high to low frequencies that is seen in the power spectrum (Fig. S2B). Thus, slowing down of North Pacific SST fluctuations is not particularly associated with the spatial pattern of PDO variability and has its own spatial pattern (Fig. S4).

• Download figure
• Open in new tab
• Download powerpoint

Fig. S3.

Slowing down in average North Pacific reconstructed SSTs,1870–2011. (A) HadISST data detrended and averaged over the North Pacific (Methods). (B) Estimated increase in AR(1) coefficient, using a window length of 850 points (half the series), without detrending (Methods), results plotted in the middle of the sliding window. (C) Estimated increase in variance. (D) Range in Kendall τ-values expected from ensembles of 10,000 realizations of a null model with fixed α = 0.95, with red vertical line denoting τ = 0.938 found in A (P < 0.0001). Sensitivity analyses of (E) the AR(1) coefficient estimate and (F) variance by testing the value of τ for a variety of window lengths and bandwidths.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S4.

Spatial pattern of slowing down in North Pacific SSTs in HadISST data when the mean annual cycle and a quadratic warming trend have been removed from each grid point. (A) Kendall τ-values are shown as well as (B) the change in α over the data set using the start and end points of a linear regression fitted to each AR(1) time series. The same analysis is carried out on HadSST3 to create Fig. 3 A and B.

We also analyzed the ERSST v3 dataset (17) and the results for average North Pacific SSTs also show a positive trend in AR(1) over time, but there are considerable differences in the early part of the two datasets and hence the corresponding AR(1) estimates (Fig. S5). Data sampling was generally sparser in the past; hence, datasets are subject to more infilling further back in time, which in turn could affect measures of autocorrelation. Hence, we analyzed the original HadSST3 dataset (18) without infilling, which has a coarser spatial resolution than HadISST. We focus on the interval 1950 onward (Fig. 2A) because data collection in the North Pacific was sparse before that. Despite the much shorter time interval, we find increasing AR(1) (τ = 0.832, Fig. 2B) with α increasing from 0.72 to 0.93, accompanied by an increase in variance (τ = 0.816, Fig. 2C). The increase in AR(1) has P = 0.0042 against the null model (Fig. 2D) and the trend in variance has P = 0.028 (Fig. S1E). Increasing trends in AR(1) and variance are generally robust to changes in sliding window length and detrending bandwidth (Fig. 2 E and F) with the exception that variance declines for bandwidth ∼200 and window length >450. A further slight decrease in variance at high frequencies (τ = −0.171 for the shortest filtering bandwidth, Fig. 2F) is not significant (P = 0.4093, Fig. S1F), commensurate with only a slight shift of power to lower frequencies (Fig. S2C). Nevertheless, slowing down of North Pacific SST fluctuations has occurred just since 1950.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Slowing down in average North Pacific SST raw data, 1950–2006. (A) HadSST3 data detrended and averaged over the North Pacific (Methods). (B) Estimated increase in AR(1) coefficient, using a window length of 350 points (half the series), without detrending (Methods), results plotted in the middle of the sliding window. (C) Estimated increase in variance. (D) Range in Kendall τ-values expected from ensembles of 10,000 realizations of a null model with fixed α = 0.95, with red vertical line denoting τ = 0.832 found in A (P = 0.0042). Sensitivity analyses of (E) the AR(1) coefficient estimate and (F) variance by testing the value of τ for a variety of window lengths and bandwidths.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S5.

Time series of average North Pacific SSTs from (A) HadISST and (B) ERSST v3 (1900 onward, with the average annual cycle and quadratic warming trend removed). (C) Estimates of the AR(1) coefficient for each time series using a window length of 675 points (approximately half the time series), plotted at the middle of the time window it is calculated on. A regression line is fitted through each of these AR(1) estimates (dotted line), and Kendall τ-values for the trends in AR(1) are given.

A grid-point by grid-point analysis of HadSST3 data was conducted to determine where the slowing down of SST fluctuations is occurring, revealing that it is widespread but not universal across the North Pacific domain (Fig. 3 A and B). Strong slowing down occurs around the basin edge, e.g., off the West Coast of North America. To try and explain this slowing down we consider a simple physical mechanism consistent with known climatic trends, namely deepening of the mixed layer (20). In the original model of Frankignoul and Hasselmann (2) the key environmental variables affecting α are mixed layer depth (h) and average wind speed (U):α=1−kU/h[2]We estimate the constant, k (s mo⁻¹), for each grid point using the mean values of U [from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis (21)], h from observations (20)], and α (from the midpoint of the fitted trend, Fig. 3B) (Methods). Trends in reanalysis (21) wind speed (U) across the domain are mixed and generally small so for simplicity we hold U constant, and consider the change in h required in Eq. 2 to explain observed trends in α. We compare this to observed changes in mixed layer depth (MLD). Due to limited availability of MLD data we start the analysis in 1960. Due to the large variability in MLD, we use linear regression to determine the overall trend in MLD and use the start and end points of the regression line to estimate the change.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Spatial analysis of slowing down in SSTs, 1960–2006 and whether it can be explained by increasing MLD. Using individual grid points from HadSST3 (1960–2006) along with MLD and wind speed data in Eq. 2 (see Methods). Analysis of SST time series at individual grid points: (A) Kendall’s τ-values for the trend in AR(1) coefficient, and (B) the estimated change in AR(1) coefficient (α) from 1960 to 2006 based on fitting a linear trend. For grid points that exhibit slowing down (increasing α): (C) the observed MLD change (1960–2005), and (D) required MLD change to explain the slowing down signal following Eq. 2. Observed (E) 1960 and (F) 2006 MLD are shown for each grid point. Observed versus estimated (required) MLD in (G) 1960 and (H) 2006. In H, outliers that are more than 2 SDs (dash-dotted lines) away from fitting the observed value are shown as red points, and the corresponding spatial locations are outlined in D.

In areas where slowing down (increasing α) is observed, the mixed layer has generally deepened since 1960, typically by 0–20 m (Fig. 3C). However, the changes in MLD required to explain the slowing down (Fig. 3D) typically exceed those observed (Fig. 3C), with a few striking regions where very large increases in MLD would be required (outlined in Fig. 3D). This result is robust to consideration of the uncertainty in MLD changes derived from the linear regression method (Fig. S6). A caveat here is that the simple model [2] chosen cannot explain the power spectrum of SST anomalies in regions strongly influenced by oceanic processes (22), including near Japan in the Kuroshio Current. However, the model [2] is valid (22) in areas of strong slowing down in the central and northeast North Pacific (Fig. 3D). A further caveat is that MLD is closely related (23) to the PDO via changes in the strength of the Aleutian Low pressure system, questioning the model assumption that MLD is an independent forcing parameter. Thus, we can partly explain observed slowing down of North Pacific SST fluctuations as due to deepening of the mixed layer (effectively giving the surface ocean a greater heat capacity), but we leave it to future work to fully explain the signal.

• Download figure
• Open in new tab
• Download powerpoint

Fig. S6.

Uncertainty in linear regression used to estimate trends in observed MLDs. SEs were used to determine the 5% and 95% confidence limits. Lower confidence limits for (A) 1960 and (B) 2006 MLDs are plotted against the estimated MLDs required to explain the slowing down signals in the grid point by grid point analysis (described in the main text). Likewise, upper limits for (C) 1960 and (D) 2006 MLDs are also plotted against the estimated MLDs. The dash-dotted lines refer to 2 SDs away from fitting the data. These plots are comparable to Fig. 3 G and H, where the mean values from the linear regression analysis of the observed MLDs are used.

To examine how slowing down of North Pacific SST fluctuations might affect the variability in marine ecosystems, we tested how two simple models behave when forced by different levels of red noise, encompassing the range of values of α found from analysis of the PDO index (Methods). The first is a bistable model with two alternative attractors, which represents the concept that some ecosystem time series react nonlinearly to climate forcing and can exhibit true regime shifts between attractors (9). The specific model chosen is generic and not based on a specific real-world system. It is set up to have no bias toward either state and to allow shifts between states to occur when α = 0.55 in the forcing time series (the lowest value considered). The second model is a double-integration model which represents the concept that some ecosystem variables have only one state and linearly track climate forcing, integrating it with their own characteristic timescale (11). The model is based on one used elsewhere to simulate Nyctiphanes simplex (a species of krill) which has a damping timescale of ∼20 mo (12) (Methods). For forcing these ecosystem models, we note that although there is an increase in variance in the climate indices we have analyzed, the increase in amplitude of variability is very modest compared with that generated by increasing α in our AR(1) model. Hence, we normalize (dividing by SD) the amplitude of the different levels of red noise (generated by different α) used to force the two idealized ecosystem models. This means that the resulting changes in the ecosystem models are due to the changing memory in climate forcing alone. The range of values of α we explore corresponds to damping timescales from ∼2.2 mo when α = 0.55, to 20 mo when α = 0.95, there being an exponential relationship between α and the timescale, and noting that α = 0.95 is often used to mimic PDO variability (10).

Increasing the autocorrelation in the climate forcing causes both the nonlinear and linear ecosystem models to become more correlated with the forcing (Fig. 4 A and B). When forcing the models with the PDO index, the correlation also increases as α increases (red lines). This is expected as the timescale of the ocean is increasing toward the assumed timescale of the ecosystems. As autocorrelation in driving SSTs increases, the SD of variations in both the nonlinear and linear ecosystem models generally increases (Fig. 4 C and D) and this is much more pronounced in the linear model (Fig. 4D). When using the PDO index to force the simple models, we again find that increases in α are generally linked to increases in SD in both ecosystem models (red lines). Examining some specific instances (Fig. S7), as autocorrelation in driving SSTs increases (Fig. S7 A–C), regime shifts in the nonlinear system become more frequent (Fig. S7 D–F), and the linear system shows increasing deviations from its single equilibrium (Fig. S7 G–I). Thus, larger ecosystem changes are associated with increasing memory in the climate forcing, especially in the linear ecosystem model, and if the ecosystem response is nonlinear, then shifts between different regimes become more frequent.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Effect of reddening climate forcing on two simple models representing marine ecosystems (Methods). Correlations between the (A) nonlinear and (B) linear model time series and the forcing time series for different values of α are shown with the 5th and 95th percentiles from 1,000 simulations at each value of α. The mean SDs of the ensemble for the (C) nonlinear and (D) linear model time series are shown with the 5th and 95th percentiles. In all four plots, the red lines show the same analysis when the original PDO index is used to force the simple models and plotting the α-value from this against the other statistics when using a moving window (Methods).

• Download figure
• Open in new tab
• Download powerpoint

Fig. S7.

Example ensemble members from the forcing of two simple model ecosystems with different levels of climatic red noise by altering the value of α (Methods). Forcing series A–C, generated from α values of 0.65, 0.85, and 0.95, respectively, are applied to the nonlinear model, creating time series D–F, and the linear model, creating time series G–I.

In summary, we detect strong slowing down (reddening) of North Pacific SST fluctuations, and of the PDO index constructed from them, over the observational record. Slowing down since 1960 can be at least partly explained by observed deepening of the ocean mixed layer. It represents a systematic change toward lower frequency, somewhat higher amplitude, North Pacific climate variability. Two of the resulting transitions in the PDO index, around 1977 and 1989, are well-known to have had significant impacts on a diverse range of ecological and climate systems (3, 4). Marine ecosystems, both those that have nonlinear dynamics (9) and those that linearly track climate forcing (11, 12), are vulnerable to large and sometimes abrupt changes in response to the low-frequency variability in the physical ocean (13). Our results suggest that ecosystem variables that respond linearly to climate variability became prone to larger changes over the observational record, as fluctuations in North Pacific SSTs slowed down. Furthermore, those ecosystem variables that respond nonlinearly to climate variability became prone to more frequent abrupt regime shifts. These results may help explain the well-known abrupt changes that occurred in North Pacific ecosystems in 1977 and 1989 (3, 4). The large size, especially of the 1977 shift, could be seen as a linear response to slowing down in climate variability, whereas the two events in relatively quick succession could be interpreted as a nonlinear response to slowing down in climate variability that was less likely to have occurred earlier in the twentieth century.

It is tempting to extrapolate forward and infer that if the trend toward increasing autocorrelation in the North Pacific ocean were to continue, the propensity for large ecosystem changes would increase. Models are generally poor at simulating observed MLD in the North Pacific (24). However, global warming is robustly expected to drive ocean stratification and a decrease in MLD over the North Pacific (25), which could potentially reverse the historical trend of slowing down.