Multi-Decadal Trends of Global Surface Temperature: A Broken Line with Alternating 30 yr Linear Segments?

We investigate global temperature data produced by the Climate Research Unit at the University of East Anglia (CRU) and the Berkeley Earth Surface Temperature consortium (BEST). We first fit the 1850-2010 data with polynomials of degrees 1 to 9. A significant ~60-yr oscillation is accounted for as soon as degree 4 is reached. This oscillation is even better modeled as a broken line, more precisely a series of ~30-yr long linear segments, with slope breaks (singularities) in ~1904, ~1940, and ~1974 (±3 yr), and a possible recent occurrence at the turn of the 20th century. Oceanic indices PDO (Pacific Decadal Oscillation) and AMO (Atlantic Multidecadal Oscillation) have undergone major changes (re-spectively of sign and slope) roughly at the same times as the temperature slope breaks. This can be interpreted with a system of oceanic non-linear coupled oscillators with abrupt mode shifts. Thus, the Earth’s climate may have entered a new mode (a new ~30-yr episode) near the turn of the 20th century: no further temperature increase, a dominantly negative PDO index and a decreasing AMO index might be expected for the next decade or two.


Introduction
Global surface temperatures are one of the parameters most commonly used to discuss the evolution of climate. Databases of instrumental temperatures that cover the past century and a half have been compiled by four main groups (Climate Research Unit at the University of East Anglia-CRU, NASA Goddard Institute for Space Studies-GISS, National Oceanic and Atmospheric Administration-NOAA, Berkeley Earth Surface Temperature consortium-BEST). This is a huge, difficult task and these databases are necessarily faced with a number of limitations: the geographical distribution of stations is far from uniform and varies with time; also, there may be fundamental difficulties in establishing meaningful global temperatures for the Earth [1]. Climate is generally defined as the ~30 year average of weather and since only 100 to 150 years of instrumental data are available, this places severe limitations on the significance of global temperature trends and multi-decadal oscillations. A monotonic (low degree) trend can be fit to all global tem-perature data sets over the period from 1850 to the present. This trend amounts to a secular rise of ~1 K over the period. Once this trend is removed, a ~60 yr oscillation stands out in most records. Such an oscillation has been discussed for instance by [2] and [3]. Lean and Rind [4,5] have built a model curve for 1889-2006 monthly mean global temperatures, in which they distinguish contributions from oceans (using the El Nino Southern Oscillation ENSO), volcanic aerosols, solar irradiance and anthropogenic forcing. Both observations and their empirical curve display the secular rise of approximately 1 K. But data are lower than model by some 0.2 K around 1910 and higher by 0.2 K around 1940: this is because of the significant ~60 year oscillation identified by previous authors, which is not accounted for by [4].
In this paper, we first check whether some of the main global temperature data sets can be reasonably fit by smooth polynomials of increasing degree from 1 (linear trend) up to 9, with particular focus on the "~60-yr oscillation". We next suggest that this oscillation may be fit by a series of ~30-yr long linear segments with rather ab-rupt changes in trend. We find correlations between the breaks in the multi-decadal trends in global surface temperature and changes in sign of the Pacific Decadal Oscillation and slope of the Atlantic Multidecadal Oscillation indices. We discuss our results in the frame of the dynamical mechanism for major climate shifts proposed by [6].

