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

J. L. Le Mouël, F. Lopes, V. Courtillot

► To cite this version:

J. L. Le Mouël, F. Lopes, V. Courtillot. Characteristic Time Scales of Decadal to Centennial Changes in Global Surface Temperatures Over the Past 150 Years. Earth and Space Science, 2020, 7, p. 99-122. 10.1029/2019EA000671. insu-03584798

HAL Id: insu-03584798

https://insu.hal.science/insu-03584798

Submitted on 24 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Earth and Space Science

## RESEARCH ARTICLE

10.1029/2019EA000671

This article was commented on in Cuypers et al. (2021), https://doi.org/10.1029/2020EA001298.

There was a reply to the comment in Le Mouel et al. (2021), https://doi.org/10.1029/2020EA001421.

**Key Points:**

* We analyze the four main data sets of global surface temperatures (1850-2017) with the powerful method of singular spectral analysis

* We find spectral periods typical of variations in solar activity

* It can be argued that much of the variability of surface temperatures could be natural and primarily controlled by the Sun

**Correspondence to:**

V. Courtillot,
courtil@ipgp.fr

**Citation:**

Le Mouël, J. L., Lopes, F., & Courtillot, V. (2020). Characteristic time scales of decadal to centennial changes in global surface temperatures over the past 150 years. *Earth and Space Science*, 7, e2019EA000671. https://doi.org/10.1029/2019EA000671

Received 12 APR 2019

Accepted 11 SEP 2019

Accepted article online 14 OCT 2019

©2019. The Authors.
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

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

J. L. Le Mouël¹, F. Lopes¹, and V. Courtillot¹ ID

¹Geomagnetism and Paleomagnetism, Institut de Physique du Globe de Paris, Paris, France

**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.

## 1. 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., 2017), in several climatic indices such as Pacific Decadal Oscillation, Madden-Julian Oscillation, and Atlantic Multidecadal Oscillation (Le Mouël et al., 2019) and in the length of day (Le Mouël et al., 2019). 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:// rudata.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. (2001), and to accounts by Vautard and Ghil (1989), Ghil et al. (2002), and Vautard et al. (1992). 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, 2019, 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).

## 2. The HadCrut Series From 1850 to 2017

Figure 1 (top) shows the full HadCrut data set from January 1850 to January 2017 at a monthly sampling rate (in red). Figure 1 (bottom) shows its Fourier spectrum, with in particular significant energy around 60 and 90 years. Figure 2 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 3–9 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 5).

