Holocene Millennial-Scale Solar Variability and the Climatic Responses on Earth

: The solar impact on Earth’s climate is both a rich and open-ended topic with intense debates. In this study, we use the reconstructed data available to investigate periodicities of solar variability (i.e., variations of sunspot numbers) and temperature changes (10 sites spread all over the Earth) as well as the statistical inter-relations between them on the millennial scale during the past 8640 years (BC 6755–AD 1885) before the modern industrial era. We ﬁnd that the variations of the Earth’s temperatures show evidence for the Eddy cycle component, i.e., the 1000-year cyclicity, which was discovered in variations of sunspot numbers and believed to be an intrinsic periodicity of solar variability. Further wavelet time-frequency analysis demonstrates that the co-variation between the millennium cycle components of solar variability and the temperature change held stable and statistically strong for ﬁve out of these 10 sites during our study interval. In addition, the Earth’s climatic response to solar forcing could be different region-by-region, and the temperatures in the southern hemisphere seemed to have an opposite changing trend compared to those in the northern hemisphere on this millennial scale. These ﬁndings reveal not only a pronounced but also a complex relationship between solar variability and climatic change on Earth on the millennial timescale. More data are needed to further verify these preliminary results in the future.


Introduction
The question of how strongly the Sun influences climate change on Earth is eternally debated, with topics ranging from the lack of evidence to contradictions. The debate involves another bigger and hotter topic, i.e., global warming, in the modern industrial period. It has been concluded that human influence, via emissions of greenhouse gases, is likely the dominant cause of the observed warming since the mid-20th century [1]. Yet it is equally clear that natural forces are still a non-negligible factor in driving climate change [2], involving co-modulating factors from both the intrinsic changes in solar irradiance and the Sun-Earth orbital conditions [3]. For example, the Earth's planetary temperature experienced some periods in recent history that were probably as warm as the present, such as the Medieval Warm Period (AD 800-AD 1300), and other periods colder than present, such as the Little Ice Age (AD 1300-AD 1850). This proves that natural forces can drive climate change of large amplitude even without a powerful human influence. Solar forcing is one important ingredient of these natural driving factors [4][5][6]. Early in 1801, Sir William Herschel found that the rainfall reduced when there were few sunspots on the Sun, and pointed out that solar variability might play a role in the change of Earth's climate [7,8]. The Maunder Minimum of anomalously low and weak sunspot activity (AD 1645-AD 1715) coincided in time with the coldest period of the Little Ice Age in many places of the world [9,10], as well as with that of the historical Ming-Qing Little Ice Age in China [11,12].
Sunspot number (SSN) is the most widely used parameter to qualitatively represent the timing and strength of solar magnetic activity, and this unique solar index also has the longest record time. Direct observations of sunspots started from the invention of the astronomical telescope, and, so far, scientists have accumulated only 400 years of data [13]. The available observational data would be insufficient in investigations exceeding centennial scales, such as millennial or even longer timescales. Fortunately, the proxy reconstructed data can be a reasonable substitute. Similarly, temperature is the most fundamental climatic index with long recording while the available instrumental thermometer data would be limited to studies of short timescales. There are many recent Holocene data sets reconstructed from various paleoclimatic proxies so far that have better resolution and chronology than earlier works [14]. For example, Marcott et al. [15] provided a wider view by reconstructing both regional and global temperatures for the past 11,300 years from 73 globally distributed sites. We found some clues about the evidence of solar modulation on the Earth's climate change on a millennial scale in our last work [16], but the adopted climatic indicators were limited only to four local places on Earth (i.e., Greenland, Lake Qinghai and Dongge cave in China, Vostok in Antarctica). In this paper, we will explore the relations between solar variability and temperature changes for a more widely distributed database on Earth in order to check and expand our previous conclusions. This is an expanded follow-up of our last work [16].