Polynomial Fits to the "~60-yr Oscillation" in Global Temperatures
We have used monthly temperature anomaly data from the CRU and BEST databases. We chose to use the "raw data" rather than other data-derived products that have been "homogenized" and have thus undergone rather extensive and not always transparent/reproducible data polishing.  [8]). All yield similar results and mainly results based on CRU-Tem4vGL are shown in this paper. We note in passing differences between the long-term linear or low-degree trends obtained in the various temperature anomaly databases (notwithstanding a possible overall difference in baseline due to the definition of the period with respect to which the anomalies are calculated); the CRU and GISS curves are about 0.3 K below the NOAA and BEST curves in the 2000-2010 decade, whereas all agree to better than 0.1 K between 1920 and 1970 [8]. We return to this in Section 5.
In order to check the significance of the "60-yr oscillation", we have very simply fit the data sets of monthly mean temperature anomalies (using least squares) with polynomials of increasing degree (Figure 1). The linear (n = 1) trend fit to the CRUTem4vGL data shows the secular rise of approximately 1 K. For n = 2 and 3 the curves fail to account for the "60-yr oscillation" (Figure  1, purple and blue curves). As soon as the degree of the polynomial reaches 4, the fit accommodates the "60-yr oscillation", after which the situation remains stable (shown up to n = 9 in Figure 1). Extrapolations of these trends outside of the data domain show quick divergence and are of course meaningless. A ~60 year oscillation would occur only 2.5 times over a 150 yr interval and could not be accounted for by a polynomial of degree less than ~4.
In order to quantify the goodness of fit, we have calculated the adjusted R-squared values not only for polynomial fits to the data with degrees 1 up to 9 but also for an exponential fit. The adjusted R-squared [9] is a statistically unbiased version (i.e. corrected for the number of degrees of freedom) of the R 2 "coefficient of determination", whose main purpose is the prediction of future outcomes on the basis of prior information. R 2 itself is 1 minus the ratio of variance of errors to variance of observations. The adjusted R-squared values take into consideration the number of free parameters of the model (the adjusted R-squared values may be negative when the parameters of the model do not improve the fit compared to the simple average value of the data, which happened in the case of one data set, CRUTem4vSH, that could not be adequately fit by an exponential). In the case of CRU-Tem4vGL, the adjusted R-squared increases significantly from 0.36 (n = 1) to 0.50 (n = 4) and remains stable thereafter (0.52 for n = 9). As soon as the degree exceeds "~60-yr fluctuation" is the topic of this section. Figure  2(a) shows another display of the CRUTem4vGL monthly data set, limited to the period from 1900 to 2010 (the previous 50 years have significantly larger variance, which is likely due to insufficient, changing station distribution and also possibly to lesser data quality). We note that a series of continuous linear segments with abrupt changes in slope may provide a good fit: such fits have appeared in the literature (e.g. Figure 8 of [10], where linear regression lines are shown for intervals 1901-1934, 1934-1979 and 1979-2010). The most recent has been proposed and discussed by [11], using methods from the dynamics of synchronized chaotic systems.
The models we are seeking belong to the class of "change-point models" [12]. Vasko and Toivonen [13] for instance have published quantitative algorithms to estimate the number of linear segments to be fit to a time series using permutation tests. We have constructed an algorithm that produces such models, paying particular attention to the management of secondary minima using non-linear simulated annealing. Simulated annealing uses the Metropolis algorithm (see [14][15][16], for description of the algorithms and definition of terms; see similar use of simulated annealing in [17]). The method is outlined briefly in the Appendix.
We have first tested our method on synthetic data consisting of N = 2 to 10 linear segments with variable durations and slopes and varying noise levels from 5% to 25% of the total amplitude range of the data. For each value of N, and for each step of the control parameter, we ran 100 attempts at locating the nodes and displayed the histograms of node locations. For 4 segments and 5% noise, the 3 nodes are identified to within a year in a 2000-month (166 yr) long test series; for 15% noise, detection of the 3 nodes is unambiguous but the histograms widen to ~100 months (±4 years uncertainties). For 25% noise, the nodes are still detected but less prominently and with larger uncertainties (~±8 years).
We have then applied the method to the 1900-2010 CRUTem4vGL monthly data, using from 3 to 10 segments. Results are shown in Figure 2(a) for the best case, i.e. with 5 segments, and in Figure 2(b) with 4 segments for comparison. Two nodes are clearly identified at ~1940 and ~1974. Two other nodes at ~1904 and ~2006 are near the edges of the data distribution and could be artefacts (end-effects), though this has not been observed with the tests on synthetic data. The ~1904 node appears to be supported by visual inspection of some data predating 1900 (but data dispersion is significantly larger before 1900 than afterwards; see Section 5). The more recent one may require at least a decade of additional data to ascertain its validity. In the case of 4 segments (Figure 2(b)), the histogram does confirm the same 4 nodes with the same dates, although of course each indi- underlines the same nodes, but with more spread in dates. As the number of segments is increased to 10, node distributions become more uniform, as segments are used to adjust to shorter (sub-decadal) fluctuations in the data. The differences in location of the nodes introduce some curvature in the "cloud" of 500 fitted broken lines, but the linear segments between the nodes are quite clear (Figures 2(a) and (b)). We confirm our preliminary hypothesis that the 1900-2010 data can be fit well by such a series of linear segments and that the data argue in favor of at least two singular events around 1940 and 1974. The slopes of the three main segments are 1. . This oscillation appears to be better represented by a sequence of linear segments with rather abrupt slope changes at the ~1904, ~1940 and ~1974 nodes (adjusted Rsquared > 0.62), rather than a sinusoidal or polynomial fit (see previous section). Climate shifts close to these dates have independently been identified by [11], using methods from the dynamics of synchronized chaotic systems (see Section 5).
In order to get a further idea of the robustness of inferences that can be made from Figure 2, we have resumed the same computations as in Figure 2(a) (5 segments, 4 nodes) for all 21 datasets that can be obtained from the CRU site http://www.cru.uea.ac.uk/cru/data/temperature/#datter (i.e. land air temperature anomalies CRUTem3, CRU-Tem3v, CRUTem4 and CRUTem4v, sea surface temperature anomalies HadSST2 and combined land and marine temperatures HadCRUT3 and HadCRU3v, each with hemispheric means for the northern and southern hemispheres and combined global series), and we have stacked all corresponding histograms, resulting in

