How much has the Sun influenced Northern Hemisphere temperature trends? An ongoing debate

Because most of the energy that keeps the Earth warmer than space comes from incoming solar radiation, i.e., TSI, it stands to reason that a multidecadal increase in TSI should cause global warming (all else being equal). Similarly, a multidecadal decrease in TSI should cause global cooling. For this reason, for centuries (and longer), researchers have speculated that changes in solar activity could be a major driver of climate change (Laut 2003; Gray et al. 2010; Lockwood 2012; Hoyt & Schatten 1997; Singer & Avery 2008; Soon et al. 2015; Maunder & Maunder 1908; Soon & Yaskell 2003; Scafetta 2010, 2014a). However, a challenging question associated with this theory is "How exactly has TSI changed over time?"

One indirect metric on which much research has focused is the examination of historical records of the numbers and types/sizes of "sunspots" that are observed on the Sun's surface over time (Beck et al. 2014; Hoppe & Rödder 2019; Sarewitz 2011; Hulme 2013; Bateman et al. 2005; Kahneman & Klein 2009; Rakow et al. 2015; Matthes et al. 2017). Sunspots are intermittent magnetic phenomena associated with the Sun's photosphere, that appear as dark blotches or blemishes on the Sun's surface when the light from the Sun is shone on a card with a telescope (to avoid the observer directly looking at the Sun). These have been observed since the earliest telescopes were invented, and Galileo Galilei and others were recording sunspots as far back as 1610 (Soon & Yaskell 2003; Hoyt & Schatten 1998; Svalgaard & Schatten 2016; Vaquero et al. 2016; Schove 1955; Usoskin et al. 2015). The Chinese even have intermittent written records since 165 B.C. of sunspots that were large enough to be seen by the naked eye (Wang & Li 2019, 2020) Moreover, an examination of the sunspot records reveals significant changes on sub-decadal to multidecadal timescales. In particular, a pronounced "sunspot cycle" exists over which the number of sunspots rises from zero during the Sunspot Minimum to a Sunspot Maximum where many sunspots occur, before decreasing again to the next Sunspot Minimum. The length of this "sunspot cycle" or "solar cycle" is typically about 11 years, but it can vary between 8 and 14 years. This 11-year cycle in sunspot behavior is part of a 22-year cycle in magnetic behavior known as the Hale Cycle. Additionally, multidecadal and even centennial trends are observed in the sunspot numbers. During the period from 1645 to 1715, known as the "Maunder Minimum" (Soon & Yaskell 2003; Hoyt & Schatten 1998; Svalgaard & Schatten 2016; Vaquero et al. 2016; Usoskin et al. 2015), sunspots were very rarely observed at all.

Clearly, these changes in sunspot activity are capturing some aspect of solar activity, and provide evidence that the Sun is not a constant star, but one whose activity shows significant variability on short and long timescales. Therefore, the sunspot records initially seem like an exciting source of information on changes in solar activity. However, as will be discussed in more detail later, it is still unclear how much of the variability in TSI is captured by the sunspot numbers. The fact that sunspot numbers are not the only important measure of solar activity (as many researchers often implicitly assume, e.g., Gil-Alana et al. (2014)) can be recognized by the simple realization that TSI does not fall to zero every ∼ 11 years during sunspot minima, even though the sunspot numbers do. Indeed, satellite measurements confirm that sunspots actually reduce solar luminosity, yet paradoxically the average TSI increases during sunspot maxima and decreases during sunspot minima (Willson & Hudson 1988; Lean & Foukal 1988; Foukal & Lean 1990; Willson & Hudson 1991).

We will discuss the current explanations for the apparently paradoxical relationship between sunspots and TSI in Sections 2.2 and 2.3. In any case, the fact that there is more to solar activity than sunspot numbers was recognized more than a century ago by Maunder & Maunder (1908) (for whom the "Maunder Minimum" is named) who wrote,

"...for sun-spots are but one symptom of the sun's activity, and, perhaps, not even the most important symptom" – Maunder & Maunder (1908), pp189-190 (Maunder & Maunder 1908);

and,

"A ['great'] spot like that of February, 1892 is enormous of itself, but it is a very small object compared to the sun; and spots of such size do not occur frequently, and last but a very short time. We have no right to expect, therefore, that a time of many sun-spots should mean any appreciable falling off in the light and heat we have from the sun. Indeed, since the surface round the spots is generally bright beyond ordinary, it may well be that a time of many spots means no falling off, but rather the reverse." – Maunder & Maunder (1908), p183 (Maunder & Maunder 1908)

At the start of the 20^th century, Langley, Abbott and others at the Smithsonian Astrophysical Observatory (SAO) recognized that a more direct estimate of the variability in TSI was needed (Langley 1904; Abbot 1911; Hoyt 1979a). From 1902 until 1962, they carried out a fairly continuous series of measurements of the "solar constant", i.e., the average rate per unit area at which energy is received at the Earth's average distance from the Sun, i.e., 1 Astronomical Unit (AU). The fact they explicitly considered the solar constant to understand climate change is apparent from the title of one of the first papers describing this project, Langley (1904), i.e., "On a possible variation of the solar radiation and its probable effect on terrestrial temperatures". However, they were also acutely aware of the inherent challenges in trying to estimate changes in solar radiation from the Earth's surface:

"The determination of the solar radiation towards the Earth, as it might be measured outside the Earth's atmosphere (called the "solar constant"), would be a comparatively easy task were it not for the almost insuperable difficulties introduced by the actual existence of such an atmosphere, above which we cannot rise, though we may attempt to calculate what would be the result if we could." (Langley 1904).

The true extent of this problem of estimating the changes in TSI from beneath the atmosphere became apparent later in the program. Initially, by comparing the first few years of data, it looked like changes in TSI of the order of 10% were occurring. However, it was later realized that, coincidentally, major (stratosphere-reaching) volcanic eruptions occurred near the start of the program: at Mt. Pelée and La Soufrière (1902) and Santa Maria (1903). Hence, the resulting stratospheric dust and aerosols from these eruptions had temporarily reduced the transmission of solar radiation through the atmosphere (Hoyt 1979a).

It was not until much later in the 20^th century that researchers overcame this ground-based limitation through the use of rocket-borne (Johnson 1954), balloon-borne (Kosters & Murcray 1979) and spacecraft measurements (Foukal et al. 1977). Ultimately, when Hoyt (1979a) systematically reviewed the entire ∼60 year-long SAO solar constant project, he found unfortunately that any potential trend in the solar constant over the record was probably less than the accuracy of the measurements (∼0.3%) (Hoyt 1979a). However, with the launch of the Nimbus 7 Earth Radiation Budget (ERB) satellite mission in 1978 and the Solar Maximum Mission (SMM) Active Cavity Radiometer Irradiance Monitor 1 (ACRIM1) satellite mission in 1980, it finally became possible to continuously and systematically monitor the incoming TSI for long periods from above the Earth's atmosphere (Willson & Hudson 1988, 1991; Hoyt et al. 1992; Willson et al. 1981).

Although each satellite mission typically provides TSI data for only 10 to 15 years, and the data can be affected by gradual long-term orbital drifts and/or instrumental errors that can be hard to identify and quantify (BenMoussa et al. 2013), there has been an almost continuous series of TSI-monitoring satellite missions since those two initial US missions, including European missions, e.g., SOVAP/Picard (Meftah et al. 2014) and Chinese missions (Fang et al. 2014; Wang et al. 2017) as well as international collaborations, e.g., VIRGO/SOHO (Fröhlich et al. 1997), and further US missions, e.g., ACRIMSAT/ACRIM3 (Willson 2014) and SORCE/TIM (Kopp 2016). Therefore, in principle, by rescaling the measurements from different parallel missions so that they have the same values during the periods of overlap, it is possible to construct a continuous time series of TSI from the late-1970s to the present.

Therefore, it might seem reasonable to assume that we should at least have a fairly reliable and objective understanding of the changes in TSI during the satellite era, i.e., 1978 to present. However, even within the satellite era, there is considerable ongoing controversy over what exactly the trends in TSI have been (Scafetta 2011; Scafetta & Willson 2014; Soon et al. 2015; Scafetta et al. 2019; Beer et al. 2000; Dudok de Wit et al. 2017; Fröhlich 2012; Gueymard 2018). There are a number of rival composite datasets, each implying different trends in TSI since the late-1970s. All composites agree that TSI exhibits a roughly 11-year cycle that matches well with the sunspot cycle discussed earlier. However, the composites differ in whether additional multidecadal trends are occurring.

The composite of the ACRIM group that was in charge of the three ACRIM satellite missions (ACRIM1, ACRIM2 and ACRIM3) suggests that TSI generally increased during the 1980s and 1990s but has slightly declined since then (Scafetta & Willson 2014; Scafetta et al. 2019; Willson 2014; Scafetta & Willson 2019). The Royal Meteorological Institute of Belgium (RMIB)'s composite implies that, aside from the sunspot cycle, TSI has remained fairly constant since at least the 1980s (Dewitte & Nevens 2016). Meanwhile, the Physikalisch- Meteorologisches Observatorium Davos (PMOD) composite implies that TSI has been steadily decreasing since at least the late-1970s (Fröhlich 2012, 2009). Additional TSI satellite composites have been produced by Scafetta (2011); Dudok de Wit et al. (2017) and Gueymard (2018).

The two main rival TSI satellite composites are ACRIM and PMOD. As we will discuss in Section 3, global temperatures steadily increased during the 1980s and 1990s but seemed to slow down since the end of the 20^th century. Therefore, the debate over these three rival TSI datasets for the satellite era is quite important. If the ACRIM dataset is correct, then it suggests that much of the global temperature trends during the satellite era could have been due to changes in TSI (Willson & Mordvinov 2003; Scafetta & West 2008b; Scafetta 2009, 2011; Scafetta & Willson 2014; Scafetta et al. 2019; Willson 2014; Scafetta & Willson 2019). However, if the PMOD dataset is correct, and we assume for simplicity a linear relationship between TSI and global temperatures, then the implied global temperature trends from changes in TSI would exhibit long-term global cooling since at least the late-1970s. Therefore, the PMOD dataset implies that none of the observed warming since the late-1970s could be due to solar variability, and that the warming must be due to other factors, e.g., increasing greenhouse gas concentrations. Moreover, it implies that the changes in TSI have been partially reducing the warming that would have otherwise occurred; if this TSI trend reverses in later decades, it might accelerate "global warming" (Fröhlich 2009; Fröhlich & Lean 2002).

The PMOD dataset is more politically advantageous to justify the ongoing considerable political and social efforts to reduce greenhouse gas emissions under the assumption that the observed global warming since the late-19^th century is mostly due to greenhouse gases. Indeed, as discussed in Soon et al. (2015), Dr. Judith Lean (of the PMOD group) acknowledged in a 2003 interview that this was one of the motivations for the PMOD group to develop a rival dataset to the ACRIM one by stating,

"The fact that some people could use Willson's [ACRIM dataset] results as an excuse to do nothing about greenhouse gas emissions is one reason we felt we needed to look at the data ourselves" – Dr. Judith Lean, interview for NASA Earth Observatory, August 2003 (Lindsey 2003)

Similarly, Zacharias (2014) argued that it was politically important to rule out the possibility of a solar role for any recent global warming,

"A conclusive TSI time series is not only desirable from the perspective of the scientific community, but also when considering the rising interest of the public in questions related to climate change issues, thus preventing climate skeptics from taking advantage of these discrepancies within the TSI community by, e.g., putting forth a presumed solar effect as an excuse for inaction on anthropogenic warming." – Zacharias (2014)

We appreciate that some readers may share the sentiments of Lean and Zacharias and others and may be tempted to use these political arguments for helping them to decide their opinion on this ongoing scientific debate. In this context, readers will find plenty of articles to use as apparent scientific justification, e.g., Refs. (Lean 2017; Meftah et al. 2014; Dudok de Wit et al. 2017; Fröhlich 2012; Dewitte & Nevens 2016; Fröhlich 2009; Fröhlich & Lean 2002; Zacharias 2014; Kopp et al. 2016; Lean 2018). It may also be worth noting that the IPCC appears to have taken the side of the PMOD group in their most recent AR5 - see section 8.4.1 of IPCC (2013a) for the key discussions. However, we would encourage all readers to carefully consider the counter-arguments offered by the ACRIM group, e.g., Refs. (Willson & Mordvinov 2003; Scafetta & Willson 2014; Scafetta et al. 2019; Willson 2014; Scafetta & Willson 2019). In our opinion, this was not satisfactorily done by the authors of the relevant section in the influential IPCC reports, i.e., section 8.4.1 of IPCC (2013a). Matthes et al. (2017)'s recommendation that their new estimate (which will be discussed below) should be the only solar activity dataset considered by the CMIP6 modeling groups (Matthes et al. 2017) for the IPCC's upcoming AR6 is even more unwise due to the substantial differences between various published TSI estimates. This is aside from the fact that Scafetta et al. (2019) have argued that the TSI proxy reconstructions preferred in Matthes et al. (2017) (i.e., NRLTSI2 and SATIRE) contradict important features observed in the ACRIM 1 and ACRIM 2 satellite measurements. We would also encourage readers to carefully read the further discussion of this debate in Soon et al. (2015).