Materials and Methods
The reconstructions of SSN are usually based on individual, sometimes statistically averaged, proxy data sets. Different from these, Wu et al. [17] combined the long-span data sets of 10 Be and 14 C in the available terrestrial archives to reconstruct, for the first time, a new consistent SSN data set. It covers 8640 years in the past (from BC 6755 to AD 1885) with the fixed time resolution of 10 years. We obtained this SSN data from their website (http://www2.mps.mpg.de/projects/sun-climate/data/SN_composite.txt).
For these 73 globally distributed temperature records used in Marcott et al. [16], we selected data for our analysis based on the following criteria: (1) covering the period BC 6755-AD 1885; (2) the time resolution higher than or about 100 years; (3) temperature coming from a single site, not composite ones. In this way we obtained temperature data of 10 sites (T 1~T10 ) distributed from the Arctic to the Antarctic, i.e., temperatures of Agassiz-Renland, Tsuolbmajavri Lake, RAPID-12-1K, OCE326-GGC30, MD98-2181, ODP_1084B, TN057-17, Homestead Scarp, Dome C, and Vostok.
The reconstructed temperature of Agassiz-Renland (T 1 ) was based on an average of uplift corrected δ 18 O data in the ice sheet from Agassiz and Renland [18]. The July temperature of Tsuolbmajavri Lake (T 2 ) was reconstructed from pollen assemblages preserved in a sediment core [19]. The reconstruction of temperature of the Atlantic inflow RAPID-12-1K (T 3 ) was based on paired Mg/Ca-δ 18 O measurements based on two species of planktonic foraminifera, Globigerina bulloides and Globorotalia inflata [20]. The alkenone-derived sea surface temperature (SST) of the slope water core OCE326-GGC30 (margin east of Nova Scotia basin, T 4 ) was determined by alkenone paleothermometry [21]. The sea temperature of the marine sediment core MD98-2181 from the Morotai Basin (T 5 ) was obtained from planktonic foraminiferal stable isotope and Mg/Ca records in the core [22].
The Mg/Ca data from ODP_1084B (Ocean Drilling Program Leg 175 Hole 1084B) in the Benguela coastal upwelling system recorded the oceanic temperature (T 6 ) [23]. The summer SST at site TN057-17 near the Polar Front in the East Atlantic Southern Ocean (T 7 ) was obtained from the down-core variability of diatom assemblages [24]. The summer temperature in the low forest at Homestead Scarp (T 8 ) (in Campbell Island in the Southern Ocean) was reconstructed from the pollen peat sequence at temperature-sensitive locations within the tree-line ecotone [25]. The temperature in the European Project for Ice Coring in Antarctica Dome C ice core (T 9 ) was obtained from a high-resolution deuterium profile in the core together with experiments performed with an atmospheric general circulation model including water isotopes [26]. Similarly, a data series of the local temperature at Vostok in East Antarctica (T 10 ) was reconstructed for the past 420,000 years based on the deuterium content of the ice (δD ice ) in Vostok's ice cores [27].
The reconstruction mechanisms and paleo-source media are greatly different for these 10 temperatures. Detailed information can be found in the related references, i.e., [18][19][20][21][22][23][24][25][26][27]. Table 1 lists the geographic information for their sites. Their locations on the Earth's surface are shown in Figure 1. We can see that these 10 sites are evenly distributed in the northern and southern hemispheres, i.e., four sites located in the northern hemisphere (T 1 , T 2 , T 3 , T 4 ), four sites in the southern hemisphere (T 7 , T 8 , T 9 , T 10 ), and two sites in the tropical region on opposite sides of the equator (T 5 , T 6 ), with similar latitudinal intervals between two neighboring sites. The averaged time resolution of each temperature is also given in Table 1. In order to align with SSN data, we computed the interpolations of T 1~T10 from BC 6755 to AD 1885 at fixed time resolutions of 10 years, and adopted them as the temperature data. mer temperature in the low forest at Homestead Scarp (T8) (in Campbell Island in the Southern Ocean) was reconstructed from the pollen peat sequence at temperature-sensitive locations within the tree-line ecotone [25]. The temperature in the European Project for Ice Coring in Antarctica Dome C ice core (T9) was obtained from a high-resolution deuterium profile in the core together with experiments performed with an atmospheric general circulation model including water isotopes [26]. Similarly, a data series of the local temperature at Vostok in East Antarctica (T10) was reconstructed for the past 420,000 years based on the deuterium content of the ice (δDice) in Vostok's ice cores [27].
The reconstruction mechanisms and paleo-source media are greatly different for these 10 temperatures. Detailed information can be found in the related references, i.e., [18] to [27]. Table 1 lists the geographic information for their sites. Their locations on the Earth's surface are shown in Figure 1. We can see that these 10 sites are evenly distributed in the northern and southern hemispheres, i.e., four sites located in the northern hemisphere (T1, T2, T3, T4), four sites in the southern hemisphere (T7, T8, T9, T10), and two sites in the tropical region on opposite sides of the equator (T5, T6), with similar latitudinal intervals between two neighboring sites. The averaged time resolution of each temperature is also given in Table 1. In order to align with SSN data, we computed the interpolations of T1~T10 from BC 6755 to AD 1885 at fixed time resolutions of 10 years, and adopted them as the temperature data.   Table 1 (see also Figure 5), and red dots Figure 1. Locations of 10 sites of the reconstructed temperature data adopted in this study. Blue dots denote locations of temperatures with positive correlations to sunspot number (SSN) reported in Table 1 (see also Figure 5), and red dots denote those with negative correlations to SSN. A rough pattern of opposing inter-hemispheric seesaw temperature responses to solar variations on the millennial timescale over the Holocene, with the one exception of T 9 , is indicated. We adopted the same analytical methods as those in our previous work [16]. We used the Lomb-Scargle periodogram method [28,29] to analyze periodicities of both solar variability and temperature change on Earth. This method has been widely used in geophysical research. In comparison with the traditional Fourier transform, this method had an evident advantage of applicability to incomplete or unevenly sampled time series.
In order to further disclose potential relations between solar variability and temperature change on Earth, we used the Wavelet Coherence (WTC) tool [30] to compute the wavelet coherence between SSN and temperatures. This wavelet coherence shows correlations of two signals in the time-frequency space and can display their relative phases, as well. Here, "frequency" can be converted to period (p) of the signal. Therefore, the corresponding wavelet coherence is expressed in the time-period space. We downloaded the software package of this tool from: http://noc.ac.uk/marine-data-products/crosswavelet-wavelet-coherence-toolbox-matlab. Figure 2 shows the Lomb-Scargle periodograms of these time series (both SSN and temperatures). The red horizontal dotted lines are the 95% confidence levels. We see that solar variability represented by SSN had a host of periodicities over the 95% confidence level, such as the 150-year cyclicity, 200-year cyclicity (Suess-de Vries cycle [31]), 350-year cyclicity, 500-year cyclicity, millennial cyclicity (990 years and 1149 years) (Eddy cycle [9,[32][33][34]). Similar results were found in our previous study [16]. On the other hand, the temperatures had relatively fewer significant periodicities. The only common cyclicity between SSN and these temperatures was the Eddy cycle, which varied from 900 years to 1100 years (shaded region in Figure 2). For SSN, T 5 and T 6 , there were two peaks over the 95% confidence level around the Eddy cycle; for other temperatures, only one peak was significant except for T 4 (its peak at 962 years was lower than the 95% confidence level). The common periodicity between SSN and temperatures indicated potential physical connections, and attracted us to make further exploration between them.

Cross Correlations at Millennial Cycle
In order to further disclose potential relations between solar variability and temperature change on Earth, we used the Wavelet Coherence (WTC) tool [30] to compute the wavelet coherence between SSN and temperatures. Figure 3 gives the squared wavelet coherence (C w ) between SSN and T 1 , T 2 , . . . , T 10 , respectively. The definition of C w can be found in Equation (8) of [30]. Most high-correlation regions appear as "isolated islands", except for the high-power region around the millennium cycle. SSN had a good correlation to temperatures at this millennial cycle. The wavelet coherences at 1000-year cyclicity were stronger for T 3 , T 4 , T 5 and T 6 than those for the temperatures at other locations. We adopted the global wavelet coherence C g (p) and angle strength S θ (p) to quantitatively estimate the cross-wavelet correlations between SSN and temperatures. They are parameters integrated along time at fixed period (p) in the wavelet coherence output by the WTC tool, and their detailed definitions can be found in [16]. Here, C g reflects an overall correlation between two series at a fixed p, and S θ depicts how strongly the phase angles between series keep stable. C g = 1 means absolute correlation, and C g = 0 means absolute non-correlation; S θ = 1 indicates absolute "phase-locked", and S θ = 0 indicates absolute "non-phase-locked." The combination of C g and S θ can reveal whether the co-variances of two analyzed series are strong and stable. We adopted the global wavelet coherence Cg(p) and angle strength Sθ(p) to quantitatively estimate the cross-wavelet correlations between SSN and temperatures. They are parameters integrated along time at fixed period (p) in the wavelet coherence output by the WTC tool, and their detailed definitions can be found in [16]. Here, Cg reflects an   Figure 4 shows the variations of C g and S θ between SSN and temperatures along p outside the COI (cone of influence). The vertical shaded region denotes the millennial cyclicity, i.e., 900 years to 1100 years. We can see that both curves (C g and S θ ) have peaks almost simultaneously at the millennial periodicity for T 1 , T 4 , T 5 , T 7 and T 10 . For example, S θ,peak (1ky) = 0.76 and C g,peak (1ky) = 0.44 for T 1 . For the other temperatures (T 2 , T 3 , T 6 , T 8 , T 9 ), their angle strengths present a peak at this millennial cycle, but the peak of their wavelet coherence is not so evident as those in angle strengths, meaning that only phases were securely locked at this millennial cycle. In general, this figure demonstrates that the wavelet coherences between SSN and temperatures kept strong and stable at the 1000-year cyclicity, relative to other cyclicities.

Solar Oscillations and Climatic Responses at Millennial Scale
For a better display of the millennial evolutions of solar variability and climate changes on Earth, we used a band-pass filter to extract millennial signals from original data. The filter was designed to be centered at 1000 years, and its two half-power points were at 1000 ± 100 years. The output signals are shown in Figure 5 (thin lines). The millennial variations of SSN had a consistent trend with those of T1, T2, T3, T4, T5, T9, and positive correlations existed between them. The linear Pearson correlation coefficients (CCs) between SSN and these temperatures ranged from 0.36 (T9) to 0.93 (T3) (see values listed in Table 1). As for T6, T7, T8, T10, their changing trends were converse to that of SSN with negative correlations; T10 had the strongest negative correlation to SSN (CC = −0.74); and T6 had the weakest negative correlations (CC = −0.11). We adopted the tool developed by Macias-Fauria et al. [35] to estimate the statistical significance of these CCs. In this tool, a combination of Monte Carlo iterations and time-series modeling in the frequency [36] or temporal domains [37] was adopted to estimate the significance of CC by comparing the correlation to the distribution of correlations between the second target series and the random-phase series, which had the same periodogram as the first

Solar Oscillations and Climatic Responses at Millennial Scale
For a better display of the millennial evolutions of solar variability and climate changes on Earth, we used a band-pass filter to extract millennial signals from original data. The filter was designed to be centered at 1000 years, and its two half-power points were at 1000 ± 100 years. The output signals are shown in Figure 5 (thin lines). The millennial variations of SSN had a consistent trend with those of T 1 , T 2 , T 3 , T 4 , T 5 , T 9 , and positive correlations existed between them. The linear Pearson correlation coefficients (CCs) between SSN and these temperatures ranged from 0.36 (T 9 ) to 0.93 (T 3 ) (see values listed in Table 1). As for T 6 , T 7 , T 8 , T 10 , their changing trends were converse to that of SSN with negative correlations; T 10 had the strongest negative correlation to SSN (CC = −0.74); and Universe 2021, 7, 36 7 of 11 T 6 had the weakest negative correlations (CC = −0.11). We adopted the tool developed by Macias-Fauria et al. [35] to estimate the statistical significance of these CCs. In this tool, a combination of Monte Carlo iterations and time-series modeling in the frequency [36] or temporal domains [37] was adopted to estimate the significance of CC by comparing the correlation to the distribution of correlations between the second target series and the random-phase series, which had the same periodogram as the first target series [36]. The CCs and their significance levels (p) between SSN and temperatures are shown in Table 2 as well. From Table 2, we find the following results for millennial variations of these series. Firstly, all temperatures in the northern hemisphere (T 1~T5 ) were positively correlated with SSN. As to temperatures in the southern hemisphere (T 6~T10 ), four of them were negatively correlated with SSN except for T 9 . The site location is marked as a blue dot in Figure 1 if the corresponding temperature was positively correlated with SSN, and marked as a red dot for the negatively correlated temperature. Secondly, we were careful to the check the significance of these correlations. The significance levels were classified into four kinds: For T 3 , p < 0.001, indicating the CC was very significant; for T 4 and T 10 , 0.1 < p < 0.2, indicating that their CCs were relatively significant; for T 1 and T 5 , 0.2 < p < 0.3, therefore, their CCs to SSN were only weakly significant. For the other temperatures (T 2 , T 6 , T 7 , T 8 , T 9 ), p > 0.3 for their CCs; we declared that their correlations to SSN were not statistically significant. Figure 5 also provides a direct comparison between SSN and temperatures for their original data (thick lines). We can see that the millennial components we derived here reflected the rough trends of their raw data on the millennial scale. The common periodicity of 1000 years (revealed in Figure 2), the generally strong and stable resonant relationships between SSN and temperatures during this periodicity (revealed in Figures 3-5) jointly indicated that solar variability may have played a role in the modulation of this millennial oscillation of the Earth's climate change, but the climatic responses on Earth are not uniform and can be different from region to region.
As far as the regional effects of climate change are concerned, an oppositely changing trend between the millennial temperatures of two hemispheres was discovered for 20 years [38][39][40][41][42][43][44], which can be termed "northern-southern hemispheric seesaw". Based on our results, all temperatures in the northern hemisphere were positively correlated with SSN (T 1 , T 3 , T 4 , T 5 ) but at different significance levels. For temperatures in the southern hemisphere, four out of five were negatively correlated with SSN, but only one (T 10 ) was significant. This indicated that temperatures in the southern hemisphere seemed to have an opposite trend in contrast with the northern hemisphere's temperatures on the millennial scale. This was consistent with previous studies on seesaw [42,45]. For example, the modeling work demonstrated that responses of regional temperatures to the solar forcing were not uniform [42]; although the millennial variations of the northern temperature kept the positive-phase consistency with the global change, a small part of the very southern regions presented a converse trend with the global variation [45].
The seesaw pattern, also revealed by other evidences, seemed to be a robust feature of the long-term climate change on Earth [40,41]. For example, Eroglu et al. [46] found a seesaw variation in the Holocene East Asian-Australian summer monsoon at millennial to sub-centennial timescales during the past 9000 years. However, our results indicated a complex seesaw pattern of the Earth's climatic responses to solar forcing, rather than a strict northern-southern seesaw: different regions within the same hemisphere could have different significant levels of variations, even an opposite changing trend. Although some regions presented evident responses to solar forcing, other regions might have been insensitive to the driving forcing. This further indicated that if we considered the regionaveraged temperatures on Earth (especially for regions in the southern hemisphere), it would be difficult to uncover the potential seesaw relation between them as the averaging process would weaken the differences between regions. This fact might explain why this millennial seesaw pattern of the Earth's climate was so highly controversial [47,48].  Table 1, Table 2, and Figure 1.
The seesaw pattern, also revealed by other evidences, seemed to be a robust feature of the long-term climate change on Earth [40,41]. For example, Eroglu et al. [46] found a seesaw variation in the Holocene East Asian-Australian summer monsoon at millennial to sub-centennial timescales during the past 9000 years. However, our results indicated a complex seesaw pattern of the Earth's climatic responses to solar forcing, rather than a strict northern-southern seesaw: different regions within the same hemisphere could have different significant levels of variations, even an opposite changing trend. Although some regions presented evident responses to solar forcing, other regions might have been insensitive to the driving forcing. This further indicated that if we considered the region-averaged temperatures on Earth (especially for regions in the southern hemisphere), it would be difficult to uncover the potential seesaw relation between them as the averaging process would weaken the differences between regions. This fact might explain why this millennial seesaw pattern of the Earth's climate was so highly controversial [47,48].

Conclusions and Discussion
In summary, we selected temperatures from 10 sites evenly distributed on the Earth's surface with time resolutions higher than 100 years to probe periodicities and cross correlations between solar variability and local temperatures on Earth. Spectral

Conclusions and Discussion
In summary, we selected temperatures from 10 sites evenly distributed on the Earth's surface with time resolutions higher than 100 years to probe periodicities and cross correlations between solar variability and local temperatures on Earth. Spectral analysis demonstrated that the millennial-scale cycle was the only common oscillation between SSN and all these local temperatures on Earth for our study interval (BC 6755-AD 1885). Further cross-wavelet analysis revealed that the millennial cycle components of SSN and temperatures had stable but varying values in the statistical correlations. These findings indicated that solar variability may have played a role in the modulation of this millennial variation of the Earth's climate, but the spatial patterns of climatic responses on the Earth are complex. For the 10 temperature records on the Earth, six gave positive correlations to SSN (five in the northern hemisphere, and one in the southern hemisphere), four gave negative correlations to SSN (all in the southern hemisphere) on the millennial time scale. Here tropical regions were classified into northern/southern hemisphere strictly. But the correlations were only statistically significant for five temperatures (four in the northern hemisphere with all positive correlations, and one in the southern hemisphere with negative correlation).
The mechanisms to explain the results obtained here involved two aspects: the solar impact mechanism and the response mechanism of the Earth's climate. For the former, the atmospheric general circulation models (GCMs) showed that when solar variability was low the changes in stratospheric ozone would trigger downward-propagating effects, which led to cooling in the high northern latitude atmosphere, a southward shift of the northern subtropical jet, and a decrease in the Northern Hadley circulation [4,49]. Bond et al. [50] pointed out that these atmospheric responses further led to an increase of the North Atlantic drift ice, cooling in Greenland, reduced monsoon activity, and reduced rainfall in low latitudes; then these atmospheric changes were amplified and transferred to other regions through their impacts on sea ice and North Atlantic thermohaline overturning [50][51][52]. Confirming the latter chain of mechanisms was far beyond the scope of this paper, but we still can list some useful clues here. The complicated oceanic circulation, especially that involving redistribution of heat through long-range advective transport, could lead to opposite trends of climate change in different regions on Earth [38,53]. We can speculate that solar variability first generates an overall positive feedback, like the atmospheric water vapor feedback [8], on the global Earth. Then this global change will lead to a negative feedback in local regions on Earth through complex oceanic and/or atmospheric processes. We note that results derived here are tentative and preliminary, and more observational evidences, rather than theoretical speculations, are needed in the next steps to generalize the findings, in addition to more credible physical mechanisms accounting for the complex responses of the Earth's climate.