Interannual oscillations and sudden shifts in observed and modeled climate

It has been proposed that time‐dependent synchronization of regional climate oscillations with interannual timescales could create multidecadal climate variability. Tsonis et al. (GRL) modeled the global climate as a coarse‐grained network with four nodes: El Niño Southern Oscillation, North Atlantic Oscillation, North Pacific Oscillation, and Pacific Decadal Oscillation, all but the last with significant spectral power at interannual timescales (∼2–7 years). This network seems to synchronize during turning points or “shifts” of multidecadal climate variability in the early 1910s and 1940s and late 1970s and 1990s. This article dissects those results and shows that while the synchronization idea is promising, the original implementation suffers from some issues and requires further modification. We present novel findings of striking irregularity in interannual variability. In climate model simulations, the rising and falling components between local minima and maxima of global annual mean temperature are normally highly correlated, but exhibit significant anticorrelation every 50–80 years. The new results are a step forward in understanding key features of internal climate variability.


| INTRODUCTION
Internal fluctuations and oscillations of the atmosphereocean system are a large part of overall climate variability. The impact of El Niño Southern Oscillation (ENSO for short from now on) on interannual global mean temperature variability is well known. Further, internal inter-or multidecadal variability was likely part of the slower rise of global mean temperature from 1998, until recently new record temperatures were again much higher than previous ones (Trenberth, 2015).
The timescale ENSO operates at is 2-7 years, which we will from now refer to as interannual timescales. Another major peak in spectra of observations as well as modeled climate is seen at multidecadal timescales of about 50-80 years (Schlesinger and Ramankutty, 1994;Knight et al., 2005;Henriksson et al., 2012). Yet, there are many unanswered questions regarding even the basic physical mechanisms behind the multidecadal variability. In the North Atlantic, the multidecadal variability is called Atlantic Multidecadal Oscillation (AMO) and it is known to be correlated with the Atlantic Meridional Overturning Circulation (AMOC) (Medhaug and Furevik, 2011). In the North Pacific, variability at multidecadal timescales is manifested in a typical pattern called Pacific Decadal Oscillation (PDO). However, the mechanisms maintaining these oscillations or fluctuations are not well known. The AMO has been linked to delayed density perturbations by advected tracers, especially salinity arriving at key downwelling areas of the Atlantic circulation (Delworth et al., 1997;Kerr, 2000;Latif et al., 2004;Vellinga and Wu, 2004;Jungclaus et al., 2005;Dima and Lohmann, 2007;Keenlyside et al., 2008;Frankcombe and Dijkstra, 2011). However, this picture seems to need modification because of irregularity of this behavior from cycle to cycle and fairly weak evidence shown in delayed correlations (Henriksson et al., 2012;Zanchettin et al., 2014).
Complex networks, an emerging subfield within climate science (e.g., Donges et al., 2009;Donner et al., 2009;Steinhaeuser et al., 2011) has been applied in attempts to explain multidecadal variability better. As a pioneering idea, Tsonis et al. (2007), Swanson and Tsonis (2009) and Wang et al. (2009) proposed to model the climate as a coarsegrained network consisting of ENSO, North Atlantic Oscillation (NAO), North Pacific Oscillation (NPO) and PDO. In their research, the global climate network experiences periodic synchronization events that lead to a "shift" in multidecadal climate variability. Couplings between ENSO variability in the tropical Pacific and variability in the North Atlantic have also been found in more physical terms (Brönnimann, 2007;López-Parages and Rodríguez-Fonseca, 2012). Additionally, couplings between interannual and multidecadal and longer timescales have been suggested in other frameworks (Palus, 2014a;2014b;García-García and Ummenhofer, 2015;Griffiths et al., 2016). Should the mode synchronization mechanism be true it would bear major applicability for climate variability attribution and perhaps potential for climate forecasting. While the hypothesis is not directly contradicting the delayed density perturbation theories, there is hardly any overlap either. Therefore hypothesis testing should be possible, even if the observational record is relatively short for this purpose. While the salinity advection and network synchronization mechanisms are seemingly very different, they could also work simultaneously and even be coupled together. The aim of this article is to evaluate the truthfulness of network-synchronization-related hypotheses in more detail.
In the first part of this article, we dissect the results of Tsonis et al. (2007) based on Niño3.4, NAO and NPO data, which are all available from 1899 onward. 1 We also add Southern Annular Mode (SAM) data from 1957 and discard PDO data for reasons that will be made clear in the next section. Additionally we will employ the COSMOS climate model (Jungclaus et al., 2010) that we evaluated, in our opinion successfully, to reproduce most of the main qualitative patterns of observed multidecadal variability in Henriksson et al. (2012). We will see that the well-founded dropping of the PDO from the considered network weakens the result somewhat and that the network synchronization mechanism hypothesis requires modifications.
In the second part, we'll therefore view the connections between interannual and multidecadal timescales from a new angle. We'll evaluate regularity and irregularity of interannual variability in different phases of multidecadal climate variability. The results shed light on hitherto unnoticed features connecting variability at the different timescales and will be important in figuring out what exactly is going on in multidecadal climate variability.
where the time-dependent distance between oscillators i and j was defined as d t ij ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 1−jr t ij j r . d(t) as a sum over all pairs of oscillators is then the aggregate network distance and, as a generalization of correlation to more than two variables, a measure of synchronization. The time-dependent (linear, Pearson) correlation r ij (t) is evaluated using a moving time window with length Δt = 11 years (with the time t being the midpoint of each interval), with small changes to window length reportedly not causing a large change. Tsonis et al. further defined more refined quantities (symbolic phase and phase prediction). Further quantities can be derived and might provide extra information for example in exotic cases, where phases correlate better than amplitudes, noise in climate data can be notoriously complex and more complicated models can become black boxes introducing new biases (see e.g., Mudelsee, 2010). More elaborate analyses could, however, attempt to extract causality not only between amplitudes of different oscillations but also between phase of one and amplitude of another (Palus, 2014a;2014b;Hlinka et al., 2017). However, by examining Figure 1 of Tsonis et al. (2007), it is clear that the main conclusions can already be drawn from the network distance and while dissecting the results we stick to the network distance.
For considerations involving Equation (1) to make sense the variables will have to be truly oscillatory. As implied already by the name, the PDO has too much energy in the low frequencies (Newman et al., 2016) and will show up as a trend in an 11-year time window, violating standard correlation analysis assumptions. Additionally, PDO is correlated with both the NPO and with ENSO, manifested in the tropical and North Pacific. The corresponding correlation coefficients are −61 and−55%, respectively, for the whole available time series (of length 116-117 years). It is well known that there are connections between PDO, NPO and ENSO, while some of the physical mechanisms and consequences are being debated (e.g., Gershunov and Barnett, 1998;Wang et al., 2008). Including such permanent correlations will merely shift the network distance to a lower value without any actual synchronization taking place. All other correlations considered are very low, with the exception of the correlation between ENSO and NPO, −37%, still statistically significant. We decide to keep the NPO as it is oscillatory at the window's timescale and still contains a lot of independent variability.
In addition to removing the PDO, we will add the SAM to the considerations to allow for possible synchronization also with the southern hemisphere beyond mere ENSO variability. With these modifications of the analysis, we will be able to reproduce and dissect the previous results and to take the research one step further. Further research could also attempt to apply state-of-the-art network theory to try to reevaluate the whole network and its ideal coarse-grained constitutes.
As shown in the results of Tsonis et al. (2007), the global network synchronizes in the early 1910s, the early 1940s, the late 1970s and the late 1990s. Global mean temperature features qualitative "shifts" with the multidecadal variability changing phase from cooling to warming (superimposed on long-term warming) and vice versa during the same times. Two other, and weaker, synchronization events were found in the 20th century, but based on additional metrics they were not considered strong synchronization events. In the following subsection, we show what the apparent synchronization events 2 are composed of, in terms of timedependent pairwise correlations of the different climate modes of variability.