The debate over which satellite composite is most accurate also has implications for assessing TSI trends in the pre-satellite era. In particular, there is ongoing debate over how closely the variability in TSI corresponds to the variability in the sunspot records. This is important, because if the match is very close, then it implies the sunspot records can be a reliable solar proxy for the pre-satellite era (after suitable scaling and calibration has been carried out), but if not other solar proxies may need to be considered.

In the 1980s and early 1990s, data from the NIMBUS7/ERB and ACRIMSAT/ACRIM1 satellite missions suggested a cyclical component to the TSI variability that was highly correlated to the sunspot cycle. That is, when sunspot numbers increased, so did TSI, and when sunspot numbers decreased, so did TSI (Willson & Hudson 1988; Lean & Foukal 1988; Foukal & Lean 1990; Willson & Hudson 1991; Hoyt et al. 1992; Wade 1995; Willson 1997). This was not known in advance, and it was also unintuitive because sunspots are "darker", and so it might be expected that more sunspots would make the Sun "less bright" and therefore lead to a lower TSI. Indeed, the first six months of data from the ACRIM1 satellite mission suggested that this might be the case because "two large decreases in irradiance of up to 0.2 percent lasting about 1 week [were] highly correlated with the development of sunspot groups" – Willson et al. (1981). However, coincidentally, it appears that the sunspot cycle is also highly correlated with changes in the number of "faculae" and in the "magnetic network", which are different types of intermittent magnetic phenomena that are also associated with the Sun's photosphere, except that these phenomena appear as "bright" spots and features. It is now recognized that the Sun is currently a "faculae-dominated star". That is, even though sunspots themselves seem to reduce TSI, when sunspot numbers increase, the number of faculae and other bright features also tend to increase, increasing TSI, and so the net result is an increase in TSI. That is, the increase in brightness from the faculae outweighs the decrease from sunspot dimming (i.e., the faculae:sunspot ratio of contributions to TSI is greater than 1). For younger and more active stars, the relative contribution is believed to be usually reversed (i.e., the ratio is less than 1) with the changes in stellar irradiance being "spot-dominated" (Lockwood et al. 2007; Hall et al. 2009; Reinhold et al. 2019).

At any rate, it is now well-established that TSI slightly increases and decreases over the sunspot cycle in tandem with the rise and fall in sunspots (which coincides with a roughly parallel rise and fall in faculae and magnetic network features) (Willson & Hudson 1988; Lean & Foukal 1988; Foukal & Lean 1990; Willson & Hudson 1991; Hoyt et al. 1992; Wade 1995; Willson 1997). Many of the current pre-satellite era TSI reconstructions are based on this observation. That is, a common approach to estimating past TSI trends includes the following three steps:

• 1  
  Estimate a function to describe the inter-relationships between sunspots, faculae and TSI during the satellite era.
• 2  
  Assume these relationships remained reasonably constant over the last few centuries at least.
