Characteristic Time Scales of Decadal to Centennial Changes in Global Surface Temperatures Over the Past 150 Years

Abstract

We apply singular spectral analysis (SSA) to series of monthly mean values of surface air temperatures T, International sunspot number (ISSN), and polar faculae PF (1850–2017 for T and ISSN). The efficiency of the SSA algorithm that we use has been regularly improved. For the T, ISSN, and PF series, the SSA eigenvalues and first components are shown with their Fourier spectrum. Components of T, ISSN, or PF share similar periods. Most are found in solar activity. The ~22‐ and ~11‐year components are modulated and drift in phase, reflecting slight differences in spectra. On the shorter‐period side, components at ~9, ~5.5, and ~4.7 years are in good agreement. They have been identified in solar activity. The 60‐year component is prominent in T. It is not immediately apparent in ISSN but can be extracted with an appropriate choice of SSA window. Other types of data allow one to explore longer periods and confirm climatic variations at ~60, ~35, and ~22 years and at 50–150 and 200–500 years. When we consider a longer ISSN series starting in 1700 and recalculate the SSA first component, the trends of solar activity and temperature over the time span from 1850 to 2017 are very similar, with slower rise before 1900 and after the late 1900s, separating a faster rise in much of the twentieth century. These trends, extracted over only 150 years, could be parts of longer, multicentennial changes in solar activity. Much of the variability of surface temperatures could be linked to the Sun.

Introduction

In a series of previous papers, we have analyzed centennial series of a number of geophysical and physical parameters using an adapted version of the powerful method of singular spectrum analysis (SSA). We found strong evidence of solar signatures in the Earth's mantle rotation pole position (Lopes et al., ), in several climatic indices such as Pacific Decadal Oscillation, Madden-Julian Oscillation, and Atlantic Multidecadal Oscillation (Le Mouël et al., ) and in the length of day (Le Mouël et al., ). Encouraged by these results, we apply in this paper the same method to the series of global surface temperatures provided by a number of research teams. More precisely, we analyze the monthly mean data from the following: (1) the Hadley Research Center and the Climate Research Unit, HadCrut (https://crudata.uea.ac.uk/cru/data/temperature/) and CrutemV4 (https://www.metoffice.gov.uk/hadobs/crutem4/data/download.html), (2) National Aeronautics and Space Administration, GISS (https://data.giss.nasa.gov/gistemp/), and (3) Berkeley, BEST (http://berkeleyearth.org/data/). Many papers have commented on the reliability of these data, which we do not attempt here. We take the published data at face value and submit them to SSA in the same way we have done in the above cited papers (see also the appendix). The original method is now rather widely known; in order to avoid repetition, we refer the readers to summaries of the method in our previous papers, to the monograph by Golyandina et al. (), and to accounts by Vautard and Ghil (), Ghil et al. (), and Vautard et al. (). Many geophysical time series are “short” in a numerical sense; that is, their length is not very much longer than some of the periodicities that they might contain. Moreover, they actually often contain pseudoperiodicities that fluctuate in both amplitude and pseudoperiod rather than truly periodical variations. SSA provides at the same time a noise reduction technique, a detrending algorithm and a way to identify oscillatory components. SSA has been applied to the irregular El Niño–Southern Oscillation phenomenon, to geopotential height data, to a number of indicators of climate variability, to solar observations, and to cosmogenic isotopes. We present here results obtained by our SSA treatment of the global temperature anomalies HadCrut, CrutemV4, GISS, and BEST. We illustrate the results mainly with the HadCrut series that have been in general confirmed by the analyses of the three other series. We have progressively improved the efficiency of SSA in the algorithms we use (see Appendix 2 of Le Mouël, Lopes, Courtillot, & Gibert, , and the appendix at the end of the present paper): In the way we implement SSA, we use large computer power. The Hankel matrix we calculate has 20 million terms. We work with 64-bit computers and double precision numerics of 8 bits each. Singular value decomposition of that matrix yields two eigenvector matrices and one matrix with eigenvalues that have the same dimensions. We therefore require an 8-Go active memory, which means we must run on a computer with at least 10 Go. This is what we do and one reason why we obtain the new results described below (for more details see the appendix).