Figure 3 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 4 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 5 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 6 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 6, 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 7). 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 8), but a ~5.5-year component is found in CrutemV4 (see section 5). 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 9): 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),
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (top) The HadCrut series of global surface mean temperature anomalies from January 1850 to January 2017 (https://rudata.uea.ac.uk/cru/data/temperature/; monthly mean values). (bottom) Its Fourier spectrum.</figcaption></figure>

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 10a we show the reconstruction obtained with only the first two components (that account for 37% of the total variance) and in Figure 10b the sum of the seven first
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The first 32 singular spectrum analysis eigenvalues and seven components of series HadCrut (see Figure 1)</figcaption></figure>

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 3–9 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 3–6, 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.

**3. 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., 2010, and section 5 below). These periods have already been encountered in our previous SSA analyses mentioned in the introduction (e.g., Le Mouël, Lopes, & Courtillot, 2019). In this section, we reanalyze briefly the monthly mean sunspot number (ISSN, http://www sidc be/silso/datafiles; Figure 11) and the monthly mean polar faculae (PF) series (Munoz-Jaramillo et al., 2012; http://www solardynamo.org/data.html; Figure 12) in the same way we analyzed the temperature series. Figure 13a shows the eigenvalues for the SSA analysis of sunspot number ISSN from 1850 to 2017, and Figure 13b 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
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (top). The first singular spectrum analysis component of series HadCrut (see Figure 1), CrutemV4, GISS, and BEST, that is, their trends. (bottom) Their Fourier transforms.</figcaption></figure>

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
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (top). The second singular spectrum analysis component of series HadCrut (see Figure 1), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.</figcaption></figure>

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 13 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 13 they cannot be separated.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (top) The fourth singular spectrum analysis component of series HadCrut (see Figure 1), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.</figcaption></figure>

**4. Comparison of Solar and Temperature SSA Components**

We now compare some components of respectively HadCrut and ISSN (and, for some, PF) in Figures 14–19.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (top) The sixth singular spectrum analysis component of series HadCrut (see Figure 1), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.</figcaption></figure>

Figure 14 shows the “60-year” component. Note that this component does not show as such explicitly in the eigenvalue spectrum of ISSN in Figure 13 (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
<figure><img src="image_9.png" /><figcaption>**Figure 7.** (top) The third singular spectrum analysis component of series HadCrut (see Figure 1), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.</figcaption></figure>

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. (2013) with rather abrupt “change points” in the temperature slopes at ~1904, ~1940, and
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (top) The fifth singular spectrum analysis component of series HadCrut (see Figure 1). (bottom) Its Fourier transform.</figcaption></figure>

~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 (2019) who detect the location of change points between 1902–1917, 1936–1945, and 1963–1976. Xie et al. (2019) 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
<figure><img src="image_11.png" /><figcaption>**Figure 9.** (top) The seventh singular spectrum analysis component of series HadCrut (see Figure 1), CrutemV4, GISS, and BEST. (bottom) Their Fourier transforms.</figcaption></figure>

association with the North Atlantic Oscillation, which they use as a predictor of the surface air temperature multidecadal variability.

Patterson et al. (2004) 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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

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.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** (top) The ISSN sunspot number series from January 1850 to January 2017. (bottom) Its Fourier spectrum.</figcaption></figure>

(2004) also detect multidecadal (Gleissberg) 50–150 years and multicentennial 200- to 500-year variations in climate and solar activity.

Perez-Peraza et al. (2012) 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 $^{10}$Be and $^{14}$C in ice cores from the
<figure><img src="image_14.png" /><figcaption>**Figure 12.** (top) The polar faculae (PF) number series from August 1906 to August 2006. (bottom) Its Fourier spectrum.</figcaption></figure>

South and North Poles. These periodicities are found in several solar and climatic indices. We return to this below.

Figure 15 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
<figure><img src="image_15.png" /><figcaption>**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).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**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.</figcaption></figure>

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., 2004). The ISSN spectrum is compatible with PF and HadCrut but slightly offset with respect to both of them. Figure 16 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 16 bottom) can also be seen as a phase drift (Figure 16 top)); unfortunately, one would need all information
<figure><img src="image_17.png" /><figcaption>**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.</figcaption></figure>

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.
<figure><img src="image_18.png" /><figcaption>**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.</figcaption></figure>

Figures 17 and 19 display the 9 and 4.7 components of Hadcrut versus ISSN in the former case and PF in the latter. In Figure 18 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, 2019; Lopes et al., 2017). The
<figure><img src="image_19.png" /><figcaption>**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.</figcaption></figure>

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
<figure><img src="image_20.png" /><figcaption>**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.</figcaption></figure>

(see Le Mouël, Lopes, & Courtillot, 2019). Note that Xie et al. (2019) do find peaks in variance at the ~6.5- and ~4-year periods.

**5. Discussion**

The origin of the 60-year component is not fully clear (Courtillot et al., 2013). It could be linked to the chaotic coupling between the major oceanic basins (Tsonis et al., 2007) or to the perturbations of the center of mass
<figure><img src="image_21.png" /><figcaption>**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.</figcaption></figure>