• 3  
  Apply these relationships to one or more of the sunspot datasets and thereby extend the TSI reconstruction back to 1874 (for sunspot areas (Foukal & Lean 1990; Wang et al. 2005)); 1700 (for sunspot numbers (Clette et al. 2014; Clette & Lefèvre 2016); or 1610 (for group sunspot numbers (Hoyt & Schatten 1998; Svalgaard & Schatten 2016)).

Although there are sometimes additional calculations and/or short-term solar proxies involved, this is the basic approach adopted by, e.g., Foukal & Lean (1990); Lean (2000); Solanki et al. (2000, 2002); Wang et al. (2005); Krivova et al. (2007, 2010). Soon et al. (2015) noted that this heavy reliance on the sunspot datasets seems to be a key reason for the similarities between many of the TSI reconstructions published in the literature (Soon et al. 2015).

However, does the relationship between faculae, sunspots and TSI remain fairly constant over multidecadal and even centennial timescales? Also, would the so-called "quiet" solar region remain perfectly constant despite multidecadal and secular variability observed in the sunspot and faculae cycles? Is it reasonable to assume that there are no other aspects of solar activity that contribute to variability in TSI? If the answers to all these questions are yes, then we could use the sunspot record as a proxy for TSI, scale it accordingly and extend the satellite record back to the 17^th century. This would make things much simpler. It would mean that, effectively, even Galileo Galilei could have been able to determine almost as much about the changes in TSI of his time with his early 16^th-century telescope as a modern (very high budget) Sun-monitoring satellite mission of today. All that he would have been missing was the appropriate scaling functions to apply to the sunspot numbers to determine TSI.

If the PMOD or similar satellite composites are correct, then it does seem that, at least for the satellite era (1978-present), the sunspot cycle is the main variability in TSI and that the relationships between faculae, sunspots and TSI have remained fairly constant. This is because the trends of the PMOD composite are highly correlated to the trends in sunspot numbers over the entire satellite record. However, while the ACRIM composite also has a component that is highly correlated to the sunspot cycle (and the faculae cycle), it implies that there are also additional multidecadal trends in the solar luminosity that are not captured by a linear relation between sunspot and faculae records. Recent modeling by Rempel (2020) is consistent with this in that his analysis suggests even a 10% change in the quiet-Sun field strength between solar cycles could lead to an additional TSI variation comparable in magnitude to that over a solar cycle (Rempel (2020)). Therefore, if the ACRIM composite is correct, then it would be necessary to consider additional proxies of solar activity that are capable of capturing these non-sunspot number-related multidecadal trends.

Over the years, several researchers have identified several time series from the records of solar observers that seem to be capturing different aspects of solar variability than the basic sunspot numbers (Hoyt & Schatten 1997, 1993; Livingston 1994; Hoyt 1979b; Friis-Christensen & Lassen 1991; Solanki & Fligge 1998; Nesme-Ribes et al. 1993; Lean et al. 1995). Examples include the average umbral/penumbral ratio of sunspots (Hoyt 1979b), the length of sunspot cycles (Solheim et al. 2012; Beer et al. 2000; Friis-Christensen & Lassen 1991; Solanki & Fligge 1998; Lassen & Friis-Christensen 1995; Soon et al. 1994; Butler 1994; Zhou et al. 1998), solar rotation rates (Nesme-Ribes et al. 1993), the "envelope" of sunspot numbers (Reid 1991; Lean et al. 1995), variability in the 10.7 cm solar microwave emissions (Labitzke & van Loon 1988; Foukal 1998a), solar plage areas (e.g., from Ca II K spectroheliograms) (Foukal 1998a,b; Foukal et al. 2009; Foukal 2012), polar faculae (Le Mouël et al. 2019a, 2020a, 2019c) and white-light faculae areas (Foukal 1993, 2015). Another related sunspot proxy that might be useful is the sunspot decay rate. Hoyt&Schatten (1993) have noted that a fast decay rate suggests an enhanced solar convection, and hence a brighter Sun, while a slower rate signifies the opposite (Hoyt & Schatten 1993). Indeed, indications suggest that the decay rate during the Maunder Minimum was very slow (Hoyt & Schatten 1993), hence implying a dimmer Sun in the mid-to-late 17^th century. Owens et al. (2017) developed a reconstruction of the solar wind back to 1617 that suggests the solar wind speed was lower by a factor of two during the Maunder Minimum (Owens et al. 2017). Researchers have also considered records of various aspects of geomagnetic activity, since the Earth's magnetic field appears to be strongly influenced by solar activity (Le Moüel et al. 2019c; Duhau & Martínez 1995; Cliver et al. 1998; Richardson et al. 2002; Le Mouël et al. 2019b; Duhau & de Jager 2012).

There are many other solar proxies that might also be capturing different aspects of the long-term solar variability, e.g., see Livingston (1994), Soon et al. (2014) and Soon et al. (2015). In particular, it is worth highlighting the examination of cosmogenic isotope records, such as ¹⁴C or ¹⁰Be (Usoskin et al. 2009), since they are used by several of the TSI reconstructions we will consider. Cosmogenic isotope records have been regarded as longterm proxies of solar activity since the 1960s (Stuiver 1961; Suess 1965; Damon 1968; Suess 1968; Eddy 1977). Cosmogenic isotopes such as ¹⁴C or ¹⁰Be are produced in the atmosphere via galactic cosmic rays (GCRs). However, when solar activity increases, the solar wind reaching the Earth also increases. This tends to reduce the flux of incoming cosmic rays, thus reducing the rate of production of these isotopes and their quantity. These isotopes can then get incorporated into various long-term records, such as tree rings through photosynthesis. Therefore, by studying the changes in the relative concentrations of these isotopes over time in, e.g., tree rings, it is possible to construct an estimate of multidecadal-to-centennial and even millennial changes in average solar activity. Because the atmosphere is fairly well-mixed, the concentration of these isotopes only slowly changes over several years, and so the 8–14 year sunspot cycle can be partially reduced in these solar proxies. However, Stefani et al. (2020) still found a very good match between ¹⁴C or ¹⁰Be solar proxies and Schove (1955) estimates of solar activity maxima, which were based on historical aurora borealis observations back to 240 AD (Stefani et al. 2020). Moreover, the records can cover much longer periods, and so are particularly intriguing for studying multidecadal, centennial and millennial variability.

We note that several studies have tended to emphasize the similarities between various solar proxies (Lockwood & Fröhlich 2007; Gray et al. 2010; Lockwood 2012; Lean 2017; Gueymard 2018; Foukal 1998a). We agree that this is important, but we argue that it is also important to contrast as well as compare. To provide some idea of the effects that solar proxy choice, as well as TSI satellite composite choice, can have on the resulting TSI reconstruction, we plot in Figure 1 several different plausible TSI reconstructions taken from the literature and/or adapted from the literature. All nine reconstructions are provided in the Supplementary Materials.

Zoom In Zoom Out Reset image size

Foukal (2012) and Foukal (2015) followed a similar approach to the Wang et al. (2005) reconstruction but relied on slightly different solar proxies for the 20^th century pre-satellite era. Foukal (2012) used 10.7 cm solar microwave emissions for the 1947–1979 period and solar plage areas (from Ca II K spectroheliograms) for the 1916–1946 period. Foukal (2015) considered the faculae areas (from white light images) for the 1916–1976 period. In contrast, Wang et al. (2005) predominantly relied on the group sunspot number series for the pre-satellite era (after scaling to account for the sunspot/faculae/TSI relationships during the satellite era). These three different reconstructions are plotted as Figure 1(a), (c) and (i) respectively. All three reconstructions have a lot in common, e.g., they all have a very pronounced ∼11 year Solar Cycle component, and they all imply a general increase in TSI from the 19^th century to the mid-20^th century, followed by a general decline to present. However, there are two key differences between them. First, the Wang et al. (2005) reconstruction implies a slightly larger increase in TSI from the 19^th century to the 20^th century. On the other hand, while Foukal (2012) and Wang et al. (2005) imply the maximum TSI occurred in 1958, the Foukal (2015) reconstruction implies a relatively low TSI in 1958 and suggested two 20th century peaks in TSI – one in the late 1930s and another in 1979, i.e., the start of the satellite era. All three reconstructions imply that none of the global warming since at least 1979 could be due to increasing TSI, and in the case of the Foukal (2012) and Wang et al. (2005) reconstructions, since at least 1958.

Meanwhile, all three of those reconstructions were based on the PMOD satellite composite rather than the ACRIM composite. Therefore, in Figure 1(b) and (d), we have modified the Foukal (2012) and Foukal (2015) reconstructions using the ACRIM series for the 1980–2012 period instead of PMOD. We did this by rescaling the ACRIM time series to have the same mean TSI over the common period of overlap, i.e., 1980–2009.

Because the ACRIM composite implies a general increase in TSI from 1980 to 2000 followed by a general decrease to present, while the PMOD composite implies a general decrease in TSI over the entire period, this significantly alters the long-term trends. The modified version of Foukal (2012) implies that the 1958 peak in TSI was followed by an equivalent second peak in 2000. This suggests that at least some of the global warming from the 1970s to 2000 could have been due to increasing TSI, i.e., contradicting a key implication of Foukal (2012). The modified Foukal (2015) reconstruction is even more distinct. It implies that TSI reached an initial peak in the late 1930s, before declining until 1958, and then increasing to a maximum in 2000. As we will discuss in Section 3, this is broadly similar to many of the Northern Hemisphere temperature estimates. Therefore, the modified Foukal (2015) is at least consistent with the possibility of TSI as a primary driver of global temperatures over the entire 20^th century.

That said, while these plausible modifications can alter the relative magnitudes and timings of the various peaks and troughs in TSI, all these reconstructions would still be what Soon et al. (2015) and Scafetta et al. (2019) refer to as "low variability" reconstructions. That is, the multidecadal trends in TSI appear to be relatively modest compared to the rising and falling over the ∼ 11-year Solar Cycle component. As we will discuss in Section 2.6, many researchers have identified evidence for a significant ∼ 11-year temperature variability in the climate of the mid-troposphere to stratosphere (Labitzke & van Loon 1988; van Loon & Labitzke 2000; Labitzke 2005; Salby & Callaghan 2000, 2004, 2006; van Loon & Shea 1999; Labitzke & Kunze 2012; Camp & Tung 2007b; Frame & Gray 2010; Zhou & Tung 2013; van Loon & Shea 2000; Hood & Soukharev 2012), which has been linked to the more pronounced ∼ 11-year variability in incoming solar ultraviolet (UV) irradiance (Haigh 2003, 1994; Lean et al. 1997; Haigh & Blackburn 2006; Rind et al. 2008; Shindell et al. 2020; Kodera & Kuroda 2002; Hood 2003, 2016; Matthes et al. 2006). However, in terms of surface temperatures, the ∼ 11-year component seems to be only of the order of 0.02-0.2° C over the course of a cycle (Shaviv 2008; White et al. 1997, 1998; Douglass & Clader 2002; Camp & Tung 2007a; Zhou & Tung 2010). Some researchers have argued that the relatively large heat capacity of the oceans could act as a "calorimeter" to integrate the incoming TSI over decadal timescales, implying that the multidecadal trends are more relevant for climate change than annual variability (Shaviv 2008; Ziskin & Shaviv 2012; Reid 2000, 1987, 1991; Soon & Legates 2013), and others have argued that these relatively small temperature variations could influence the climate indirectly through e.g., altering atmospheric circulation patterns (van Loon & Meehl 2011; van Loon et al. 2012; Roy 2014, 2018; Christoforou & Hameed 1997; Meehl et al. 2009). However, this observation appears to have convinced many researchers (including the IPCC reports (IPCC 2013a)) relying on "low variability" reconstructions that TSI cannot explain more than a few tenths of a ° C of the observed surface warming since the 19th century, e.g., (IPCC 2013a; Crowley 2000; Laut 2003; Haigh 2003; Damon & Laut 2004; Foukal et al. 2006; Bard & Frank 2006; Benestad & Schmidt 2009; Gray et al. 2010; Lockwood 2012; Gil-Alana et al. 2014; Lean 2017). We will discuss these competing hypotheses and ongoing debates (which several co-authors of this paper are actively involved in) in Section 2.6, as these become particularly important if the true TSI reconstruction is indeed "low variability", i.e., dominated by the ∼ 11-year cycle.

On the other hand, let us consider the possibility that the true TSI reconstruction should be "high variability". In Figure 1(e)-(h), we consider four such "high variability" combinations, and we will discuss more in Section 2.4. All four of these reconstructions include a ∼11-year Solar Cycle component like the "low variability" reconstructions, but they imply that this quasi-cyclical component is accompanied by substantial multidecadal trends. Typically, the ∼11-year cycle mostly arises from the solar proxy components derived from the sunspot number datasets (as in the low variability reconstructions), while the multidecadal trends mostly arise from other solar proxy components.

Solanki & Fligge (1998) considered two alternative proxies for their multidecadal component and treated the envelope described by the two individual components as a single reconstruction with error bars. Solanki & Fligge (1999) also suggested that this reconstruction could be extended back to 1610 using the Group Sunspot Number time series of Hoyt & Schatten (1998) as a solar proxy for the pre-1874 period. However, in Figure 1(e) and (f), we treated both components as separate reconstructions, which we digitized from Solanki & Fligge (1998)'s figure 3, and extended up to 2012 with the updated ACRIM satellite composite. Both reconstructions are quite similar and, unlike the low variability estimates, imply a substantial increase in TSI from the end of the 19^th century to the end of the 20^th century. They also both imply that this long-term increase was interrupted by a decline in TSI from a mid-20^th century peak to the mid- 1960s. However, the reconstruction using Ca II K plage areas (Fig. 1e) implies that the mid-20^th century peak occurred in 1957, while the reconstruction incorporating Solar Cycle Lengths (Fig. 1f) implies the mid-20^th century peak occurred in the late 1930s and that TSI was declining in the 1940s up to 1965. In terms of the timing of the mid- 20^th century peak, it is worth noting that Scafetta (2012a) found a minimum in mid-latitude aurora frequencies in the mid-1940s, which is indicative of increased solar activity (Scafetta 2012a).

Figure 1(g) plots the updated Hoyt & Schatten (1993) TSI reconstruction. Although the original Hoyt & Schatten (1993) reconstruction was calibrated to the satellite era by relying on the NIMBUS7/ERB time series as compiled by Hoyt et al. (1992), it has since been updated by Scafetta & Willson (2014) and more recently by Scafetta et al. (2019) using the ACRIM composite until 2013 and the VIRGO and SORCE/TIM records up to present. The Hoyt & Schatten (1993) reconstruction is quite similar to the two Solanki & Fligge (1998) reconstructions, except that it implies a greater decrease in TSI from the mid-20^th century to the 1960s, and that the mid-20^th century peak occurred in 1947.

We note that there appear to be some misunderstandings in the literature over the Hoyt & Schatten (1993) reconstruction, e.g., Fröhlich & Lean (2002) mistakenly reported that "...Hoyt & Schatten (1993) is based on solar cycle length whereas the others are using the cycle amplitude" (Fröhlich & Lean 2002). Therefore, we should stress that, like Lean et al. (1995), the Hoyt & Schatten (1993) reconstruction did include both the sunspot numbers and the envelope of sunspot numbers, but unlike most of the other reconstructions, they also included multiple additional solar proxies (Hoyt & Schatten 1993). We should also emphasize that the Hoyt & Schatten (1998) paper describing the widely-used "Group Sunspot Number" dataset is a completely separate analysis, although it was partially motivated by Hoyt & Schatten (1993).

The Lean et al. (1995) reconstruction of Figure 1(h) also implies a long-term increase in TSI since the 19^th century and a mid-20^th century initial peak – this time at 1957, i.e., similar to Figure 1(e). The Lean et al. (1995) reconstruction was based on the Foukal & Lean (1990) reconstruction (Foukal & Lean 1990), and itself evolved into Lean (2000), which evolved into Wang et al. (2005), which in turn has evolved into the Coddington et al. (2016) reconstruction, which as we will discuss in Section 2.4 is a major component of the recent Matthes et al. (2017) reconstruction. However, Soon et al. (2015) noted empirically (in their fig. 9) that the main net effect of the evolution from Lean et al. (1995) to Lean (2000) to Wang et al. (2005) has been to reduce the magnitude of the multidecadal trends, i.e., to transition towards a "low variability" reconstruction. We note that Coddington et al. (2016) and Matthes et al. (2017) has continued this trend. Another change in this family of reconstructions is that the more recent ones have used the PMOD satellite composite instead of ACRIM (which is perhaps not surprising given that Lean was one of the PMOD team, as mentioned in Section 2.2, as well as being a co-author of all of that family of reconstructions).

Therefore, a lot of the debate over whether the high or low variability reconstructions are more accurate relates to the question of whether or not there are multidecadal trends that are not captured by the ∼11-year Solar Cycle described by the sunspot numbers. This overlaps somewhat with the ACRIM vs. PMOD debate since the PMOD composite implies that TSI is very highly correlated to the sunspot number records (via the correlation between sunspots and faculae over the satellite era), whereas the ACRIM composite is consistent with the possibility of additional multidecadal trends between solar cycles (Scafetta 2011; Scafetta & Willson 2014; Soon et al. 2015; Scafetta et al. 2019; Dudok de Wit et al. 2017; Fröhlich 2012; Gueymard 2018).

This has been a surprisingly challenging problem to resolve. As explained earlier, the ∼11-year cyclical variations in TSI over the satellite era are clearly well correlated to the trends in the areas of faculae, plages, as well as sunspots over similar timescales (Willson & Hudson 1988; Lean & Foukal 1988; Foukal & Lean 1990; Willson & Hudson 1991; Hoyt et al. 1992; Wade 1995; Willson 1997). However, on shorter timescales, TSI is actually anti-correlated to sunspot area (Willson et al. 1981; Eddy et al. 1982). Therefore, the ∼11-year rise and fall in TSI in tandem with sunspot numbers cannot be due to the sunspot numbers themselves, but appears to be a consequence of the rise and fall of sunspot numbers being commensally correlated to those of faculae and plages. However, Kuhn et al. (1988) argued that "...solar cycle variations in the spots and faculae alone cannot account for the total [TSI] variability" and that "...a third component is needed to account for the total variability" (Kuhn et al. 1988). Therefore, while some researchers have assumed, like Lean et al. (1998), that there is "...[no] need for an additional component other than spots or faculae" (Lean et al. 1998), Kuhn et al. (Kuhn & Libbrecht 1991; Kuhn et al. 1998, 1999; Kuhn 2004; Armstrong 2004) continued instead to argue that "sunspots and active region faculae do not [on their own] explain the observed irradiance variations over the solar cycle(Kuhn 2004) and that there is probably a "third component of the irradiance variation" that is a "nonfacular and nonsunspot contribution" (Kuhn & Libbrecht 1991). Work by Li, Xu et al. is consistent with Kuhn et al.'s assessment, e.g., Refs. (Li et al. 2012, 2016; Xu et al. 2017; Li et al. 2020a), in that they have shown: that TSI variability can be decomposed into multiple frequency components (Li et al. 2012); that the relationships are different between different solar activity indices and TSI (Li et al. 2016, 2020a); and that the relationship between sunspot numbers and TSI varies between cycles (Xu et al. 2017). Indeed, in order to accurately reproduce the observed TSI variability over the two most recent solar cycles using solar disk images from ground-based astronomical observatories, Fontenla & Landi (2018) needed to consider nine different solar features rather than the simple sunspot and faculae model described earlier.

In summary, there are several key debates ongoing before we can establish which TSI reconstructions are most accurate:

• 1  
  Which satellite composite is most accurate? In particular, is PMOD correct in implying that TSI has generally decreased over the satellite era, or is ACRIM correct in implying that TSI increased during the 1980s and 1990s before decreasing?
• 2  
  Is it more realistic to use a high variability or low variability reconstruction? Or, alternatively, has the TSI variability been dominated by the ∼11-year Solar Cycle, or have there also been significant multidecadal trends between cycles?
• 3  
  When did the mid-20^th century peak occur, and how much and for how long did TSI decline after that peak?

The answers to these questions can substantially alter our understanding of how TSI has varied over time. For instance, Velasco Herrera et al. (2015) used machine learning and four different TSI reconstructions as training sets to extrapolate forward to 2100 AD and backwards to 1000 AD (Velasco Herrera et al. 2015). The results they obtained had much in common, but also depended on whether they used PMOD or ACRIM as well as whether they used a high or low variability reconstruction. As an aside, the forecasts from each of these combinations implied a new solar minimum starting in 2002-2004 and ending in 2063-2075. If these forecasts are correct, then in addition to the potential influence on future climate change, such a deficit in solar energy during the 21^st century could have serious implications for food production; health; in the use of solar-dependent resources; and more broadly could affect many human activities (Velasco Herrera et al. 2015).

Soon et al. (2015) identified eight different TSI reconstructions (see fig. 8 in that paper) (Soon et al. 2015). Only four of these reconstructions were considered by the CMIP5 modeling groups for the hindcasts that were submitted to the IPCC AR5: Wang et al. (2005) described above, as well as Krivova et al. (2007); Steinhilber et al. (2009); and Vieira et al. (2011). Coincidentally, all four implied very little solar variability (and also a general decrease in TSI since the 1950s). However, Soon et al. (2015) also identified another four TSI reconstructions that were at least as plausible – including the Hoyt & Schatten (1993) and Lean et al. (1995) reconstructions described above. Remarkably, all four implied much greater solar variability. These two sets are the "high solar variability" and "low solar variability" reconstructions discussed in Section 2.3 which both Soon et al. (2015) and more recently Scafetta et al. (2019) have referred to.

Since then, eight additional estimates have been proposed – four low variability and four high variability. Coddington et al. (2016) have developed a new version of the Wang et al. (2005) estimate that has reduced the solar variability even further (it relies on a sunspot/faculae model based on PMOD). Recently, Matthes et al. (2017) took the mean of the Coddington et al. (2016) estimate and the (similarly low variability) Krivova et al. (2007, 2010) estimate, and proposed this as a new estimate. Moreover, Matthes et al. (2017) recommended that their new estimate should be the only solar activity dataset considered by the CMIP6 modeling groups (Matthes et al. 2017). Clearly, Matthes et al.'s (2017) recommendation to the CMIP6 groups goes against the competing recommendation of Soon et al. (2015) to consider a more comprehensive range of TSI reconstructions.

In Figure 2, we plot the four "low solar variability" reconstructions from Soon et al. (2015) as well as these two new "low variability" estimates along with another two estimates by Dr. Leif Svalgaard (Stanford University, USA), which have not yet been described in the peer-reviewed literature but are available from Svalgaard's website [https://leif.org/research/download-data.htm, last accessed on 2020/03/27], and have been the subject of some discussion on internet forums.

Zoom In Zoom Out Reset image size

Recently, Egorova et al. (2018) proposed four new "high variability" estimates that built on the earlier Shapiro et al. (2011) estimate. The Shapiro et al. (2011) estimate generated some critical discussion (Schmidt et al. 2012; Judge et al. 2012; Shapiro et al. 2013) (see Sect. 2.5.4). Egorova et al. (2018) have taken this discussion into account and proposed four new estimates utilizing a modified version of the Shapiro et al. (2011) methodology. Therefore, in Figure 3, we plot the four "high solar variability" reconstructions from Soon et al. (2015) as well as these four new "high variability" estimates. This provides us with a total of 16 different TSI reconstructions. Further details are provided in Table 1 and in the Supplementary Materials. For interested readers, we have also provided the four additional TSI reconstructions discussed in Figure 1 in the Supplementary Materials.

Zoom In Zoom Out Reset image size

The primary focus of the new analysis in this paper (Sect. 5) is on evaluating the simple hypothesis that there is a direct linear relationship between incoming TSI and Northern Hemisphere surface air temperatures. As will be seen, even for this simple hypothesis, a remarkably wide range of answers is still plausible. However, before we discuss in Section 3 what we currently know about Northern Hemisphere surface air temperature trends since the 19^th century (and earlier), it may be helpful to briefly review some of the other frameworks within which researchers have been debating potential Sun/climate relationships.

The gamut of scientific literature which encompasses the debates summarized in the following subsections (2.5 and 2.6) can be quite intimidating, especially since many of the articles cited often come to diametrically opposed conclusions that are often stated with striking certainty. With that in mind, in these two subsections, we have merely tried to summarize the main competing hypotheses in the literature, so that readers interested in one particular aspect can use this as a starting point for further research. Also, several of the co-authors on this paper have been active participants in some of the debates we will be reviewing. Hence, there is a risk that our personal assessments of these debates might be subjective. Therefore, we have especially endeavored to avoid forming definitive conclusions, although many of us have strong opinions on several of the debates we will discuss here.

The various debates that we consider in this subsection (2.5) can be broadly summarized as being over whether variations in solar activity have been a major climatic driver in the past. We stress that a positive answer does not in itself tell us how much of a role solar activity has played in recent climate change. For instance, several researchers have argued that solar activity was a major climatic driver until relatively recently, but that anthropogenic factors (chiefly anthropogenic CO₂ emissions) have come to dominate in recent decades (Crowley 2000; Lockwood & Fröhlich 2007; Hegerl et al. 2007; Lean & Rind 2008; Benestad & Schmidt 2009; Gray et al. 2010; Lean 2017; Beer et al. 2000; de Jager et al. 2010; Lean et al. 1995). However, others counter that if solar activity was a major climatic driver in the past, then it is plausible that it has also been a major climatic driver in recent climate change. Moreover, if the role of solar activity in past climate change has been substantially underestimated, then it follows that its role in recent climate change may also have been underestimated (Bond et al. 2001; Scafetta & West 2006b; Svensmark 2007; Courtillot et al. 2007, 2008; Singer & Avery 2008; Lüning & Vahrenholt 2015, 2016; Mörner et al. 2020; Friis-Christensen & Lassen 1991; Lassen & Friis-Christensen 1995; Soon et al. 1994; Scafetta 2013, 2020; Scafetta et al. 2016; Loehle&Singer 2010; Shaviv & Veizer 2003; Judge et al. 2020; Baliunas & Jastrow 1990; Zhang et al. 1994; Zhao et al. 2020).

2.5.1. Evidence for long-term variability in both solar activity and climate

2.5.2. Similarity in frequencies of solar activity metrics and climate changes

2.5.3. Sun-Planetary Interactions as a plausible mechanism for long-term solar variability

2.5.4. Analogies of solar variability with the variability of other "Sun-like" stars

If you consider all of the TSI reconstructions among the "low variability estimates" (Fig. 2), except for the Steinhilber et al. (2009) reconstruction which is based on cosmogenic isotope proxies, it could appear that the most significant feature is the short-term maximum-minimum sunspot cycle fluctuations which occur with a roughly-11 year period (i.e., the "Schwabe cycle"). Therefore, initially, it might be supposed that the influence of solar variability on the Earth's climate should be most obvious over the course of each sunspot cycle. This applies even more so if you treat the raw Sunspot Number (SSN) record as a proxy for TSI, since the SSN falls to zero during every cycle (Gil-Alana et al. 2014).

This has been a puzzle for the community since the beginning of modern research into possible Sun/climate connections, since the fluctuations in global surface air temperature (for instance) over a sunspot cycle are relatively small at best (Gray et al. 2010; Scafetta 2009), and often quite ambiguous (Pittock 1983).

Typically, the peak-to-trough variability in global surface air temperatures over a sunspot cycle is estimated empirically at about 0.1° C (Gray et al. 2010; Scafetta 2009), although Scafetta (2009) notes that the estimates of this "11-year solar cycle signature" in the literature vary from about 0.05° C to about 0.2° C (Scafetta 2009). He also notes that typical climate models are unable to simulate even this modest temperature variability over a solar cycle, with some climate models predicting the solar cycle signature to be as low as 0.02 – 0.04° C (Scafetta 2009). Partly on this basis, he suggests that there are "...reasons to believe that traditional climate models cannot faithfully reconstruct the solar signature on climate and are significantly underestimating it" (Scafetta 2009).

In any case, if you assume that (a) the low variability TSI estimates are more reliable than the high-variability estimates, and (b) there is a linear relationship between TSI and global (or hemispheric) surface air temperatures, these relatively low 11-year Solar Cycle signature estimates would initially appear to put a very modest upper bound on the maximum contribution of solar variability to the Northern Hemisphere surface temperature trends since the 19^th century. In Section 5, we will compare and contrast the linear fits using the high and low-variability TSI estimates, i.e., we will be implicitly evaluating the first assumption. However, there is also a considerable body of literature critically evaluating the second assumption from several different avenues. Therefore, in this section, we will briefly review some of the main attempts to resolve this apparent "11-year paradox".

The apparent paradox could indicate that the Sun affects the climate by other covarying aspects of solar variability (changes in the UV component, GCR fluxes, etc.) rather than just TSI changes. Indeed, much of the literature over the last few decades has suggested that we should not be only looking for a direct linear relationship between TSI and global surface air temperatures, but rather considering the possibility of more indirect and/or subtle Sun/climate relationships. Some of the main hypotheses are summarized schematically in Figure 4.

Zoom In Zoom Out Reset image size

2.6.1. "Top-down" vs. "bottom-up" mechanisms

2.6.2. "The ocean as a buffer": Ocean heat capacity as a "filter capacitor"–based buffering mechanism

2.6.3. Sun-climate effects are more pronounced in certain regions

2.6.4. Galactic Cosmic Ray-driven amplification mechanisms

2.6.5. Short-term orbital effects

2.6.5.1.  The difference between the average Earth-Sun distance (1 AU) and the daily Earth-Sun distance

It is also worth briefly distinguishing between changes in the solar irradiance leaving the Sun and changes in the solar irradiance reaching the Earth. Much of the interest in understanding the changes in solar activity have focused on the former. As a result, the TSI is typically described in terms of the output reaching 1 AU, i.e., the average distance of the Earth from the Sun. This applies to most of the TSI reconstructions discussed in this paper. However, as discussed in Soon et al. (2015), because the Earth's orbit of the Sun is elliptical rather than circular, the physical distance of the Earth from the Sun varies quite a bit over the course of the year. Currently, the Earth receives 6.5% more TSI (88 W m⁻²) in January (i.e., the Northern Hemisphere winter) during the perihelion (the Earth's closest point to the Sun) than in July (i.e., the Northern Hemisphere summer) during the aphelion (the Earth's furthest point from the Sun). Therefore, if we are interested in the effects of changes in the solar activity on the Earth's climate, then arguably we are more interested in the changes in TSI actually reaching the Earth, rather than the changes in TSI reaching a distance of 1 AU.

Because the seasonal cycles of the Earth's orbit are almost identical each year, it might initially be supposed that the annual averages of the TSI reaching the Earth and that reaching 1 AU should be perfectly correlated with each other. If so, then this would mean that, when averaged over the year, this difference between the TSI reaching the Earth compared to that at 1 AU would be trivial. In that case, using 1 AU estimates for evaluating the potential effects of varying TSI on the Earth's climate would be easier, since the variations in TSI from year-to-year are easier to see in the TSI data at 1 AU compared to that at Earth's distance (as will be seen below).

This appears to be an implicit assumption within much of the literature evaluating the effects of varying TSI on the Earth's climate. However, we note that there are subtle, but often substantial, differences that can arise between the annual averages at 1 AU versus Earth's distance since most of the trends in TSI (including the ∼11-year cycle) are on timescales that are not exact integer multiples of the calendar year.

We appreciate that mentally visualizing the differences that can arise between the two different estimates is quite tricky, even if you have a high degree of visual spatial intelligence. Therefore, in order to demonstrate that there are indeed subtle but significant differences between the annual averages of TSI at 1 AU versus the TSI reaching the Earth, we compare and contrast the SORCE/TIM TSI datasets for 1 AU and Earth's distance in Figure 5. The SORCE/TIM datasets are the results from a single satellite that operated from 2003 to 2020, and it is particularly relevant because the data are reported for both 1 AU and at the in-situ Earth's orbit. This period covers roughly 1.5 Solar Cycles – Figure 5(b).

Zoom In Zoom Out Reset image size

Before describing the results in Figure 5, we should note some technical points on this analysis. We downloaded the daily-resolved datasets for both versions from http://lasp.colorado.edu/data/sorce/tsi_data/ (accessed 26/06/2020). Unfortunately, the full SORCE TSI data records from 2003 February 25 through 2020 February 25 do not have continuous daily measured values for either the in-situ Earth's orbit or the 1 AU distance adjusted data. On average there are 11 missing days for all years except 2013, 2003 and 2020. 2013 was missing 160 daily values, and discounting 2013, the average for the full 2004–2019 period is 22.5 missing days per year. To fill in the missing daily TSI values, we applied a more sophisticated method than just a linear interpolation. For the 1 AU TSI values, we rely on an artificial intelligence algorithm that matched not only the amplitude of the measured daily TSI values but also the spectral properties of the measured SORCE TSI record. In addition, PMOD TSI data were utilized and were calibrated to SORCE's TIM daily series applying the method proposed by Soon et al. (2019a). The PMOD data between 2003 and 2017 were standardized with characteristics of the TIM and then added to the TIM records. The data between 2018 and 2020 that were missing in the TIM series were consecutively less than 5 days and therefore were estimated using the method of Radial Basis Function Artificial Neural Networks (RBFANN). Our RBFANN has three layers of neurons: one input (objective data in this case from TIM), one hidden and one output (matching the high- or low- frequency spectral properties). For the daily TSI values at Earth's orbit, another set of RBFANN was constructed to fill in the missing daily values. Therefore, by applying these advanced and more elegant techniques than interpolation per se, we can generate the complete daily SORCE/TIM's TSI composite time series for both Earth orbit and at 1 AU perspectives.

The daily results are shown in Figures 5(b) and 5(c), with the interpolated points indicated with dashed lines (and slightly different colors). The annual means of both versions are compared in Figure 5(d). Note that the y-axes each have a different range. This is because, as noted above, the seasonal cycle in the TSI reaching the Earth is of the order of 90 W m⁻², while the variability in TSI over the solar cycle is only of the order of a few W m⁻². As a result, the pronounced ∼11 year Solar Cycles that can be seen in the 1 AU plot of Figure 5(b) are barely noticeable when viewed on the scale of Figure 5(c).

That said, when the annual averages of both time series are calculated, this seasonal cycle is no longer an issue, and the two time series can be directly compared, as in Figure 5(d). However, as can be seen from Figure 5(e), while the two annual time series are broadly similar, albeit with the 1 AU values being slightly lower than that at Earth distance, the differences between the two time series vary slightly from year-to-year. Over the 16 year period, the differences between the two annual averages varied from +0.35 W m⁻² (2013) to –0.06 W m⁻² (2014), i.e., a range of 0.41 W m⁻². It could be argued that the 2014 estimate is anomalous in that this was the year with the most data interpolation. But, even neglecting that year, the differences between the two averages varied between 0.34 and 0.15 W m⁻², i.e., a range of nearly 0.2 W m⁻². For comparison, the difference between the maximum and minimum annual TSI at 1 AU over the same period was 0.9 W m⁻². So, while small as a percentage of the absolute TSI, these subtle differences are not insignificant.

These differences between the annual averages for each year might initially be surprising. For the analysis in Figure 5, we are assuming that the SORCE datasets at both 1 AU and the Earth's distance are reliable. We also are assuming that the interpolations we have carried out (described earlier) are reasonable. However, even if either of those assumptions are problematic, we should stress that the fact that there are differences between the annual means for both versions that vary from year to year is actually to be expected on statistical grounds. The general principle can be understood once we recognize (a) the elliptical shape of the Earth's orbit and (b) that the ∼11-year Solar Cycle does not fall exactly on the calendar year. This means that the times of the year during which a given rise or fall in TSI occurs can make a difference. For instance, if the maximum of a solar cycle occurred during January, then this will lead to a greater mean TSI for the year than if it had occurred during July of the same year. This is because the Earth is currently closer to the Sun in January than in July.

To clarify, if the trends in TSI over a given calendar year are reasonably linear statistically, i.e., it has a constant slope, then this seasonality should not make much difference to the annual mean TSI. This is regardless of the slope itself, i.e., whether the trend is rising, falling or near-zero. However, if the trends for that year are non-constant, then the annual mean may be slightly higher or lower depending on whether the Earth is closer to perihelion or aphelion when the changes in the trends in TSI occur.

More generally, the annual average TSI reaching the Earth depends not just on the changes in TSI, but the times of the year over which those changes occur. The SORCE data in Figure 5 only covers roughly 1.5 Solar Cycles, but in principle, the same could apply to any other multidecadal trends which might be occurring in addition to the ∼11-year cycle.

Since all of the TSI reconstructions discussed in Sections 2.2–2.5 are calculated in terms of the annual averages at 1 AU, for this paper we will limit our analysis to this. However, we encourage researchers who have until now limited their analysis of TSI variability to that at 1 AU to consider this extra complication in future research.

As can be seen from Figure 5(e), the changing differences between the annual average TSI at 1 AU versus that reaching the Earth are subtle, but non-trivial. In addition to this complication in terms of the annual averages, a comprehensive analysis of the effects of TSI on the Earth's climate should consider the seasonal changes in the different latitudinal distributions of the incoming TSI, due to the seasonal orbital variation of the Earth. There are several different aspects to this, but for simplicity they are collectively referred to as "orbital forcings".

2.6.5.2.  Comparison with long-term orbital forcing

The theory that changes in atmospheric CO₂ are a primary driver of climate change was originally developed by Arrhenius in the late 19^th century as a proposed explanation for the transitions between glacial and interglacial periods during the ice ages (Arrhenius 1896). The existence of these dramatic climatic changes on multi-millennial timescales was only established in the 19^th century and was one of the great scientific puzzles of the time. [As an aside we note that glaciologically speaking, an "ice age" is usually defined as a period where large permanent ice sheets are present in both hemispheres. These ice sheets can substantially expand during "glacial periods" and retreat during "interglacial periods". As Greenland and Antarctica both currently have large ice sheets, we are currently in an interglacial period (the "Holocene") within an ice age, even though colloquially, the term "ice age" is popularly used just to describe the "glacial periods".

This CO₂ driven explanation for the glacial/interglacial transitions was later criticized by, e.g., Ångström (1901) and Simpson (1929). However, it was later revived by Callendar (1938) who extended the theory to suggest that anthropogenic CO₂ emissions were also the primary driver of the warming from the late-19^th century to mid-1930s (Callendar 1938), and Plass (1956) who speculatively proposed (anticipating that this would encourage scientific debate) that atmospheric CO₂ was the primary driver of climate change on most timescales (Plass 1956).

A competing hypothesis that several 19^th-century researchers proposed, e.g., Adhmar, and later Croll (Fedorov 2019a; Imbrie 1982; Bol'shakov et al. 2012; Sugden 2014), was that long-term cyclical changes in the Earth's orbit around the Sun were the driver of the glacial/interglacial transitions. In the early 20^th century, Milankovitch carried out an extensive series of calculations that demonstrated that there are several important cyclical variations in the Earth's orbit that vary over tens of thousands of years and that these influence the incoming solar radiation at different latitudes for each of the seasons (Cvijanovic et al. 2020; Szarka et al. 2021).

In the 1970s, relatively high precision estimates of the timings of the glacial/interglacial transitions from ocean sediment cores appeared to vindicate the Milankovitch ice age theory (Imbrie 1982; Hays et al. 1976). That is, when frequency analyses like those in Sections 2.5.2 and 2.5.3 were carried out on the deep-sea sediment cores, they suggested that the past climate changes of the last few hundred thousand years were dominated by periodicities of ∼90 000 – 120 000 years and to a lesser extent 40 000–42 000 years, and also had peaks at 22 000–24 000 years and 18 000–20 000 years. These peaks were approximately similar to the main astronomical cycles calculated by Milankovitch: 41 000 years; 23 000 years; 19 000 years and to a lesser extent ∼100 000 years. As a result, analogous to the arguments described in Sections 2.5.2–2.5.3, it was argued that the glacial/interglacial transitions were indeed driven by orbital forcings (Imbrie 1982; Hays et al. 1976). Later analysis in terms of ice core measurements also appeared to confirm this theory (Lorius et al. 1985, 1992; Petit et al. 1999).

This appears to have convinced most of the scientific community that the Milankovitch orbital-driven explanation for the glacial/interglacial transitions is correct, and this currently seems to be the dominant paradigm within the literature (Cvijanovic et al. 2020; Szarka et al. 2021; Petit et al. 1999; Kawamura et al. 2007; Maslin & Ridgwell 2005; Roe 2006; Imbrie et al. 1992; Lisiecki & Raymo 2005) including the IPCC reports (IPCC 2013a). Ironically, this means that the original CO₂-driven theory for the glacial/interglacial transitions as proposed by Arrhenius (1896); Callendar (1938); Plass (1956) has been largely discarded even though the current theory that recent climate change has been largely driven by changes in CO₂ was developed from that early theory. That said, we emphasize that this is not necessarily a contradiction in that several researchers argue that changes in atmospheric CO₂ concentrations caused by orbitally-driven warming or cooling might act as a positive feedback mechanism (IPCC 2013a; Cvijanovic et al. 2020; Petit et al. 1999; Kawamura et al. 2007; Maslin & Ridgwell 2005; Roe 2006; Imbrie et al. 1992; Lisiecki & Raymo 2005).

We note that there are possible problems with the theory that the Milankovitch orbitals are the primary driver of the glacial/interglacial transitions. To clarify, the Milankovitch orbital variations are clearly climatically significant as discussed above. Also, the idea that some combination of these variations could provide the explanation for the glacial/interglacial transitions seems plausible and intuitive. Indeed, the approximate similarity in the timings of both phenomena is intriguing. However, as will be discussed below, there is still considerable debate over the exact causal mechanisms and over which specific aspects of the Milankovitch orbital variations would drive such dramatic long-term climate changes and why. That said, each of these problems has been countered, and the current consensus among the scientific community is indeed that the glacial/interglacial transitions are driven by orbital forcings.

With that in mind, we will briefly touch on some of the concerns that have been raised over the Milankovitch theory of ice ages and also the corresponding defences that have been made. As for the SPI theories, many of the main concerns with the theory have been highlighted by proponents of the theory. For example, Hays et al. (1976) noted that the dominance of the ∼100 000 year periodicities in the climate records was unexpected since the relevant Milankovitch ∼100 000 year orbital cycles (associated with changes in eccentricity) had not been expected to have much significance for climate change, and Milankovitch and others had assumed the ∼41 000 year obliquity cycle would have been more climatically significant in terms of glacial/interglacials (Hays et al. 1976).

Imbrie (1982) noted another problem: because there are so many different aspects of the timings of the various orbitals in terms of different latitudes and seasons, there is a danger that researchers could cherry-pick the metric that best fits the data for their hypothesis:

"There has also been a tendency for investigators to believe they could model the response of the system from a radiation curve representing the input at a single latitude and season [...] Since no one could be sure which insolation curve, if any, was the crucial one, investigators had great freedom to choose a curve that resembled a particular set of data. Understandably the resulting ambiguity did much to undermine confidence in the validity of the time domain prediction." – Imbrie (1982), p413 (Imbrie 1982)

An additional puzzle is that, while the dominant peak on orbital timescales for the last ∼1 million years appears to be ∼100 000 years, for most of the Pleistocene Ice Age until then, the dominant peak appears to have been the (expected) 41 000 years. This has been dubbed the "mid-Pleistocene transition", and there has been much debate over its explanation (Maslin & Ridgwell 2005; Raymo et al. 2006; Ashkenazy & Tziperman 2004; Tziperman et al. 2006; Rial et al. 2013; Huybers & Wunsch 2005).

Another potential concern is the so-called "Devils Hole" record. Winograd et al. (1992) obtained a continuous 500 000 year climate record from a core in Devils Hole, NV (USA), which matched quite well with both the Antarctic ice cores and various ocean sediment cores (Winograd et al. 1992). However, the Devils Hole record had a key difference from the other estimates – it implied that the previous interglacial period ("Termination II") began several millennia before Milankovitch theory predicted it should have (Winograd et al. 1992; Broecker 1992; Karner & Muller 2000). Therefore, if this estimate transpires to be accurate, then it creates a significant problem for the Milankovitch ice age theory. That said, proponents of the Milankovitch explanation have criticized the reliability of the earlier Termination II timing implied by the Devils Hole record (Shackleton 1993; Imbrie et al. 1993; Shakun et al. 2011; Moseley et al. 2016a,b). However, rebuttals have been offered (Winograd & Landwehr 1993; Ludwig et al. 1993; Winograd et al. 2006; Coplen 2016; Winograd 2016).

If the glacial/interglacial transitions are not primarily-driven by the Milankovitch orbital cycles, then it raises the question as to why the glacial/interglacial transitions occur. Several different explanations have been proposed (Sharma 2002; Ehrlich 2007; Roe & Allen 1999; Wunsch 2004; Muller & MacDonald 1997b,a; Kirkby et al. 2004; Puetz et al. 2016; Ellis & Palmer 2016; Marsh 2014). However we note that all of these cited alternative hypotheses are openly speculative, and that several of these studies (and others) argue that Milankovitch orbital cycles are at least partially involved (Ashkenazy & Tziperman 2004; Tziperman et al. 2006; Huybers & Wunsch 2005; Roe & Allen 1999; Wunsch 2004; Ellis & Palmer 2016).

At any rate, regardless of the role the Milankovitch orbital cycles play in the glacial/interglacial transitions, we emphasize that the changes in the latitudinal and seasonal variability of the incoming TSI over these cycles are clearly climatically important. Indeed, while much of the research until now has focused on the gradual variations over millennial timescales or longer, e.g., Huybers & Denton (2008); Davis & Brewer (2009); Berger et al. (2010), in recent years, some groups have begun emphasizing the significance of "short-term orbital forcing" (STOF) (Cionco & Soon 2017), i.e., secular drifts on multidecadal to centennial timescales in the average daily insolation at different latitudinal bands for different seasons (Cionco & Soon 2017; Fedorov 2019a; Szarka et al. 2021; Cionco et al. 2018, 2020; Fedorov 2019b; Fedorov & Kostin 2020).

A challenge that arises when considering the Earth's orbital variability on these shorter timescales is that it is no longer sufficient to consider the Earth/Sun relationship in isolation. The perturbations of the Earth's orbit around the Sun by the Moon and by the other planets also need to be considered. For this reason, there is some overlap between much of the ongoing research into STOF (Cionco & Soon 2017; Cionco et al. 2018, 2020; Fedorov 2019b,b; Fedorov & Kostin 2020) and the calculations within the SPI theories described in Section 2.5.3. However, we stress that the two fields of research are distinct. The former field is interested in how the perturbations of the other planetary bodies on the Earth's orbit around the Sun influence the distribution of the incoming TSI reaching the Earth. The latter field is concerned with what effects (if any) the orbits of the various planetary bodies might have on TSI (or more broadly, solar activity) itself. Scafetta et al. (2020) adds a third possible planetary mechanism that could be influencing the Earth's climate, by arguing that planetary configurations could directly modulate the flux of the interplanetary dust reaching the Earth's atmosphere thereby potentially influencing cloud cover (Scafetta et al. 2020).

In any case, for simplicity, most discussions on the climatic significance of orbital variations have tended to pick a particular latitudinal band and season, e.g., the average insolation during the Northern Hemisphere summer at 65°N. However, as discussed in the quote above from Imbrie (1982), identifying the particular latitude and season whose curve is "most important" is quite subjective. For instance, Huybers & Denton (2008) argue that "the duration of Southern Hemisphere summer is more likely to control Antarctic climate than the intensity of Northern Hemisphere summer with which it (often misleadingly) covaries" (Huybers & Denton 2008). Indeed, Davis & Brewer (2009) argued that a more climatically relevant factor is the "latitudinal insolation gradient (LIG)", which in turn "... creates the Earth's latitudinal temperature gradient (LTG) that drives the atmospheric and ocean circulation" (Davis & Brewer 2009). Note that this latter LTG concept is equivalent to the "EPTG" parameter considered by Lindzen (1994) and Soon & Legates (2013) and comparable to the "meridional heat transfer (MHT)" concept of Fedorov (2019b). Another related index is the zonal index pressure gradient, which Mazzarella & Scafetta (2018) argue is linked to the NAO, LOD and atmospheric temperatures, as well as directly influencing atmospheric circulation patterns (Mazzarella & Scafetta 2018).

Although Davis & Brewer (2009)'s focus was on timescales of millennia (Davis & Brewer 2009), Soon (2009) noted that the variability of the LIG can also be significant on multidecadal-to-centennial timescales (Soon 2009). Therefore, recent work into studying the variability of the LIG (or Fedorov (2019b)'s related "meridional insolation gradient (MIG)" (Fedorov 2019b)) on these shorter timescales (Cionco et al. 2018, 2020; Fedorov 2019b; Davis & Brewer 2011) could be a particularly important "STOF" for understanding recent climate changes.

Many researchers have relied on weather station records to estimate Northern Hemisphere (and global) land surface air temperature trends since at least the late-19^th/early-20^th century (Soon et al. 2015; Lugina et al. 2006; Hansen et al. 2010; Lenssen et al. 2019; Lawrimore et al. 2011; Muller et al. 2014; Jones et al. 2012; Cowtan & Way 2014; Xu et al. 2018). The results are typically similar and imply an average warming trend of ∼1° C century⁻¹ since the late-19^th century. Thus, most studies assumed these estimates are robust and well-replicated. Recently, however, Soon et al. (2015) constructed a new estimate using only rural (or mostly rural) stations taken from four regions with a high density of rural stations, and they found that this new estimate yielded noticeably different results.

Soon et al. (2015)'s estimate incorporated more than 70% of the available rural stations with records covering both the early and late 20^th century. "Rural" was defined in terms of both low night-brightness and low associated population using the urban ratings provided by the dataset. If they are correct that urbanization bias has not been adequately corrected for in the standard estimates, then the Soon et al. (2015) estimate should be more representative of the Northern Hemisphere land surface temperature trends than the standard estimates. McKitrick & Nierenberg (2010), Scafetta & Ouyang (2019) and more recently Scafetta (2021) and Zhang et al. (2021) offer some support for this.

On the other hand, several researchers have claimed that the standard estimates are not majorly affected by urbanization bias. Some argue that urbanization bias is only a small problem for global and hemispheric temperature trends, e.g., Jones et al. (1990), Parker (2006), Wickham et al. (2013). Others concede that urbanization bias is a concern for the raw station data, but argue that after statistical homogenization techniques (usually automated) have been applied to the data, most of the nonclimatic biases (including urbanization bias) are removed or substantially reduced, e.g., Peterson et al. (1999), Menne & Williams (2009), Hausfather et al. (2013), Li & Yang (2019), Li et al. (2020b). If either of these two groups is correct, then considering all available station records (urban and rural) should be more representative of the Northern Hemisphere land surface temperature trends.

That said, we note that the methodology and/or justification of each of the above studies (and others reaching similar conclusions) has been disputed, e.g., see Connolly & Connolly (2014a) for a critique of 13 such studies (Connolly & Connolly 2014a). We also recommend reading the review comments of McKitrick who was a reviewer of an earlier version of Wickham et al. (2013) and has made his reviews publicly available at: https://www.rossmckitrick.com/temperature-data-quality.html. Also, some of us have argued that these statistical homogenization techniques inadvertently lead to "urban blending" whereby some of the urbanization bias of the urban stations is aliased onto the trends of the rural stations. This means that, counterintuitively, the homogenized station records (whether rural or urban) may often end up less representative of climatic trends than the raw rural station records (Soon et al. 2018, 2019b).

As in Soon et al. (2015), our analysis in this section is based on version 3 of NOAA's Global Historical Climatology Network (GHCN) monthly temperature dataset (Lawrimore et al. 2011). This GHCN dataset was also used by the Hansen et al. (2010) and Lawrimore et al. (2011) estimates, and it became a major component of the Jones et al. (2012), Xu et al. (2018) and Muller et al. (2014) estimates, which will be considered in Section 3.2. However, all these estimates incorporated urban as well as rural stations.

Version 4 of the GHCN monthly temperature dataset Soon et al. (2015). This newer version provides a much larger number of stations (∼20 000 versus 7 200 in version 3), although most of the new stations have short station records (only a few decades) and ∼40% of the stations cover only the contiguous United States. However, unlike version 3, this new dataset does not yet include any estimates of how urbanized each of the stations are. Some of us have begun developing new estimates of the urbanization of the version 4 stations on a regional basis – see Soon et al. (2018) for the results for China. This research is ongoing and has not yet been completed. In the meantime, since it has been claimed that the changes in long-term trends introduced by the switch to version 4 are relatively modest (Menne et al. 2018), we believe it is useful to update and consider in more detail the Soon et al. (2015) "mostly rural" Northern Hemisphere temperature time series using Version 3.

We downloaded a recent update of the Version 3 dataset from https://www.ncdc.noaa.gov/ghcnm/ [version 3.3.0.20190821, accessed 21/08/2019]. This dataset has two variants – an unadjusted dataset that has only quality control corrections applied and an adjusted dataset that has been homogenized using the automated Menne & Williams (2009) algorithm. A major component of the GHCN dataset is a subset for the contiguous United States called the US Historical Climatology Network (USHCN). As part of our analysis, we use USHCN stations that have been corrected for Time-of-Observation Biases (Karl et al. 1986) but have not undergone the additional Menne & Williams (2009) homogenization process.We downloaded this intermediate USHCN dataset from https://www. ncdc.noaa.gov/ushcn/introduction [version 2.5.5.20190912, accessed 12/09/2019]. In this study, we are considering the annual average temperature (as opposed to monthly or seasonal averages), and so we introduce the requirement of 12 complete months of data for a given year. Also, to define how the annual average temperatures vary for a given station, we adopt the popular approach of converting each temperature record into an "anomaly time series" relative to a constant baseline period of 1961-1990 (which is the 30 year baseline period with maximum station coverage). After our final hemispheric series have been constructed, they are rescaled from the 1961-1990 baseline period to a 1901-2000 baseline period, i.e., relative to the 20^th-century average. We require a station to have a minimum of at least 15 complete years of data during this 1961-1990 time period to be incorporated into our analysis. These restrictions, although relatively modest, reduce the total number of stations in the dataset from 7 280 to 4 822.

Version 3 of the GHCN dataset includes two different estimates of how urbanized each station is. Each station is assigned a flag with one of three possible values ("Rural", "Semi-urban" or "Urban") depending on the approximate population associated with the station location. The stations are also assigned a second flag (again with one of three possible values) according to the brightness of the average night-light intensity associated with the station location – see Peterson et al. (1999) for details.

In isolation, both flags are rather crude and somewhat out-dated estimates of the station's degree of urbanization. Nonetheless, stations categorized as "rural" according to both flags are usually relatively rural, while those categorized as "urban" according to both flags are usually highly urbanized.

Therefore, by using both flags together, it is possible to identify most of the highly rural and highly urban stations. As mentioned above, some of us are currently developing more sophisticated and more up-to-date estimates of urbanization for version 4 of the GHCN dataset and the related ISTI dataset, e.g., Soon et al. (2018). However, for this study (which utilized GHCN version 3), we follow the approach proposed in Soon et al. (2015) and divide the stations into three categories:

• –  
  "Rural" = rural according to both flags (1278 of the 4822 stations, i.e., 27%)
• –  
  "Urban" = urban according to both flags (1129 of the 4822 stations, i.e., 23%)
• –  
  "Semi-urban" = the remainder of the stations (2415 of the 4822 stations, i.e., 50%).

In this paper, we are studying the Northern Hemisphere, which has much greater data coverage. However, for context, in Figure 6, we also compare the data availability for both hemispheres. The total number of stations (in either hemisphere) available for each year is plotted in Figure 6(a). The total number of stations that meet the "Rural" requirements are plotted in Figure 6(b). It can be seen that both plots reach a maximum during the 1961–1990 period, but the available data are fewer outside this period, and much lower for the early 20^th century and earlier. This has already been noted before, e.g., Lawrimore et al. (2011). However, as Soon et al. (2015) pointed out, the problem is exacerbated when we consider the rural subset, i.e., Figure 6(b). While 27% of the GHCN stations with at least 15 complete years in the optimum 1961-1990 baseline period are rural, most of these rural stations tend to have relatively short and/or incomplete station records.

Zoom In Zoom Out Reset image size

These findings are not surprising since it is more challenging to staff and maintain a long and continuous multidecadal record for an isolated rural location than for a well-populated urban location (Soon et al. 2015; Soon et al. 2018; Connolly & Connolly 2014c; Ren & Ren 2011; Ren et al. 2015). This was especially the case before the relatively recent invention of automated weather stations. As a result, the percentage of available stations that are rural is much lower for the earlier periods. For instance, only 454 of the 2163 stations (21%) with data for 1931 and only 300 of the 1665 stations (18%) with data for 1901 are rural, according to our categorization.

A further difficulty is that many of the most rural stations with relatively long station records do not have complete records. The station records might nominally cover a relatively long period, but there will often be large gaps of several years or even decades, and the average annual temperature of the station will often be quite different before and after the gap (Connolly & Connolly 2014a,c,b). For instance, using the same categorization as here, there are eight rural stations in Version 3 of the GHCN dataset for India and all of them cover a relatively long period (i.e., five have at least some data for the late 19^th century, and two of the station records begin in the first decade of the 20^th century) (Connolly & Connolly 2014b). However, most of the records contain substantial data gaps – often accompanied by substantial "jumps" in the average annual temperature that suggest a station move or some other non-climatic bias. None of the records is complete enough to continuously describe temperature trends from the late 19^th century to present.

If station histories (often called "station metadata") that indicate documented changes in station location, instrumentation, time of observation, etc., exist then, in many cases, it may be possible to correct the station records for these non-climatic biases – especially if parallel measurements associated with the station change are available (as is becoming more common in recent years). However, unfortunately, station histories are not currently provided with the GHCN dataset.

As mentioned earlier, several groups have argued that by applying statistically-based "homogenization" techniques to the dataset, the homogenization algorithms will accurately detect and correct for the main biases (Peterson et al. 1999; Menne & Williams 2009; Hausfather et al. 2013; Li & Yang 2019; Li et al. 2020b). However, when applied to rural stations utilizing urban neighbors, these techniques are prone to "urban blending", which has a tendency to "alias" urbanization bias onto the rural station records. That is, the homogenization process can contaminate the station records of nominally "rural" stations even if they had been unaffected by urbanization bias before homogenization (Soon et al. 2015; Soon et al. 2018, 2019b; Connolly & Connolly 2014a,c,b).

For those reasons, Soon et al. (2015) constructed their estimate, relying on the non-homogenized dataset, but only rural (or mostly rural) stations taken from four regions with a high density of rural stations and/or where they had relevant station history information. Here, we will adopt this approach, but updating and slightly modifying the four regions as follows:

• –  
  Rural Arctic. The gridded mean average of all 110 rural GHCN stations north of 60°N. Soon et al. (2015) limited their analysis to the Arctic Circle (i.e., north of 66.7°N), which included fewer stations, and, therefore, they also included some urban and semi-urban stations. However, following Connolly et al. (2017), we include only rural stations and have expanded the "Arctic" region to include any stations north of 60°N.
• –  
  Rural Ireland. Soon et al. (2015) demonstrated that the adjustments applied by NOAA to the longest rural record for Ireland, Valentia Observatory, were remarkably inconsistent and also failed to identify the –0.3° C cooling bias introduced by the station's February 2001 move. For this reason, we use the unadjusted rural records for Ireland and manually apply a correction to account for this non-climatic bias. The Soon et al. (2015) composite regarded this manually homogenized Valentia Observatory record as representative of rural Ireland. However, they noted that the other four rural Irish stations in the GHCN dataset also implied similar trends over their period of overlap. Therefore, in this study, we use the gridded mean of all five rural Irish stations.
• –  
  Rural United States. Our rural US series is the same as in Soon et al. (2015), except updated to 2018.
• –  
  Rural China. Our Chinese series is also the same as in Soon et al. (2015), except updated to 2018. As there are very few rural stations in China with records covering more than ∼70 years, this component of our analysis does include some stations that are currently urban or semi-urban, but whose station records have been explicitly adjusted to correct for urbanization bias as described in Soon et al. (2015).

The locations of the stations in the four regions are displayed in Figure 6(c). The regions are distributed at different locations in the Northern Hemisphere and include sub-tropical regions (lower United States and China) as well as polar regions (the Arctic). However, none of the regions are in the Southern Hemisphere; thus, this is purely a Northern Hemisphere estimate.

If accurate estimates of the urbanization bias associated with individual station records could be determined and corrected for (as for the Chinese subset), then in principle, this analysis could be expanded to include more of the available GHCN data. This could be particularly promising for some of the longer station records, especially if they have only been modestly affected by urbanization bias. For instance, Coughlin & Butler (1998) estimated that the total urbanization bias at the long (1796-present) Armagh Observatory station in Northern Ireland (UK) was probably still less than 0:2° C by 1996 (Coughlin & Butler 1998). Similarly, Moberg & Bergström (1997) were able to develop urbanization bias corrections for two long records in Sweden (Uppsala and Stockholm) (Moberg & Bergström 1997). Given the variability in urbanization biases between stations, we suggest that attempts to correct records for urbanization bias may require case-by-case assessments along these lines, rather than the automated statistical homogenization techniques currently favored. We encourage more research into expanding the available data in this way, and some of us are currently working on doing so for some regions with relatively high densities of stations. However, in the meantime, we suggest that the approach described in this paper of only considering rural data is an important first step in overcoming the urbanization bias problem.

As can be seen from Figure 6(d)-(f), the four regions alone account for more than 80% of the rural data for the early 20^th century from either hemisphere. Therefore, this estimate is more likely to be representative of rural Northern Hemisphere land surface temperature trends than the standard estimates described in Section 3.2 since, in those estimates, most of the additional data comes from stations that are more urbanized. Nonetheless, because the stations are confined to four regions rather than distributed evenly throughout the hemisphere, it is unclear what is the most suitable method for weighting the data from each station. In Soon et al. (2015), weighting was carried out in two stages. For each region, the stations were assigned to 5° × 5° grid boxes according to latitude and longitude. For each year where data for at least one station exist for a grid box, the temperature anomaly for that grid box was the mean of all the temperature anomalies of the stations in that grid box. However, the surface area of a grid box decreases with latitude according to the cosine of the latitude. Therefore, when calculating the regional average for a given year, the grid boxes with data were weighted by the cosine of the mid-latitude for that box.

Soon et al. (2015) argued that each of the four regions sampled a different part of the Northern Hemisphere and therefore the average of all four regions was more representative of hemispheric trends than any of the individual regional estimates. However, since the samples covered different regions of the Northern Hemisphere, they also weighted the four regions according to the cosine of the mid-latitude of the region. They also confined their time series to the period from 1881 onward, i.e., when all four regions had data. The updated version of this weighting approach is shown in Figure 7(a).

Zoom In Zoom Out Reset image size

Other approaches for weighting the data could be followed instead. One approach would be to give all four regions equal weighting. The results from this "equal weighting scheme" are depicted in Figure 7(c), and we extend this series back to the earliest year in which we have data, i.e., 1841. Another approach is perhaps more conventional – instead of calculating each regional estimate and averaging them together, the hemispheric averages are calculated directly from all available gridded averages (from all four regions) for each year. This is the "standard weighting scheme" shown in Figure 7(b).

It might be argued that one of our four regions could be unusual for some reason and, thereby, might be unrepresentative of the hemispheric trends. For that reason, we also calculate four additional versions of the "standard weighting scheme" in which we only use three of the four regions. These four estimates are shown in Figure 7(d)-(g).

Slight differences exist between each of the estimates. Therefore, for our analysis in this paper, we utilize the mean as well as the upper and lower bounds of all seven estimates. These bounds are calculated as the mean of the seven series for each year ±2σ. Nonetheless, all seven estimates are broadly similar to each other, implying that the effects of different weighting schemes are relatively minor. Specifically, all estimates suggest that the Northern Hemisphere warmed from the late-19^th century to the mid-1940s; cooled to the mid-1970s; and then warmed until present. According to the longer estimates, a relatively warm period also existed in the mid-19^th century. Therefore, according to the rural data, the current warm period is comparable to the earlier warm periods in the mid-1940s and possibly mid-19^th century. Although the long-term linear trend is of warming, i.e., there has been global warming since the late-19^th century, this seems to be mainly because the late-19^th century was relatively cold. This is quite different from the estimates that IPCC (2013) considered, although interestingly it is quite similar to Lansner & Pepke Pedersen (2018)'s analysis based on "sheltered" stations (Lansner & Pepke Pedersen 2018). Such a result might be surprising since many studies have claimed that urbanization bias is not a problem for estimating Northern Hemisphere air temperature trends and that all land-based global temperature estimates imply almost identical results (Lugina et al. 2006; Hansen et al. 2010; Lenssen et al. 2019; Lawrimore et al. 2011; Muller et al. 2014; Jones et al. 2012; Cowtan & Way 2014; Xu et al. 2018; Jones et al. 1990; Parker 2006; Wickham et al. 2013; Peterson et al. 1999; Hausfather et al. 2013; Li & Yang 2019; Li et al. 2020b; Menne et al. 2018; Jones 2016; Hawkins & Jones 2013). Therefore, in the next subsection, possible reasons for these differences will be assessed.

Peterson et al. (1999) compared trends from a subset of only rural stations with the full GHCN dataset, but they argued that the trends were equivalent and concluded that urbanization bias was a negligible problem. However, the Peterson et al. (1999) study relied on a version of the GHCN that had been adjusted utilizing an automated statistically-based homogenization algorithm in an attempt to reduce the effects of non-climatic biases. On the other hand, Soon et al. have argued that not only do the current homogenization algorithms perform poorly when applied to the non-climatic biases in the GHCN, but they also inadvertently lead to a blending of the non-climatic biases in different stations (Soon et al. 2015; Soon et al. 2018, 2019b). Hence, after homogenization, urban blending would have transferred much of the urbanization bias of the neighboring urban stations into the "rural" stations used by Peterson et al. (1999).

To identify the reasons for the differences between our rural-only estimates and the standard estimates, in Figure 8, we compare our rural-only estimate, Figure 8(a), to three alternative estimates that are comparable to the standard estimates. The estimate in Figure 8(b) was calculated considering all stations (rural, semi-urban or urban) for the same four regions and relying on the version of the GHCN dataset that has been homogenized using the automated Menne & Williams (2009) homogenization algorithm (Menne & Williams 2009). As in our ruralonly estimate, USHCN stations are also corrected for changes in time-of-observation in the homogenized GHCN dataset (Karl et al. 1986). However, no additional attempts to correct for non-climatic biases are applied other than NOAA's Menne & Williams (2009) homogenization. Note that this series has slightly more rural stations for the United States component than our rural-only series (295 instead of 272) because there are some non-USHCN stations for the contiguous US in the GHCN dataset.

Zoom In Zoom Out Reset image size

The estimate in Figure 8(c) was calculated incorporating all Northern Hemisphere GHCN stations (rural, semi-urban or urban), and only using the Menne & Williams (2009) homogenized GHCN dataset. Finally, the estimate in Figure 8(d) is calculated as for Figure 8(c), except it is a global estimate including all stations from either hemisphere.

Let us now consider five plausible objections which might be raised as to the reliability of the rural-only series relative to the standard estimates:

• 1.  
  Could the four regions be unrepresentative of global trends? Comparing Figure 8(a) to Figure 8(c) and Figure 8(d) shows that the rural-only estimate implies much less warming than the standard estimates. Specifically, the rural-only estimate implies a long-term 1841–2018 linear trend of only 0.41° C century⁻¹ for the Northern Hemisphere, which is less than half the equivalent trends following the standard approach, i.e., 0.86° C century⁻¹ for the Northern Hemisphere, Figure 8(c), and 0.83° C century⁻¹ globally, Figure 8(d). One possibility might be that these four regions are coincidentally underestimating the warming for the rest of the Northern Hemisphere (and globe). However, as can be seen from Figure 8(b), when the standard approach is applied to those four regions, the linear trend is even greater (0.94° C century⁻¹). In other words, if anything, the four regions are ones that slightly overestimate the warming of the rest of the world.
• 2.  
  Why was an automated statistically-based homogenization algorithm like the standard estimates not used to reduce the effects of non-climatic biases? Various non-climatic biases in the data need to be corrected. However, as described earlier, Soon et al. have demonstrated that such algorithms perform poorly when applied to the non-climatic biases in the GHCN and frequently result in inappropriate adjustments to the data (Soon et al. 2015; Soon et al. 2018, 2019b). Moreover, they showed that the algorithms inadvertently lead to a blending of the non-climatic biases in different stations. Therefore, we argue that empirically-based homogenization adjustments (ideally using station history information) are more reliable than the statistically-based homogenization adjustments currently utilized by most groups, e.g., Refs. (Hansen et al. 2010; Lenssen et al. 2019; Lawrimore et al. 2011; Jones et al. 2012; Cowtan & Way 2014; Xu et al. 2018; Menne et al. 2018; Hawkins & Jones 2013).
• 3.  
  Why were the rural stations from outside these four regions not included? As can be ascertained from Figure 6, there are quite a few rural stations in the GHCN dataset that were outside of the four regions. Specifically, there are 503 rural stations in our four regions, but 455 rural stations in the rest of the Northern Hemisphere and 269 in the Southern Hemisphere. Therefore, initially it might appear that our four regions only considered 50% of the total Northern Hemisphere rural data and 39% of the global rural data. However, as can be seen from Figure 6(d)-(f), our four regions comprise the vast majority of the available rural data with relatively long station records. For instance, the four regions account for more than 80% of the rural stations from either hemisphere with early 20^th century records. This is not accidental, as Soon et al. (2015) had specifically identified those four regions as being ones that contained a relatively high density of rural stations or for which they had station history information for applying empirically-based homogenization adjustments (Soon et al. 2015). Therefore, while it might be surprising to many readers, these four regions account for the vast majority of the long-term rural station records.
• 4.  
  Were too few stations used? The rural-only estimate was constructed from a total of 554 stations (this figure is slightly higher than the 503 mentioned above, due to the non-rural Chinese stations utilized for some of the Chinese component), whereas the Northern Hemisphere estimate in Figure 8(c) relied on 4176 stations and the global estimate in Figure 8(d) incorporated 4822 stations. Moreover, several of the estimates by other groups purport to use even more stations. Indeed, the Berkeley Earth group claims to use ∼40 000 stations (Muller et al. 2014) (although most of these stations have station records with fewer than 30 years of data). Therefore, some might argue that the rural-only series differs from the standard estimates because the sample size was too small. However, as discussed in Soon et al. (2015) (see in particular their table 6), and argued by Hawkins & Jones (2013), many of the earlier attempts to estimate Northern Hemisphere temperature trends utilized even fewer stations, yet obtained fairly similar results to the latest estimates following the standard approach. In particular, Mitchell (1961) relied on only 119 stations; Callendar (1961) utilized only ∼ 350; and Lugina et al. (2006) used only 384. Moreover, Jones et al. (1997) calculated that about 50 well-distributed stations should be adequate for determining annual global air temperature trends. Since our composite is for the Northern Hemisphere only, even fewer stations should be needed. Therefore, the differences in trends between the rural-only series and the standard estimates are not a result of the smaller number of stations.
• 5.  
  Why are the differences most pronounced for the earlier period? Given that the rate of urbanization has accelerated in recent decades, initially it might be assumed that the differences between urban and rural time series would be greatest for recent decades. However, we remind readers of two factors:
• (i)  
  All the time series in Figure 8 are relative to the 20^th-century average. This partially reduces the apparent differences between the time series during the 20^th century.
• (ii)  
  From Figure 6, it can be seen that the GHCN dataset has a relatively high rural composition during the 1951–1990 period. The shortage of rural stations in the standard estimates is most pronounced before (and, to a lesser extent, after) this period. Therefore, counterintuitively, the differences between "rural only" and "urban and rural" tend to be larger before 1951, due to changes in data availability

Nonetheless, some readers might still prefer the standard estimates that incorporate both urban and rural stations. We will discuss these in the next section.

Figure 9 compares the Northern Hemisphere land-only air temperature trends of seven different estimates calculated following the standard approach, i.e., including both urban and rural stations but applying homogenization procedures such as the Menne & Williams (2009) automated statistical homogenization algorithm or Berkeley Earth's "scalpel" procedure (Muller et al. 2014). Although slight differences exist between the approaches taken by each group, the results are remarkably similar. For comparison, the "all Northern Hemisphere stations" estimate of Figure 8(c) is replicated here in Figure 9(a), albeit extended back to 1800, and it is very similar to the other six estimates. The striking similarity of all these estimates has been noted by most of the groups when describing their individual estimates, e.g., CRUTEM (Jones et al. 2012), i.e., Figure 9(b); Cowtan and Way (Cowtan & Way 2014), i.e., Figure 9(c); NOAA NCEI (Lawrimore et al. 2011), i.e., Figure 9(d); NASA GISS (Hansen et al. 2010; Lenssen et al. 2019), i.e., Figure 9(e); Berkeley Earth (Muller et al. 2014), i.e., Figure 9(f); and the Chinese Meteorological Administration (Xu et al. 2018), i.e., Figure 9(g).

Zoom In Zoom Out Reset image size

The fact that all of these estimates are so similar appears to have led to a lot of confidence within the community that these estimates are very reliable and accurate. As discussed in the previous section, we believe this confidence is unjustified and that the new rural-only estimate is more reliable. Nonetheless, for the sake of argument, and because we appreciate that many readers may disagree with us, we will also carry out our analysis for these "urban and rural"-based estimates.

Although the seven estimates are all remarkably similar, as noted above, slight differences still exist between each of the estimates. Therefore, as before, we use the mean with the upper and lower bounds of all seven estimates. These bounds are calculated as the mean of the seven series for each year ±2σ – see Figure 9(h).

One way to potentially bypass the debate over urbanization bias might be to consider ocean surface temperature trends instead of land surface temperature trends. Unfortunately, considerable uncertainties are also associated with the ocean surface temperature data, especially before the 1950s. Before the International Geophysical Year of 1957/1958, available data are quite sparse for the Northern Hemisphere and very sparse for the Southern Hemisphere (Kennedy et al. 2019). Moreover, the methods and instrumentation that were utilized for measuring ocean surface temperatures varied from ship to ship and are often poorly documented (Kent et al. 2017). For this reason, it has been argued that the ocean surface temperature data are less reliable than the land surface temperature data, e.g., Jones (2016). Nonetheless, the ocean data are largely independent of the land surface temperature data – even if some groups have partially reduced this independence by adjusting the ocean surface temperature data to better match the land surface temperature records, e.g., Cowtan et al. (2018).

Two competing sets of measurements for the ocean surface temperature exist. The Marine Air Temperatures (MATs) are based on the average air temperature recorded on the decks of ships. The SSTs are based on the average water temperature near the surface. In some sense, the MAT is more comparable to the Land Surface Air Temperature since they are both measurements of the air temperature near the surface. However, the available MAT records are more limited than the SST records, and it has been suggested that the MATs are more likely to be affected by non-climatic biases, especially the daytime measurements, e.g., Rayner et al. (2003); Kent et al. (2013).

For this reason, most of the analysis of ocean surface temperature trends has focused on the SST data. For brevity, our ocean surface temperature analysis will be confined to the SST data, but we encourage further research examining both types of dataset and note that they can yield different results (Rubino et al. 2020). Yet, even with the SST data, it is widely acknowledged that nonclimatic biases are likely to exist in the data, especially for the pre-1950s period (Jones 2016; Kennedy et al. 2019; Kent et al. 2017; Cowtan et al. 2018; Kennedy 2014; Davis et al. 2019).

It is likely that if the true magnitudes and signs of these biases could be satisfactorily resolved, as Kent et al. (2017) have called for, this could also help us to resolve the debates over the land surface temperature data. For instance, Davis et al. (2019) noted that if you separately analyze the SST measurements that were taken via "engine-room intake" and those taken by "bucket" measurements, this implies quite different global temperature trends over the period 1950-1975. Specifically, the "engine-room intake" measurements imply a substantial global cooling trend from 1950 to 1975, commensurate with our rural estimates for land surface temperatures in Section 3.1. Interestingly, the "sheltered" stations subset of Lansner & Pepke Pedersen (2018) suggests the same. By contrast, the "bucket" measurements imply a slight increase over this same period – consistent with our "urban and rural" estimates described in Section 3.2 and also the "non-sheltered" subset of Lansner & Pepke Pedersen (2018).

Currently, three different versions of SSTs are available from the Hadley Centre (HadISST (Rayner et al. 2003), HadSST3 (Kennedy 2014; Kennedy et al. 2011b,a) and HadSST4 (Kennedy et al. 2019)) and three from NOAA (ERSST v3 (Smith et al. 2008), v4 (Huang et al. 2015) and v5 (Huang et al. 2017)). The Japanese Meteorological Agency also has a version (COBE SST2 (Hirahara et al. 2014)), but here, our analysis will be limited to the more widely-used NOAA and Hadley Centre datasets.

Figure 10 compares the various Northern Hemisphere SST estimates. In recognition of the considerable uncertainties associated with estimating sea surface temperature trends, the Hadley Centre has taken to providing ensembles of multiple different plausible estimates (100 each) for their two most recent versions, i.e., HadSST3 (Kennedy 2014; Kennedy et al. 2011b,a) and HadSST4 (Kennedy et al. 2019). However, the previous version, HadISST (Rayner et al. 2003), and the three NOAA versions did not take this approach. Therefore, if we were to treat all HadSST3 and HadSST4 ensemble members as separate estimates, then the other four datasets would not make much contribution to our analysis. On the other hand, if we only considered the "median" estimates from the HadSST3 and HadSST4 datasets, then we would be underestimating the uncertainties associated with these datasets by the Hadley Centre. With that in mind, for both of these datasets, we have extracted the "lower realization" and "upper realization", i.e., the lower and upper bounds, and treated these as four separate estimates that are representative of the uncertainties implicit within the combined 200 "ensemble" estimates of the two datasets – see Figure 10(b)-(e). We then treat the Hadley Centre's other version, HadISST – see Figure 10(a) - and the three NOAA estimates – see Figure 10(f)-(h) – as an additional four distinct estimates.

Zoom In Zoom Out Reset image size

This provides us with a total of eight different estimates. For our SST analysis, we use the mean and the upper and lower bounds of all eight estimates. These bounds are calculated as the mean of the eight series for each year ±2σ. See Figure 10(i).

To assess the role of the Sun in Northern Hemisphere temperature trends before the mid-19^th century, analysis must be mostly limited to indirect estimates from paleoclimate reconstructions developed utilizing temperature proxy series, e.g., tree rings (Esper et al. 2018).

Esper et al. (2018) recently compiled and reviewed six different tree-ring based millennial-length reconstructions for the Northern Hemisphere. Figure 11 compares the six different Northern Hemisphere tree-ring temperature proxy-based reconstructions described by Esper et al. (2018). Although all these reconstructions share some of the same underlying tree-ring data, and they imply broadly similar long-term temperature trends, Esper et al. (2018) stressed that each of the reconstructions gave slightly different results. They emphasized that it was inappropriate to simply average the reconstructions together – a point which St. George & Esper (2019) reiterated. Therefore, here, each of these six reconstructions should be treated as a different estimate. However, for consistency with the rest of our estimates, our analysis will be based on the mean along with the upper and lower bounds, which are calculated as the mean of the six series for each year ±2σ. See Figure 11(g) for the entire period spanned by these estimates, which begin in the 9^th century, and Figure 11(h) for a close-up of the estimates since the start of the 19^th century.

Zoom In Zoom Out Reset image size

Leclercq & Oerlemans 2012 constructed a dataset of all the available multidecadal glacier length records for both hemispheres. They then regarded the changes in glacier length at each of these glaciers as a proxy for local temperature. Leclercq & Oerlemans nominally provided separate estimates for global temperatures as well as each hemisphere. However, as can be seen from Figure 12(b), most of the data was for the Northern Hemisphere, which is the hemisphere our analysis focuses on. It can also be seen that most of the glacier length records began in the 20^th century.

Zoom In Zoom Out Reset image size

Nonetheless, they had some data for both hemispheres to cover the entire period 1600–2000, although the error bars provided by Leclercq & Oerlemans (2012) are quite substantial, especially for the pre-20^th century period. We will treat the mean with upper bound and lower bound from the Leclercq & Oerlemans (2012) Northern Hemisphere reconstruction as our estimates – see Figure 12(a).

In Figure 13, we compare and contrast the five different estimates of Northern Hemisphere temperature trends since the 19^th century (we only present the post-18^th century data for our two longer proxy-based estimates in this figure for clarity). Most of the estimates share key similarities, specifically,

• –  
  They all imply that current temperatures are higher than at the end of the 19^th century, i.e., that some warming has occurred since the end of the 19^th century.
• –  
  There was an "Early 20^th Century Warm Period" (ECWP) in the mid-20^th century - peaking sometime in the 1930s–1950s – followed by several decades of either cooling or a lack of warming.
• –  
  There was a definite period of warming from the late-1970s to the end of the 20^th century. Some estimates suggest that this warming has continued to the present day, although there is some debate in the literature over whether there has been a "warming hiatus" after the end of the 20^th century and if so whether this hiatus has ended.

Zoom In Zoom Out Reset image size

However, subtle differences exist between the various estimates, including,

• –  
  The timings of the various warming and cooling periods, e.g., most of the estimates imply the ECWP peaked during the 1940s, while the tree ring-based proxy estimates imply a later peak during the 1950s.
• –  
  The exact magnitudes of the various warming and cooling periods, e.g., the estimates using both urban and rural stations, imply a larger long-term warming trend and that the warming has been almost entirely continuous, while the other estimates imply a more nuanced alternation between multidecadal warming and cooling periods.
• –  
  The two proxy-based estimates imply a more dampened temperature multidecadal variability than the instrumental estimates. This is not surprising considering that they are based on indirect estimates of temperature variability, rather than direct temperature measurements. For instance, while glacier lengths are certainly influenced by the local air temperature during melt season (i.e., summer), they are also influenced by winter precipitation (i.e., how much snow accumulates during the growth season). Indeed, Roe and O'Neal (2009) note that glaciers can, in principle, advance or retreat by kilometers without any significant climate change, but simply due to year-to-year variability in local weather (Roe & O'Neal 2009). Similarly, García-Suárez et al. (2009) argue that tree-ring proxies are often proxies for multiple different climatic variables as well as temperature, e.g., sunshine, precipitation, soil moisture, etc., and also mostly reflect conditions during the growing season (García-Suárez et al. 2009). Therefore it can be difficult to extract a pure "temperature" signal (García-Suárez et al. 2009). Meanwhile, Loehle (2009) notes that even with tree-ring proxies from "temperaturelimited regions" (usually either high latitude or high elevation sites), the growth rates may still be influenced by other factors, and that the relationship between growing season temperatures and tree-ring growth may be nonlinear (Loehle 2009).

It is plausible that some (or even all) of these differences arise from non-climatic biases in one or more (and perhaps all) of the time series. That said, some (or even all) may arise because each time series is capturing a slightly different aspect of the true climatic variability, e.g., the oceanic surface temperature variability may be slightly different from the land surface temperature variability.

At any rate, these nuanced differences are important to consider for our analysis, because we will be estimating the influence of the Sun's variability on Northern Hemisphere temperature trends using a relatively simple linear least-squares fit between a given temperature estimate and each of our 16 estimates of solar variability, which we discussed in Section 2.4. The successes of the fits will depend largely on the linear correlations between the timings and magnitudes of the various rises and falls in the time series being fitted. For this reason, we propose to analyze the different categories of Northern Hemisphere temperature variability estimates separately. Our fits will be repeated separately for each of the 16 solar variability estimates to both the upper and lower bounds. This then provides us with an uncertainty range for each of our estimates.