| Dissecting the network in observations
The pairwise correlations of the different climate indices as well as network distance and an extra plot with monthly correlations are shown in Figure 1. The synchronization events of the 1910s and 1940s can also be seen clearly in the analysis discarding the PDO. The 1970s and 1990s synchronization events are clearly weaker when the PDO has been dropped.
We see that the 1910s synchronization event involves a positive correlation between ENSO and NAO and negative correlations between NPO with both ENSO and NAO. The 1940s synchronization events involve negative correlations of ENSO with both NAO and NPO and a positive correlation between NAO and NPO. The least clear event of the 1970s involves all three considered correlations being negative first with the NAO and NPO correlation suddenly becoming positive. The late 1990s event is not as strong as the 1910s and 1940s events either, suggesting that there might be more to the story. The intermonthly correlations valuated year by year between ENSO and NAO are similar to the interannual results and suggests that the results are consistent over another timescale. Adding the SAM to the considerations does not change the picture and the network distance is almost the same. In the top subplot in Figure 1, the correlation of SAM with NAO and NPO are omitted for clarity, but they contribute to the network distance.
In summary, we reproduce most main features of Tsonis et al. (2007), Swanson and Tsonis (2009) and Wang et al. (2009), but in a weaker fashion. Dropping the decadal PDO mode weakens the results and suggests the mechanism for multidecadal variability is not quite as simple as four oscillators synchronizing. In general, dropping the PDO reveals that the network is not as synchronized as thought before, even though the limit for reaching statistical significance also increases when accounting for PDO's correlations. The 95% statistical significance value in Figure 1 is produced through a four-dimensional AR-1 process with cross-index correlated Gaussian noise following Tsonis et al. (2007). Reproducing the statistical significance analysis assuming correlated AR(1) noise with the correlation matrix calculated for the full time series-without the PDO-provides a slightly higher 95% significance limit of about 1.01, which would again indicate higher synchronization compared to random incidents, except for the 1970s and 1990s events. On the other hand, by varying assumptions as allowed by the short time series data in statistical bootstrapping with longer memory, including allowing for oscillatory memory, the statistical significance limit can be lower too, down to about 0.95 (results not shown).
The interannual irregularity considerations add some information about the "shift" events. The results shown here also bring to light the important follow-up question of what determines the sign of the different correlations during different events. There is no consistency regarding different turning points of multidecadal variability and the signs change from event to event. Additionally, there is no physical explanation for the synchronization mechanism yet, even though physical coupling mechanisms between ENSO and NAO have been proposed (see e.g., Wu and Lin, 2012 and references therein). Further possibilities also remain: some causality might run in the opposite direction from multidecadal to interannual timescales or there could be complex phase-amplitude couplings such as proposed at faster timescales in an information-theoretic framework by Palus (2014a;2014b) and Hlinka et al. (2017). The earlier reference (Palus and Vejmelka, 2007) concludes on a general level that limited datasets are more suitable for detecting synchronization than causality between timescales. They additionally present a seven-item list for challenges and pitfalls in such models tracking causality and its direction, for which reason we leave such considerations outside the scope of this article and propose it as future research for information theory specialists.