of the planetary system (primarily the Sun) under the influence of the giant planets (Scafetta, 2010; 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
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Singular spectrum analysis trends of the four temperature series shown in Figure 3 (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.</figcaption></figure>

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., 1969, to Masson-Delmotte et al., 2004, Moberg et al., 2005, or Wu et al., 2018). 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 20 and compared to the trends of the four temperature series (already seen based on a shorter ISSN series in Figure 3 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., 2017). Its contribution to the 1850–2018 ISSN series longer-period components cannot be directly evaluated with that shorter record.

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 (2010, 2012, 2016). 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 2) but weaker in the sunspot
series (Figure 13). Scafetta (2012) 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, 2016; Figures 12 and 13) are good, those based on SSA fit the data at least equally well (Figure 10). 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.

**6. 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, 2019). For the temperature series, the SSA eigenvalues are shown in Figure 2 and the seven first components are shown with their Fourier spectrum in Figures 3–9. For the ISSN and PF series, the SSA eigenvalues are shown in Figure 13. Comparisons of components with similar periods from temperatures on one hand and ISSN or PF on the other hand are shown in Figures 14–19. Many of the periods that are found are known to be characteristic of solar activity (Gray et al., 2010).

The 22- and 11-year components (Figures 15 and 16) 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 17–19) are in good to excellent agreement. They have already been identified in solar activity by several authors (see above). Scafetta (2012) suggests a soli-lunar tidal oscillation for the 9.1-year component.

The 60-year component (Figure 14) is particularly interesting. It has been recognized in various ways in the temperature series (Figure 2; e.g., Scafetta, 2010; Courtillot et al., 2013; Yu & Ruggieri, 2019; Xie et al., 2019). The 60-year component is not immediately apparent in the solar series (Figures 11–13), though we have seen that it could be extracted with an appropriate choice of the *L* window value (Figure 14). 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. (2004) 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 20 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 2 and 10): 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 3–9 top), including frequency and amplitude changes (modulation), and to compare components with similar pseudoperiods (Figures 14–20 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, 1986a, 1986b; Mañé, 1981; Takens, 1981). 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 \le L \le 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, 1971), 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 2). That step in a basic form of SSA is close to principal component analysis in the frame of multivariate SSA (Golyandina et al., 2013) or to Karhunen-Loeve decomposition used in the case of stationary signals (Karhunen, 1947; Loève, 1978). 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 $\mathbf{X}_1$ and $\mathbf{X}_2$. Separability means that **X** can be decomposed in $\mathbf{X}_1$ and $\mathbf{X}_2$. 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,
2013 p 48). This is often called multi-stage SSA or reiterated SSA (Golyandina et al., 2001). 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, 1997); 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 independent component analysis; Hyvärinen & Oja, 2000) or varimax rotation (Kaiser, 1958). 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.

References

Broomhead, D. S., & King, G. P. (1986a). Extracting qualitative dynamics from experimental data. *Physica D, 20*, 217–236.

Broomhead, D. S., & King, G. P. (1986b). On the qualitative analysis of experimental dynamical systems. In S. Sarkar (Ed.), *Nonlinear Phenomena and Chaos* (pp. 113–144). Bristol: Adam Hilger.

Courtillot, V., Le Mouël, J. L., Kossobokov, V., Gibert, D., & Lopes, F. (2013). Multi-decadal trends of global surface temperature: A broken line with alternating ~30 yr linear segments? *Atmospheric and Climate Sciences, 3*, 364–371. https://doi.org/10.4236/acs.2013.33038

Dansgaard, W., Johnsen, S. J., Möller, J., & Langway, C. C. (1969). One thousand centuries of climatic record from Camp Century on the Greenland Ice Sheet. *Science, 166*, 377–381.

Efron, B., & Tibshirani, R. (1997). Improvements on cross-validation: The 632+ bootstrap method. *Journal of the American Statistical Association, 92*(438), 548–560.