Comparing Trends in Global Temperature and the PDO and AMO Oceanic Indices
We have applied the same method to an oceanic proxy, the Pacific Decadal Oscillation (PDO is the leading principal component of monthly sea surface temperature anomalies in the North Pacific Ocean, poleward of 20˚N). The optimal number of segments is 6 (Figure 2(c)). Two prominent doublets of nodes separated by ~33 years are seen, one at ~1940 and ~1948, the other at ~1973 and ~1982. A less prominent node is located at ~1919 (29 years before the 1948 cluster). The fitting exercises above suggest a comparison, which is illustrated in Figures 4(b) and (c), where the monthly values of PDO and global temperature anomaly are shown on the same time scale. Periods when the multi-decadal slopes of linear segments fitted to the temperature data, as outlined in Figures 2 and 3, were positive and negative have been colored, respectively, in pale red and blue. There appears to be a strong correlation between the sign of the slope of T (i.e. its multidecadal time derivative) and the dominant sign of PDO. There is good correspondence between the ~1940 and ~1974 nodes of CRUTem4vGL and the main changes in PDO sign, with PDO shifting from mainly positive (between ~1930 and ~1940) to mainly negative (between ~1950 and ~1976), then back to mainly positive (between ~1976 and ~2000).
The Atlantic Multidecadal Oscillation index (AMO characterizes the average sea surface temperature over the northern Atlantic; it is defined as the SST averaged over 0˚N -60˚N, 0˚W -80˚W, minus SST averaged over 60˚S -60˚N) also correlates well with global mean surface temperature anomalies (as for instance already pointed out by [8]), and therefore also with the dominant multidecadal sign changes of PDO, as seen in Figure 4(d): a decreasing linear segment prior to ~1920, an increasing one from ~1920 to ~1940, a decreasing one from ~1950 to ~1975, an increasing one after ~1975, generating a strong "~60-yr oscillation".