| The network in the COSMOS control simulation
The COSMOS earth system model, consisting of the ECHAM5 atmospheric model, the MPI-OM ocean model and other components (Henriksson et al., 2012) was evaluated for its internal multidecadal variability against observations. Considering the limitation that internal variability never happens in isolation from externally forced variability in the physical climate, the model performed very well. We can examine the physics behind multidecadal variability in the model to form hypotheses to be tested in observed climate. The following results are produced with data from the same 1201-year unforced simulation as the results in Henriksson et al. (2012). At the end, we also perform an analysis on the forced 1205-year simulation presented in that article. The forced simulation includes all known important forcings of the 20th century: greenhouse gases, anthropogenic aerosols, volcanic aerosols, changes in incoming solar radiation, and land use. The available time series are an order of magnitude longer than the observed ones. The climate variability itself is described in more detail in the reference. Figure 2 shows a 300-year time slice of the correlations between the different climate indices examined above for the observations. We see synchronization at least around 1690, 1720, 1760, 1800, 1830, 1890, 1910, 1940 and 1970. All events are accompanied by a turn in multidecadal variability, except the 1940 one. 3 However, the coupling strength and its temporal development vary from event to event, which is why we take a completely new angle on examining shift events in the next section.

| SHIFT EVENTS AND INTERANNUAL IRREGULARITY
The global annual mean temperature time series contains local minima and maxima. A maximum occurs at time t when the temperature is higher than at its adjacent points in time: T(t − 1) ≤ T(t) ≥ T(t + 1). It is possible to form a new time series consisting of the rising parts of this time series and the falling. If t n max is the nth maximum and t n min is the nth minimum, then f n ð Þ ¼ T t n max À Á −T t n min À Á and r n ð Þ ¼ T t n min À Á −T t n + 1 max À Á , when falls are denoted by f(n) and rises by r(n). 4 Figure 3 shows the correlation between the time series consisting of the falls with that of the rises. This measures the regularity of the oscillation, whether a largerthan-average fall in temperature gets followed by a corresponding larger-than-average rise or the opposite. A moving time window of seven local maxima and minima was chosen, making the length of the time window about 23 years on average.
In our 1201-year simulation, we have 366 local maxima and minima. We can see that the correlation between rises and falls is close to unity during the major part of the simulation with sporadic dips to even large negative values occurring. By setting a threshold for correlation, here chosen to be +0.2, we can identify these rare episodes, shown in the central subplot of Figure 3. Where the correlation stayed below Unforced COSMOS simulation data, 300-year time slice. Top: correlation between falls and rises between local maxima and minima of global annual mean temperature. Middle: episodes, where correlation dips below +0.2. Bottom: same correlation as in top plot, global annual mean temperature with multidecadal-pass filtered smoothened fit (K). Strongest anticorrelations between falls and rises marked by vertical lines the threshold for some time, these time steps were clustered together and the midpoint of each episode plotted. We observe a total of 16 such events, meaning the average separation is about 75 years. This is close to the dominant timescales of 50-80 years for multidecadal climate variability both in the model simulation and in observations, as concluded from spectral analysis (Henriksson et al., 2012). When we compare the timing of the negative correlations to the multidecadal medium-pass filtered timeseries in the lowest subplot, we see that most such events are followed by a turn in the multidecadal variability. On the other hand, not every turn is associated with an anticorrelation, a feature which should be studied further. Nevertheless, the connection between anticorrelations in interannual variability and turns in multidecadal variability is a novel finding that could lead us forward toward breakthroughs in understanding internal climate variability.
Finally, let us take a look at the same data in observations and in the long forced simulation presented in Henriksson et al. (2012) as plotted in Figure 4. In the observations, the correlation between the falls and rises between local maxima and minima has clear drops during the "shifts" of the 1880s, 1940s and 1990s. As in the unforced simulation, not every shift is associated with a turn in the multidecadal variability. The forced simulation shows similar traits as the unforced one with anticorrelations occurring before shifts in multidecadal variability. Differences with measurements have to do with the fact that the phasing of internal variability is not aligned in the two datasets. 5 Additionally, ENSO is too strong in the used model version (Jungclaus et al., 2006). Nevertheless, qualitatively the picture is the same in observations and both simulations.
As a final note, let us recall from Henriksson et al., 2012 that the estimated amplitude of the multidecadal variability is about 0.05-0.15 K, which is a large part of observed variability between decades, but small compared to long-term greenhouse gas warming.