Ghil, M., Allen, M. R., Dettinger, M. D., Ide, K., Kondrashov, D., Mann, Michael, et al. (2002). Advanced spectral methods for climatic time series. *Reviews of Geophysics, 40*(1), 1003. https://doi.org/10.1029/RG000092

Golub, G. H., & Reinsch, C. (1971). *Singular value decomposition and least squares solutions, Linear Algebra*. Springer, Berlin, Heidelberg: Springer, Linear Algebra (pp.134–151).

Golyandina, N., Korobeynikov, A., Shlemov, A., & Usevich, K. (2013). Multivariate and 2D extensions of singular spectrum analysis with the Rssa package, arXiv.

Golyandina, N., Nekrutkin, V., & Zhigljasky, A. (2001). *Analysis of time series structure: SSA and related techniques*. Monographs on Statistics and Applied Probability, (Vol. 90). Boca Raton, 295pp: Chapman and Hall. https://doi.org/10.1201/9781420035841

Golyandina, N., & Zhigljavsky, A. (2013). *Singular spectrum analysis for time series*. Heidelberg, New York, Dordrecht, London: Springer Briefs in Statistics 119pp. https://doi.org/10.1007/978-3-642-34913-3

Gray, L. J., Beer, J., Geller, M., Haigh, J. D., Lockwood, M., Matthes, K., et al. (2010). Solar influences on climate. *Reviews of Geophysics, 48*, RG4001. https://doi.org/10.1029/2009RG000282

Hyvärinen, A., & Oja, E. (2000). Independent component analysis: algorithms and applications. *Neural Networks, 13*(4-5), 411–430. https://doi.org/10.1016/S0893-6080(00)00026-5

Kaiser, H. (1958). The varimax criterion for analytic rotation in factor analysis. *Psychometrika, 23*(3), 187–200. https://doi.org/10.1007/BF02289233

Karhunen, K. (1947). Über lineare Methoden in der Wahrscheinlichkeitsrechnung. *Annales Academiæ Scientiarum Fennicæ Mathematica, 37*, 1–79.

Le Mouël, J. L., Lopes, F., & Courtillot, V. (2017). Identification of Gleissberg cycles and a rising trend in a 315-year-long series of sunspot numbers. *Solar Physics, 292*(3), 43. https://doi.org/10.1007/s11207-017-1067-6

Le Mouël, J. L., Lopes, F., & Courtillot, V. (2019). A Solar signature in several climatic indices. *Journal of Geophysical Research: Atmospheres, 124*(5), 2600–2619. https://doi.org/10.1029/2018JD028939

Le Mouël, J. L., Lopes, F., Courtillot, V., & Gibert, D. (2019). On Forcings of Length of Day Changes: From 9-day to 18.6-year oscillations. *Physics of the Earth and Planetary Interiors, 292*, 1–11. https://doi.org/10.1016/j.pepi.2019.04.006

Loève, M. (1978). Probability theory, Vol. II, 4th ed. Graduate Texts in Mathematics, 46 Springer-Verlag, ISBN 978-0-387-90262-3.

Lopes, F., Le Mouël, J. L., & Gibert, D. (2017). The mantle rotation pole position: A solar component. *Comptes Rendus Geoscience, 349*(4), 159–164. https://doi.org/10.1016/j.crte.2017.06.001

Mañé, R. (1981). On the dimension of the compact invariant sets of certain nonlinear maps. In D. A. Rand, & L. S. Young (Eds.), *Dynamical systems and turbulence*, Lecture Notes in Mathematics, (Vol. 898, pp. 230–242). Springer, Berlin, Heidelberg: Springer-Verlag.

Masson-Delmotte, V., Stenni, B., & Jouzel, J. (2004). Common millennial-scale variability of Antarctic and Southern Ocean temperatures during the past 5000 years reconstructed from the EPICA Dome C ice core. *The Holocene, 14*(2), 145–151. https://doi.org/10.1191/0959683604h1697ft