The HadCrut Series From 1850 to 2017

Figure  (top) shows the full HadCrut data set from January 1850 to January 2017 at a monthly sampling rate (in red). Figure  (bottom) shows its Fourier spectrum, with in particular significant energy around 60 and 90 years. Figure  shows the SSA eigenvalues ranked in order of decreasing amplitude and the following figures show the seven first SSA components in the order of decreasing pseudoperiod, together with their Fourier spectrum. The amplitudes of eigenvalues versus the main frequencies of eigenvectors can be considered as estimates of the spectral density, for monotonic spectral densities. Figures  show these eigenvectors and allow one to calculate their power, the frequency being evaluated from the peak value of the component's Fourier transform (or mean in the case of two peaks, which happens only once, in the case of BEST; Figure ).

Figure 1. (top) The HadCrut series of global surface mean temperature anomalies from January 1850 to January 2017 (https://crudata.uea.ac.uk/cru/data/temperature/; monthly mean values). (bottom) Its Fourier spectrum.

Figure 2. The first 32 singular spectrum analysis eigenvalues and seven components of series HadCrut (see Figure )

Figure 3. (top). The first singular spectrum analysis component of series HadCrut (see Figure ), CrutemV4, GISS, and BEST, that is, their trends. (bottom) Their Fourier transforms.

Figure 4. (top). The second singular spectrum analysis component of series HadCrut (see Figure ), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.

Figure 5. (top) The fourth singular spectrum analysis component of series HadCrut (see Figure ), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.

Figure 6. (top) The sixth singular spectrum analysis component of series HadCrut (see Figure ), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.

Figure 7. (top) The third singular spectrum analysis component of series HadCrut (see Figure ), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.

Figure 8. (top) The fifth singular spectrum analysis component of series HadCrut (see Figure ). (bottom) Its Fourier transform.

Figure 9. (top) The seventh singular spectrum analysis component of series HadCrut (see Figure ), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.

Figure  shows not only the trend for the HadCrut data but also the trends found for the three other time series GISS, BEST, and CrutemV4. The dispersion of the four curves gives a measure of the reliability (or lack of reliability) of the estimates, given that they should all refer to the same physical parameter, that is, global surface temperature. All curves show roughly a flat trend before 1900, a positive trend from 1900 to 2000 and a flatter trend since then (an indication of a “plateau”). In the main rising part of the curves, GISS has a steeper slope than the other three curves and a rather different spectrum above 40 years in period. Figure  shows the next component of HadCrut, with a pseudoperiod of 68.7 ± 23 years. BEST provides a similar value at 67.2 ± 27.9 years, whereas GISS and CrutemV4 are offset at 62.6 ± 16.9 and 73.5 ± 17 years (given the estimated uncertainties that are given by the half-width at half-peak value of the spectral peaks, these values cannot be distinguished robustly). The CrutemV4 component is quite different from the three other curves, with a steep, large-amplitude variation in the last decades of the twentieth century. Figure  shows the fourth HadCrut component, with a pseudoperiod of 22.1 ± 2.4 years. It is strongly modulated and decreases in amplitude from 0.1 °C in 1900 to 0.03 °C in 2000. The other series show very similar components ranging from 20.8 to 22.3 years in pseudoperiod. However, BEST has a double-peaked wider spectral content between 20 and 30 years in period. Figure  shows the sixth HadCrut component, with a pseudoperiod of 11.42 ± 0.47 years. BEST and CrutemV4 show essentially the same peak as HadCrut, whereas it is not found in the first 22 eigenvalues of GISS. The HadCrut and CrutemV4 components are perfectly in phase but the latter has double the amplitude of the former; the BEST component agrees with Hadcrut in amplitude and modulation but is slightly phase lagged. In other terms, this corresponds to a slightly (significantly?) different pseudoperiod (Figure , bottom). The reasons for such phase drift are unclear. One would have to know the details of all steps of data processing in the temperature curves from distinct institutions.

The third component is found close by, at 9.10 ± 0.35 years (Figure ). SSA with this sampling rate and data length allows to separate with significance the 11 from the 9-year component. Hadcrut, GISS, and CrutemV4 show the same spectral peak (to two significant digits) but Hadcrut has a strong modulation from almost nothing in 1850 to 0.1 °C in 2000. The 6.5-year component is found only in the Hadcrut series (Figure ), but a ~5.5-year component is found in CrutemV4 (see section ). The 6.5-year component is modulated and decreases over the 150 years from 0.06 to 0.02 °C in amplitude. On the other end, the three 4.7-year components of HadCrut, GISS, and CrutemV4 are in excellent phase and period agreement but differ in amplitude (Figure ): HadCrut is modulated with a period about 40 to 50 years, whereas CrutemV4 has 2 to 3 times the amplitude of HadCrut and decreases almost monotonously over the century and a half.

The amplitudes of the SSA components of HadCrut are as follows: 0.65 °C (trend, 24% of the total variance), 0.24 °C (68 years, 13% of the total variance), 0.09 °C (22 years, 7% of the total variance), 0.04 °C (11 years, 6% of the total variance), 0.06 °C (9 years, 4% of the total variance), 0.04 °C (6.5 years, 4% of the total variance), and 0.04 °C (4.7 years, 5% of the total variance). SSA allows one to reconstruct the original signal from a reduced number of components. In Figure a we show the reconstruction obtained with only the first two components (that account for 37% of the total variance) and in Figure b the sum of the seven first components shown above (from the trend to the 4.7-year component, amounting to 63% of the total variance). Note again that the four temperature series should represent essentially the same quantity. Hence, the differences between the components of the four series in each of Figures  are evidence of uncertainties in the temperature anomaly curves. For instance, the GISS components are different from the other three and even sometimes missing (Figures , trend and periods of 60, 22, and 11 years) and GISS must refer to a temperature construct somewhat different from the other three. On the other hand, it is reassuring that there is significant agreement between these three.

Figure 10. (top) Reconstruction of the HadCrut series using only the first two eigenvalues/components from the singular spectrum analysis (trend and “60 years”). (bottom) Reconstruction of the HadCrut series using the first seven eigenvalues/components from the singular spectrum analysis.

SSA Analysis of Solar Indices (Sunspot and Faculae Series)

The presence of the 22- and 11-year components found in the temperature anomaly series reveals characteristic time scales of solar activity variations (see Gray et al., , and section below). These periods have already been encountered in our previous SSA analyses mentioned in the introduction (e.g., Le Mouël, Lopes, & Courtillot, ). In this section, we reanalyze briefly the monthly mean sunspot number (ISSN, http://www.sidc.be/silso/datafiles; Figure ) and the monthly mean polar faculae (PF) series (Munoz-Jaramillo et al., ; http://www.solardynamo.org/data.html; Figure ) in the same way we analyzed the temperature series. Figure a shows the eigenvalues for the SSA analysis of sunspot number ISSN from 1850 to 2017, and Figure b shows the same for the faculae series PF from 1906 to 2006. ISSN reveals a trend and successive components at 11, 12.1, 90, 40, 8.5, 14, 8.1, 5.5, 6.2, 22, and 9.1 years. As expected, the 11-year period dominates; the 12.1, 5.5, 6.2 years could be part of the same “11-year related” packets of energy. The 22-year cycle is present and so is the 9.1-year component and the 90-year Gleissberg cycle. As far as PF is concerned, the first component as expected is the solar “magnetic” cycle at 22 years followed by components at 13, 40, 11, 5.5, 8.1, 4.8, and 3.6 years. The components at 40, 22, 11, and 5.5 are common to both ISSN and PF. In PF, components at 13 and 3.6 years could be related to the solar sunspot cycle (the former is statistically close to 11 years and the latter is the second harmonic 11/3). In PF, there are a series of mixed components (eigenvalues of rank 13 to 19, Figure  bottom). As can be checked by changing the window length L, they result from mixing of oscillatory, modulated components at 7.4-, 9.1-, and 33-year periods. With the L value we use in Figure  they cannot be separated.

Figure 11. (top) The ISSN sunspot number series from January 1850 to January 2017. (bottom) Its Fourier spectrum.

Figure 12. (top) The polar faculae (PF) number series from August 1906 to August 2006. (bottom) Its Fourier spectrum.

Figure 13. (top) The first 32 eigenvalues of the sunspot series ISSN (1850–2017) and the corresponding first 12 singular spectrum analysis components. (bottom) Same for the polar faculae series polar faculae (1906–2006).

Comparison of Solar and Temperature SSA Components

We now compare some components of respectively HadCrut and ISSN (and, for some, PF) in Figures .

Figure 14. (top) Comparison of the “60”-year singular spectrum analysis component of the Hadcrut temperature anomaly series (second component) with the corresponding component of the ISSN sunspot series. The latter involves a further singular spectrum analysis of the trend with adapted window length; it is its fourth component. (bottom) Their Fourier transforms.

Figure 15. (top) Comparison of the “22”-year singular spectrum analysis component of the Hadcrut temperature anomaly series (sixth component) with the corresponding components of the ISSN (sunspot, eleventh component) and PF (polar faculae, first component) series. (bottom) Their Fourier transforms.

Figure 16. (top) Comparison of the “11”-year singular spectrum analysis component of the Hadcrut temperature anomaly series (fourth component) with the corresponding components of the ISSN (sunspot, second component) series. (bottom) Their Fourier transforms.

Figure 17. (top) Comparison of the “9”-year singular spectrum analysis component of the Hadcrut temperature anomaly series (third component) with the corresponding component of the ISSN (sunspot, twelfth component) series. (bottom) Their Fourier transforms.

Figure 18. (top) Comparison of the “5.5”-year singular spectrum analysis component of the CruTemV4 temperature anomaly series (fifth component) with the corresponding component of the ISSN series. (bottom) Their Fourier transforms.

Figure 19. (top) Comparison of the “4.8”-year singular spectrum analysis component of the Hadcrut temperature anomaly series (seventh component) with the corresponding component of the PF (polar faculae, seventh component) series. (bottom) Their Fourier transforms.

Figure  shows the “60-year” component. Note that this component does not show as such explicitly in the eigenvalue spectrum of ISSN in Figure  (top), which displays important components at 90 and 40 years. It is actually extracted from a reanalysis of the first component, that is, the trend, again by changing the L value to one that is better adapted to long periods but not to shorter ones. The 60-year component is actually at 57 years for ISSN and 69 years for HadCrut; these could be similar given the large uncertainties (respectively 12 and 23 years). The two components are in phase before 1950 and in phase opposition in the 2000s, ISSN with an amplitude of 0.4, and HadCrut with an amplitude of 0.2 °C. A strong 60-year pattern was identified by Courtillot et al. () with rather abrupt “change points” in the temperature slopes at ~1904, ~1940, and ~1974 (±3 years) and changes at the same dates in Atlantic Multidecadal Oscillation and Pacific Decadal Oscillation. This has been confirmed by Yu and Ruggieri () who detect the location of change points between 1902–1917, 1936–1945, and 1963–1976. Xie et al. () detect peaks in variance of the temperature series at ~70 years and also at ~6.5 and ~4 years, as we do (see below). They find a strong association with the North Atlantic Oscillation, which they use as a predictor of the surface air temperature multidecadal variability.

Patterson et al. () study a sedimentary archive of climate variability, for the interval ~1440–4485 yBP, at annual to millennial scales. Using among other methods wavelet transforms, these authors find cyclicity at 50–85, 33–36, and 22–29 years. They conclude that the sediments have recorded a response to “abrupt changes and long-term variability in climate that can be linked to external forcing (e.g., solar and cosmic irradiance)”. We can associate their “intense, but nonpersistent cyclicity” with our modulated components at “60” and “22” years (and a 30-year harmonic of the “60”-year component). Patterson et al. () also detect multidecadal (Gleissberg) 50–150 years and multicentennial 200- to 500-year variations in climate and solar activity.

Perez-Peraza et al. () give evidence of periodicities around 240, 120, 60, and 30 years in what they believe are cosmic ray fluctuations, based on analysis of cosmogenic isotopes ¹⁰Be and ¹⁴C in ice cores from the South and North Poles. These periodicities are found in several solar and climatic indices. We return to this below.

Figure  compares the “22-year” components of HadCrut, ISSN, and PF. The main peak in both spectra is virtually identical. From 1940 to 2006, the HadCrut component is perfectly in phase with PF. The amplitudes of the two components are respectively 50 for PF and 0.1 °C for HadCrut. The presence of a “22-year” (Hale) component in both the Sun and climate indices is well recognized (e.g., Patterson et al., ). The ISSN spectrum is compatible with PF and HadCrut but slightly offset with respect to both of them. Figure  displays in a similar format the “11-year” components (THE solar cycle) of HadCrut and ISSN, whose spectral peaks overlap and are statistically similar (the slight difference in pseudoperiods (Figure  bottom) can also be seen as a phase drift (Figure 16 top); unfortunately, one would need all information on the complex treatment of data bases by each respective institution). The amplitudes of the two components are respectively 160 for ISSN and 0.04 °C for HadCrut. There is little doubt that the components at 11 and 22 years in the global temperature series are in a way or other due to some form of forcing by the variability in solar activity. The 22-year cycle generates an even larger component than the sunspot cycle.

Figures  and display the 9 and 4.7 components of Hadcrut versus ISSN in the former case and PF in the latter. In Figure  the best agreement is between the 5.25 ± 0.10-year component of Crutem V4 and the “5.5” (actually 5.45 ± 0.30)-year component of ISSN. This 5.5 component is particularly important in a number of geophysical and climate indices (Le Mouël, Lopes, & Courtillot, ; Lopes et al., ). The agreement of the spectral peaks is quite good. The (small) uncertainties in the peak frequency of the sharp spectral peaks are compatible with either the same frequencies for each couple of components or with slight phase drifts between them. But recall that SSA does not have (as Fourier analysis does) to extract orthogonal sinusoidal components but that it adapts to the best shapes, allowing for frequency and amplitude variations (see Le Mouël, Lopes, & Courtillot, ). Note that Xie et al. () do find peaks in variance at the ~6.5- and ~4-year periods.

Discussion

The origin of the 60-year component is not fully clear (Courtillot et al., ). It could be linked to the chaotic coupling between the major oceanic basins (Tsonis et al., ) or to the perturbations of the center of mass of the planetary system (primarily the Sun) under the influence of the giant planets (Scafetta, ; see below). We propose in this paper that the 60-year component of Earth surface temperature variability, the largest one after the trend, is in any case a direct consequence of variations in solar activity.

The origin of the trend in temperature is also a topic of debate. It is often considered as being due to the increase of greenhouse gases released by humanity since the beginning of the industrial age. It could also be the rising part of the multicentennial pseudocycle that caused the Medieval Climatic Optimum followed by the Little Ice Age; the atmosphere has been coming out of the Little Ice Age since the second half of the nineteenth century. It is clear that the period of a multicentennial phenomenon cannot be determined with only 150 years of data. But evidence that the Medieval Climatic Optimum is global and that there were successive occurrences of cooler and warmer intervals on centennial to millennial time scales has been accumulating (e.g., from Dansgaard et al., , to Masson-Delmotte et al., , Moberg et al., , or Wu et al., ). We have recalled above that the Hale, Gleissberg (and we can add de Vries) multidecadal to multicentennial components are all present in both solar and climate indices. As a final test, we have performed a SSA analysis of the more complete series of sunspot numbers ISSN from 1700 to today. The trend is shown in Figure  and compared to the trends of the four temperature series (already seen based on a shorter ISSN series in Figure  top). Although the number of degrees of freedom is small, it is clear that all trends, which can be represented by a polynomial of degree 3, feature the same S-shaped curve from the end of the Little Ice Age to the beginning of the 21st century. This is consistent with a slowdown of the long-term rate of global warming in the first decades of the 21st century. Of course this correlation is suggestive but certainly not a proof of causality by itself. The SSA analysis of the 1700–2018 ISSN series (http://www.sidc.be/silso/datafiles) also displays a prominent ~90-year modulated but quasiperiodic component (not shown in this paper) corresponding to the Gleissberg cycle (e.g., Le Mouël et al., ). Its contribution to the 1850–2018 ISSN series longer-period components cannot be directly evaluated with that shorter record.

Figure 20. Singular spectrum analysis trends of the four temperature series shown in Figure  (top) compared with the trend of the 1700–2018 ISSN series (shown only from 1850 onward). The latter actually corresponds to the addition of the first (trend or de Vries?) and second (Gleissberg) singular spectrum analysis components.

Several authors have attempted to reconstruct the temperature series with a small number of sinusoidal components. A recent analysis tentatively involving the gravitational effects of the giant planets Jupiter and Saturn on the Sun has been put forward by Scafetta (, , ). This author extracts from the global surface temperature records “a warming secular trend plus a set of major oscillations with periods of about 9.1, 10.4, 20, and 60 years.” Scafetta argues for a planetary origin of the 9-, 10- to 12-, 15-, 20- to 30-, and 60-year astronomical harmonics; he does not find the 9.1-year oscillation in the analyzed solar record, but we do, though it is quite strong in the temperature series (Figure ) but weaker in the sunspot series (Figure ). Scafetta () argues that the 9.1-year climatic oscillation is induced by a soli-lunar tidal oscillation (he notes that 9.1 years is about half of the 18.3-year nodal cycle, about half the 18-year Saros eclipse cycle, and about the 8.85-year lunar precession cycle). Future research may determine whether this strong 9.1-year climate oscillation is more related to solar activity variations (as we suggest may be the case) or to a soli-lunar tidal cycle (as Scafetta argues).

Note that whereas Scafetta isolates purely sinusoidal components, using maximum entropy spectral analysis, SSA allows further to identify and isolate the components with their slight frequency changes and their sometimes very significant modulation: Both may carry significant information. His reconstructions (e.g., Scafetta, ; Figures  and ) are good, those based on SSA fit the data at least equally well (Figure ). The fact that SSA is able to extract components with varying amplitudes and pseudoperiods should provide more useful constraints for physical models of solar-terrestrial interactions.

Conclusion

In this paper, we have applied the powerful SSA to long series of monthly mean values of the four most used surface air temperatures (HadCrut, CrutemV4, GISS, and BEST), of sunspot numbers ISSN and PF, from 1850 to 2017 (except for PF available from 1906 to 2006). SSA identifies in the series the trends, pseudoperiodic, and oscillatory components. The Hadley results are used in most comparisons of these SSA results and the three other series are used to show that they converge on several results though not all, which gives an idea of the uncertainties and errors that must still be embedded in these data. We choose to use them at face value and not to fiddle with the details of how these series are built.

The efficiency of the SSA algorithm that we use has been regularly improved (see, e.g., Le Mouël, Lopes, & Courtillot, ). For the temperature series, the SSA eigenvalues are shown in Figure  and the seven first components are shown with their Fourier spectrum in Figures . For the ISSN and PF series, the SSA eigenvalues are shown in Figure . Comparisons of components with similar periods from temperatures on one hand and ISSN or PF on the other hand are shown in Figures . Many of the periods that are found are known to be characteristic of solar activity (Gray et al., ).

The 22- and 11-year components (Figures  and ) are of course the first expected ones in the period range we analyze. These components are often modulated; some of them slowly drift in phase, as reflected by the slight differences in their spectral peak values. It is a strength of SSA that can isolate components that are pseudoperiodic yet rather strongly varying in amplitude and phase. Of course, it is reasonable to associate such components with solar variability.

On the shorter-period side, the common components at 9, 5.5, and 4.7 years (Figures ) are in good to excellent agreement. They have already been identified in solar activity by several authors (see above). Scafetta () suggests a soli-lunar tidal oscillation for the 9.1-year component.

The 60-year component (Figure ) is particularly interesting. It has been recognized in various ways in the temperature series (Figure ; e.g., Scafetta, ; Courtillot et al., ; Yu & Ruggieri, ; Xie et al., ). The 60-year component is not immediately apparent in the solar series (Figures ), though we have seen that it could be extracted with an appropriate choice of the L window value (Figure ). A 60-year pseudoperiod is on the edge of what can reliably be extracted from a 150-year-long time series. But, as we have recalled above, other types of data, such as ice or sedimentary cores, allow one to explore longer periods. Patterson et al. () confirm the existence of climatic variations at ~60, ~35, and ~22 years and find significant variance at 50–150 and 200–500 years. At these periods, the reasonable starting assumption is that they are forced by changes in solar activity. There is abundant literature on these multiple century time scales. Figure  shows that when we consider a longer ISSN time series starting in 1700 and recalculate its SSA first and second components (that are the trend and Gleissberg components) over the time span from 1850 to 2017, the trends of solar activity and temperature are almost identical, with a slower rise before 1900, a faster rise over much of the twentieth century, and a slowdown in the late 1900s (corresponding to what some call the hiatus or plateau). These trends over 150 years can be parts of longer, multicentury changes in solar activity.

It could therefore be argued, as an application of Occam's razor, that most of the variability of surface temperatures could be natural and primarily controlled by the Sun (Figures  and ): Components at 22, 11, 5.5, and 3.6 years must be solar (the latter two being the first two harmonics of 11 years). As we have seen above, the components at 9 and 4.7 years and also 60 years could also be solar generated. The trend itself could, at least in part, be a segment of a much longer, multicentennial solar period. In conclusion, SSA has revealed modulated components that can reasonably be assigned to solar variability for most and even possibly all of them. Naturally, mechanisms should be sought to account for these observations and hypotheses. This is beyond the scope of the present paper. Its main message is to show how interesting it is to actually display the main SSA components of the time series under study (Figures  top), including frequency and amplitude changes (modulation), and to compare components with similar pseudoperiods (Figures  top). When components are shared by two different time series, it is reasonable to assume that there is a common source to these variations, without having to know the full theory of that source.

Our main aim in this paper is not to argue in detail about the characteristics of our implementation of SSA (though, see the appendix). Being geophysicists, we are primarily interested in the results, which is the recognition of similar components in ISSN and temperature series. We believe these results are compelling and exciting. We believe they have not yet been reported with any other spectral method. Would they be in the future, we would see this as a welcome confirmation of our findings.

Appendix

In this appendix, we briefly describe important aspects of the SSA method that we paid particular attention to, in order to improve its efficiency compared to other versions of the method (at the cost of significant computer time and memory). The SSA decomposition of a time series s_N of length N consists in two steps; the first is a decomposition, the second a reconstruction. If all goes well, it allows one to separate the trend, periodical or pseudoperiodical oscillations (amplitude and frequency modulation) and noise in s_N. The first step is the most sensitive, for this is when one builds from s_N the trajectory matrix X that needs to be decomposed. This is when 1-D information contained in s_N is transformed into a 2-D “cross-lagged correlation” matrix X that is diagonal. This operation is called “embedding” (Broomhead & King, , ; Mañé, ; Takens, ). As is done in a sliding Fourier transform, s_N is divided in (possibly overlapping) segments of length L. In the more basic form, the number of lines in X is equal to L, and therefore, X has K columns such that K = N − L − 1 (2 ≤ L ≤ N − 1). L is one of the important “tuning factors” of the method. The smaller L, the closer one comes to a fractal analysis; the larger L, the closer one comes to reconstructing the initial time series. Choosing L is not obvious: It depends on the amount and quality of information carried in s_N and therefore present in X, which itself depends on the conditioning of the matrix. In order to be analyzed, X must be decomposed or more precisely projected on an orthonormal basis. This is what is done in Fourier analysis, where one projects the series on a basis of sinusoidal functions. The same goes for wavelet decomposition. In SSA, the basis is not imposed but is created when X (that depends stricto sensu on L) is decomposed. That is an ad hoc basis for the signal decomposed K times: The signal imposes the basis. Therefore, one must “construct the good matrix for the good signal.” A natural way to obtain that basis is to use singular value decomposition (Golub & Reinsch, ), which gives eigenvalues (that can be assimilated with the coefficients of a Fourier series) and eigenvectors (whose function is that of sine functions in Fourier decomposition). The triplets of eigenvalues and eigenvectors must then be coupled with common sense, in order to properly reconstruct either the trend, or the oscillations, or the noise (see how eigenvalues are grouped in components in Figure ). That step in a basic form of SSA is close to principal component analysis in the frame of multivariate SSA (Golyandina et al., ) or to Karhunen-Loeve decomposition used in the case of stationary signals (Karhunen, ; Loève, ). The main difference between these methods is in the Hankel diagonal matrices used by SSA; the lines and columns have the same meaning in terms of time. Independently from the criteria or method, the grouping of triplets is called (either weak or strong) separability. Let us take a simple example, when the series s_N is the sum of two components s_1N and s_2N, each with trajectory matrix X₁ and X₂. Separability means that X can be decomposed in X₁ and X₂. If there are intersections between the two sets of eigenvalues, separability is weak. Otherwise it is strong, of course the more favorable case. In summary, SSA is sensitive to the nature of the time series s_N (stationary or not, in the strict, wide, ergodic sense, signal to noise ratio,...), to the construction of X (and therefore to the choice of L) and to the separability of eigenvalues (how one groups triplets is a question of its own). Several authors propose simply (and rightly so) to repeat the SSA with different L values: “Whatever the circumstances, it is always a good idea to repeat SSA analysis several times using different values of L” (Golyandina & Zhigljavsky, p 48). This is often called multi-stage SSA or reiterated SSA (Golyandina et al., ). We can also circumscribe a number of problems by optimizing the choice of X, by reducing its conditioning (the ratio of the squares of the first to last eigenvalues must be between 1 and 10; if it is larger, lines and columns can be shifted to reorganize the Hankel matrix into submatrices that obey the conditioning condition and satisfy separability). One can refine information by using statistical inference methods such as the bootstrap (Efron & Tibshirani, ); one builds larger X matrices (ensuring that if L and K are the dimensions of X, then K and L are indeed those of X^T). This implies more powerful computers, time, and memory, but this is becoming possible with present computing power. Also, one can improve separability by implementing matrix rotation as in SSA-ICA (for independant component analysis; Hyvärinen & Oja, ) or varimax rotation (Kaiser, ). These are the tools and solutions we use in our previous papers and the present one. We have in each case to carefully adapt to the data, and this can be different from one series to the other. There is a test phase when we determine the best approaches as a function of the signal that one analyzes. To end with a repeat, one must find the right matrix for the right problem.

Acknowledgments

We thank the Editor and two anonymous reviewers whose questions led to significant improvements. The data used in the paper are accessible at the following sites: the Hadley Research Center and the Climate Research Unit, HadCrut (https://crudata.uea.ac.uk/cru/data/temperature/) and CrutemV4 (https://www.metoffice.gov.uk/hadobs/crutem4/data/download.html), National Aeronautics and Space Administration, GISS (https://data.giss.nasa.gov/gistemp/), Berkeley, BEST (http://berkeleyearth.org/data/), ISSN http://www.sidc.be/silso/datafiles; and PF http://www.solardynamo.org/data.html. IPGP Contribution 4069.