| CONCLUSIONS
We dissected the climate oscillation network consisting of ENSO, NAO, NPO and SAM. The oscillations seem to have some degree of synchronization during turning points of multidecadal climate variability, even though it is weak for the late 1970s and 1990s, while clearer for the 1910s and the 1940s. It seems that the (Tsonis et al., 2007) results indicated that network synchronization events were stronger than they really are, because of the inclusion of PDO that distorts the analysis through low-frequency behavior and strong correlation with other oscillations. Additionally, it is unclear what determines the sign of the correlations as they change from event to event. The synchronization hypothesis still seems to need refining. Top: correlation between falls and rises between local maxima and minima of global annual mean temperature as well as global mean temperature (K), observations. Bottom: forced COSMOS simulation data. Correlation between falls and rises and global annual mean temperature with multidecadal-pass filtered smoothened fit (K). Vertical line showing the strongest anticorrelation between fall and rise A novel finding was the irregularity in rises and falls between local minima and maxima of global annual mean temperature every 50-80 years. Most of the times this irregularity also occurred just before "shift" events, where multidecadal cooling and warming alter. These irregularities did not occur strongly during every "shift"-as is the case for network synchronization-but both still occur consistently enough throughout the datasets to warrant further attention.
As outlined in this article, there remains plenty of scope for further research. The physical manifestation of the mechanisms is one obvious topic to study. Further, more elaborate methods such as those of Palus (2014a;2014b) can be applied while keeping in mind the many complexities and pitfalls involved. The frameworks applied so far all seem to extract some essential information while also leaving some ambiguities.
Understanding the connections between interannual and multidecadal climate variability have the potential to lead to breakthroughs regarding internal climate variability. The physics of the reported findings present plenty of opportunities for interesting further research.