Moberg, A., Sonechkin, D. M., Holmgren, K., Datsenko, N. M., & Karlen, W. (2005). Highly variable Northern Hemisphere temperatures reconstructed from low- and high-resolution proxy data. *Nature, 433*(7026), 613–617. https://doi.org/10.1038/nature03265

Munoz-Jaramillo, A., Sheeley, N., Zhang, J., & DeLuca, E. (2012). Calibrating 100 years of polar faculae measurements: Implications for the evolution of the heliospheric magnetic field. *The Astrophysical Journal, 753*(2), 146. https://doi.org/10.1088/0004-637X/753/2/146

Patterson, R. T., Prokoph, A., & Chang, A. (2004). Late Holocene sedimentary response to solar and cosmic ray activity influenced climate variability in the NE Pacific. *Sedimentary Geology, 172*(1-2), 67–84. https://doi.org/10.1016/j.sedgeo.2004.07.007

Perez-Peraza, J., Velasco, V., Libin, I. Y., & Yudakhin, K. E. (2012). Thirty-year periodicity of cosmic rays. *Advances in Astronomy, 2012*, 1–11. https://doi.org/10.1155/2012/691408
Scafetta, N. (2010). Empirical evidence for a celectial origin of the climate oscillations and its implications. *Journal of Atmospheric and Solar - Terrestrial Physics*, 72(13), 951-970. https://doi.org/10.1016/j.jastp.2010.04.015

Scafetta, N. (2012). Testing an astronomically based decadal-scale empirical harmonic climate model versus the IPCC (2007) general circulation climate models. *Journal of Atmospheric and Solar - Terrestrial Physics*, 80, 124-137. https://doi.org/10.1016/j.jastp.2011.12.005

Scafetta, N. (2016). Problems in modeling and forecasting climate change: CMIP5 general circulation models versus a semi-empirical model based on natural oscillations. *International Journal of Heat and Technology*, 34, Spec. Iss. 34(S2), S435-S442. https://doi.org/10.18280/ijht.34S235

Takens, F. (1981). Detecting strange attractors in turbulence. In D. A. Rand, & L. S. Young (Eds.), *Dynamical systems and turbulence*, Lecture Notes in Mathematics (Vol. 898, pp. 366-381). Springer, Berlin, Heidelberg: Springer-Verlag.

Tsonis, A. A., Swanson, K., & Kravtsov, S. (2007). A new dynamical mechanism for major climate shifts. *Geophysical Research Letters*, 34, L13705. https://doi.org/10.1029/2007GL030288

Vautard, R., & Ghil, M. (1989). Singular spectrum analysis in nonlinear dynamics, with applications to paleoclimatic time series. *Physica D: Nonlinear Phenomena*, 35(3), 395-424. https://doi.org/10.1016/0167-2789(89)90077-8

Vautard, R., Yiou, P., & Ghil, M. (1992). Singular-spectrum analysis: A toolkit for short, noisy chaotic signals. *Physica D: Nonlinear Phenomena*, 58(1-4), 95-126. https://doi.org/10.1016/0167-2789(92)90103-T

Wu, C. J., Krivova, N. A., Solanki, S. A., & Usoskin, I. G. (2018). Solar total and spectral irradiance reconstruction over the last 9000 years. *Astronomy and Astrophysics*, 620, A120. https://doi.org/10.1051/0004-6361/201832956

Xie, Tiejun, Li, Jianping, Sun, Cheng, Ding, Ruiqiang, Wang, Kaicun, Zhao, Chuanfeng, & Feng, Juan (2019). NAO implicated as a predictor of the surface air temperature multidecadal variability over East Asia. *Climate Dynamics*, 53(1-2), 895-905. https://doi.org/10.1007/s00382-019-04624-4

Yu, M., & Ruggieri, E. (2019). Change point analysis of global temperature records. *International Journal of Climatology*, 39(8), 3679-3688. https://doi.org/10.1002/joc.6042