Discussion
oceanic indices PDO (and AMO) and in global surface temperature since the early 20th century, more precisely a causal connection between changes in the modes of oceanic circulation and changes in multi-decadal (linear) trends of global surface temperature. Given the respective masses, impedances and time constants involved, it is reasonable to argue that the oceanic system forces the atmospheric system on these time scales rather than the opposite.
Tsonis et al. [6] have proposed a dynamical system approach of major climate shifts, such as the ones identified by the nodes at ~1940 and ~1974. These authors define a network of oceanic indices (PDO, the El Nino Southern Oscillation ENSO, the North Atlantic Oscillation NAO and the North Pacific Oscillation NPO). They find that over the period 1900-2000, this network synchronized several times. [6] define the coupling strength between indices and find that "where the synchronous state was followed by a steady increase in the coupling strength (···), the previous synchronous state was destroyed, after which a new climate state emerged". The dates of the events they find are near 1910, 1940 and 1975, i.e. essentially identical (given data variance and uncertainties in node locations) to the dates we have identified for CRUTem4vGL temperature anomalies and also most other temperature data sets (Figure 3).
[6] find the same features in models of "systems with nonlinear coupled oscillators, caused by bifurcations as the coupling parameter changes". They cannot conclude (nor can we) whether changes are triggered by some external forcing or are generated within the chaotic system itself. Their figure 2(d) illustrates the shifts in the global-SST ENSO index: this is in essence the same "boxcar" correlation that can be observed in our Figure 4 between T and PDO; in the present paper, we believe we provide stronger and more quantitative observational evidence of the reality of these shifts and when they occur.
A recent (~2006) shift of the temperature anomaly slope is suggested by Figures 2 and 3. In the early 2000s, the PDO turned back to a negative-dominated mode (Figure 4(c)) and the slope of the multidecadal linear trend in AMO may have become flat, as had happened towards the end of the ~1920-1940 episode (Figure 4(d)). Also, the succession of ~30 yr long quasi-linear segments in the temperature anomaly curve (Figure 4(b): possibly 1860-1880 (?) and 1880-1910 (?), then 1910-1940, 1940-1970 and 1970-2000) could itself suggest that the ~60 yr oscillation, whatever its source, may continue as a ~2000-2030 segment.
We point out the puzzling fact that an earlier CRU global temperature dataset, HadCRUT3GL, displayed the alternating segments better than the latest version CRU-Tem4vGL (Figures 4(a) vs (b)). The annual variability of the monthly averages increased by a factor of 2 in the recent revision. So did the amplitude of the temperature anomaly since 1975, that changed from 0.5 K to about 1.0 K. The segment with negative slope prior to ~1910-1920, and more importantly the segment with flat or even negative slope after 2000, are particularly conspicuous on HadCRUT3GL (Figure 4(a)). [8] discuss the fact that the four available global temperature anomaly curves start diverging in the past two decades (their Figure 8): there is close agreement from 1900 to the late 1980s, but after that GISS and HadCRU display a plateau, whereas the NOAA and BEST curves continue to rise to values that are 0.25 K above the GISS and HadCRU plateau. It is awkward that it is the most recent part of the data compilation that displays this divergence. The plateau after the late 1990s has been the topic of much recent discussion; the previous HadCRUT3 data display this plateau much more clearly than CRUTem4vGL that seems to have suppressed it.
The shift in dominant sign of PDO (Figure 4(c)), the flattening in AMO (Figure 4(d)), and the plateau in several global temperature anomaly data (e.g. Figure 4(a)) taken together would support the hypothesis of a regime change having occurred near the end of the 20th century or early in the 21st century. A ~2000 regime change would complement the series of events at ~30-yr intervals (~1910-1920, ~1940, ~1970, ~2000). This change, if real, could imply that a ~15 year long period with no further warming lies ahead. Klyashtorin and Lyubushin [2], de Jager and Duhau [18], Scafetta [3] and Russell et al. [19] are among several authors who have suggested decades of future cooling: these authors argue for external forcings due to solar or planetary effects (see also the recent paper by [20]), when the mechanism of [6] envisions the chaotic result of interaction between oceanic non-linear dynamical systems. Tsonis et al. [6] conclude that the 1970 to 2000 warming may not be (wholly) due to the radiative effect of greenhouse gases overcoming shortwave reflection effects due to aerosols and that the climate may indeed have shifted to a different state.
Our analysis strongly argues for the presence of a ~60yr oscillation in the climate system (at least over the limited time interval-100 to 150 years-covered by reliable instrumental observations). More precisely, global temperature data can be interpreted as a series of linear segments interrupted by rather fast changes in slope, i.e. abrupt changes in regimes. Each episode or regime maintains itself for approximately 30 years. We favor the oceanic system as a driver of atmospheric temperatures on these multi-decadal time scales. It remains to be seen whether abrupt changes in climate mode are a result of internal chaotic dynamics of the ocean system only or could be forced by external factors. Copyright

Acknowledgements
We acknowledge financial support from IPGP as part of the IEPT RAS-IPGP cooperation. IPGP Contribution NS 3391.