Evidence for Solar Modulation on the Millennial-Scale Climate Change of Earth

In this study, we use available reconstructed data to investigate periodicities of solar activity (i.e., sunspot number) and the Earth’s climate change (temperatures of Lake Qinghai in China and Vostok in Antarctica, the GISP δ18O climate record of Greenland, and the stalagmite δ18O monsoon records of Dongge Cave in China) as well as their cross-wavelet coherences on millennial scale. We find that the variations of the Earth’s climate indices exhibited the 1000-year cyclicity, which was recently discovered in solar activity (called Eddy cycle). The cross-wavelet correlations between the millennium-cycle components of sunspot number and the Earth’s climate change remains both strong and stable during the past 8640 years (BC 6755–AD 1885). The millennial variation of sunspot number keeps in-phase with variations of Lake Qinghai temperature, Greenland temperature, and East Asian Monsoon, but anti-phase with the variation of Antarctica temperature. The strong and stable resonant relationships between sunspot numbers and these climate indices indicate that solar variability may have played a role in modulation on this millennial seesaw pattern of the Earth’s climate change before the modern industrial era.


Introduction
As the unique star in the solar system, the Sun is the main source of external energy for the Earth's climate system. Therefore, solar variability, through the intrinsic magnetic activity on the Sun, is believed to be one of the important natural driving factors for the Earth's climate change, e.g., [1][2][3]. The variations in solar activity are reflected by sunspot number (SSN), which is the solar parameter with the longest record time. Similarly, temperature is the most basic index to reflect climate change. Therefore, studies on the Sun-climate relationship often pay attention to the relations between SSN and temperatures on Earth. For example, Eddy [4] noticed that strong solar activity with numerous sunspots coincided with warm periods on Earth and weak solar activity with few sunspots coincided with cold periods. Reid [5] discovered common features and similarities between the secular variation of the globally averaged sea surface temperature, and the envelope of the 11-year running mean of SSN during the past 120 years. Friis-Christensen & Lassen [6] found the connection between the northern hemisphere temperature and the solar cycle length computed from SSN.
The role of solar variability on the Earth's climate change is one of the unresolved questions in the area of solar-terrestrial relations due to the complexity of the Earth's climate system [7]. In addition, the directly measured data for these parameters only cover timescales of centuries. Specifically, the directly observed SSN data cover only 400 years up until now (see http://www.sidc.be/silso/). The continuous measurement of the globally-averaged temperature of Earth started from AD 1850-1880, and the corresponding data covers only about 140-170 years until today. These are evidently insufficient for studying the long-term change of the Earth's climate over a millennium or even longer intervals.
However, the reconstruction methods can provide data over millennial timescales in both solar physics and paleoclimate studies. The SSN data over millennia could be computed indirectly from cosmogenic radionuclide proxy records ( 10 Be and 14 C) in terrestrial archives [8] because the contents of these cosmogenic elements are determined by the shielding effects of the Sun's charged particles output (i.e., solar wind), and the latter is related to the solar activity level (i.e., SSN). In paleoclimate, there are many ways to reconstruct the ancient climatic signals from various media, such as stable-isotope records from ice cores and stalagmites, alkenones in lake sediments, and so on.
In this paper, we will use the reconstructed data to explore the relations between the long-term changes of solar activity and the Earth's climate over a millennial timescale. Meanwhile, we aim to study the natural variations of the climate system before any potential or significant interferences from human factors. Therefore, we pay attention to the historical period before the modern industrial revolution.

Data and Methods
Wu et al. [8] obtained a new multi-proxy reconstruction of SSN, covering 8640 years from BC 6755 to AD 1885, by combining all available long-span datasets of 10 Be and 14 C in terrestrial archives. This is the most comprehensive and reasonable long-time reconstructed SSN data and can be download from http://www2.mps.mpg.de/projects/sun-climate/data/SN_composite.txt. Petit et al. [9] provided a data series of the local temperature at Vostok in East Antarctica (abbreviated as T A ; elevation of 3489 m) during the past 420,000 years based on the deuterium content of the ice (δD ice ) in Vostok's ice cores. We download this data from ftp://ftp.ncdc.noaa.gov/pub/data/paleo/icecore/antarctica/ vostok/deutnat.txt. Hou et al. [10] presented a reconstructed summer water temperature of Lake Qinghai (on the northeastern Tibetan Plateau, abbreviated as T Q ; elevation of 3260 m) in China. This reconstruction was based on alkenones in the lake's sediments. This data was downloaded from NOAA paleoclimate website (https://www.ncdc.noaa.gov/data-access/paleoclimatology-data). Stable-isotope records from ice cores and stalagmites are also widely used to reveal past climatic changes. The Greenland Ice Sheet Project 2 (GISP2) provided the δ 18 O climate record, the index of temperature, of the past 16,500 years from ice cores (abbreviated as δ 18 O G ) [11] that can be obtained from http://depts.washington.edu/qil/datasets/gisp2_main.html. The stalagmite δ 18 O records of Dongge Cave reconstructed the history of the East Asian summer monsoon activity over the past 9000 years (abbreviated as δ 18 O D ) [12] that can be downloaded from NOAA paleoclimate website (https://www.ncdc.noaa.gov/data-access/paleoclimatology-data). Due to different source data and reconstruction methods, the ranges of data are also different. The SSN series is given in year, while the other four climate indices (T A , T Q , δ 18 O G , δ 18 O G ) are given in age; here all ages are in years BP (before present, i.e., AD 1950). In this study, we use the data under the uniform time resolution of 10 years, and select the common section from BC 6755 to AD 1885 as the time period under study.
The Lomb-Scargle periodogram method [13,14] is adopted to analyze the periodicities of both SSN and these climate indices of Earth. This method extracts useful signals from potentially noisy backgrounds and can be applied to incomplete or unevenly sampled time series. The wavelet coherence (WTC) tool [15,16] is used in this study to examine the cross-correlations and phase angles between solar activity and the Earth's climate change in the 2D space of period/frequency and time. For the two series to be physically related, we could expect both a high-power coherence and a consistent phase relation between them in the cross-wavelet coherence analysis [16]. This is a significant advantage over the general Pearson correlation analysis. The software package was downloaded from http://noc.ac.uk/marine-data-products/cross-wavelet-wavelet-coherence-toolbox-matlab.
Universe 2020, 6, 153 3 of 9 Figure 1 shows the Lomb-Scargle periodograms of these time series. It can be seen that the variations of SSN have some cyclicities higher than the 95% confidence level (red horizontal dotted line), i.e., 150 year, 208 year (Suess or de Vries cycle [17]), 350 year, 500 year, 720 year, 990 year (Eddy cycle [4,[18][19][20]), and so on. This is consistent with previous studies [21]. As a contrast, the variations of these climatic indices are smooth with less cyclicities. For T Q , only a cyclicity at 962 year is significant during our study interval; For T A , there are some significant cyclicities, such as 188 year, 295 year, 483 year, 730 year, 1020 year; For δ 18 O G , the significant cyclicities include 532 year, 585 year, 719 year, 885 year, 990 year, and 1149 year; While for δ 18 O D , no cyclicity is found significant (all below the 95% confidence level). The only common period between SSN and these climate indices (temperatures and δ 18 O) is the millennial-scale cycle (known as the Eddy cycle), which is 990-year period for SSN and δ 18 O G , 962-year period for T Q and δ 18 O D , 1020-year period for T A . All of them are over 95% confidence level except for δ 18 O D . The common periodicity between them suggests that there might be some physical relationship between solar activity and the Earth's climate change on this millennial scale.

Eddy Cycle in SSN and Climate Indices
Universe 2020, 6, x FOR PEER REVIEW 3 of 10 over the general Pearson correlation analysis. The software package was downloaded from http://noc.ac.uk/marine-data-products/cross-wavelet-wavelet-coherence-toolbox-matlab. Figure 1 shows the Lomb-Scargle periodograms of these time series. It can be seen that the variations of SSN have some cyclicities higher than the 95% confidence level (red horizontal dotted line), i.e.,: 150 year, 208 year (Suess or de Vries cycle [17]), 350 year, 500 year, 720 year, 990 year (Eddy cycle [4,[18][19][20]), and so on. This is consistent with previous studies [21]. As a contrast, the variations of these climatic indices are smooth with less cyclicities. For TQ, only a cyclicity at 962 year is significant during our study interval; For TA, there are some significant cyclicities, such as 188 year, 295 year, 483 year, 730 year, 1020 year; For δ 18 OG, the significant cyclicities include 532 year, 585 year, 719 year, 885 year, 990 year, and 1149 year; While for δ 18 OD, no cyclicity is found significant (all below the 95% confidence level). The only common period between SSN and these climate indices (temperatures and δ 18 O) is the millennial-scale cycle (known as the Eddy cycle), which is 990-year period for SSN and δ 18 OG, 962-year period for TQ and δ 18 OD, 1020-year period for TA. All of them are over 95% confidence level except for δ 18 OD. The common periodicity between them suggests that there might be some physical relationship between solar activity and the Earth's climate change on this millennial scale.

Cross-Wavelet Analysis
We use the WTC tool to compute the wavelet coherence between SSN and these climate indices. Figure 2 gives the squared wavelet coherence (C w ) between SSN and climate indices. Here, the squared wavelet coherence (C w = C w (t,p)) is defined in the time (t)-period (p) space. This parameter discloses correlations of two signals in the t-p space. We adopt the global wavelet coherence Cg( ) and angle strength Sθ( ) at a given to measure the cross-correlations between SSN and each climate index quantitatively. Cg( ) is defined as t1, t2, …tn are discrete moments in the study interval, and Cg( ) depicts the averaged correlation between SSN and climate signal at a given p. However, only high power of coherence is not sufficient to reveal the causal relationship between two variables. Their phase angles should be stable in addition to a high correlation if two variables are physically related. The phase angle strength Sθ( ) is defined as where θi( ) (I = 1, 2, 3,…, n) are the phase angles in the cross-wavelet transform at time t1, t2, t3…tn. Sθ = 1 implies strictly stable phase angles, and Sθ=0 implies non-stable phase angles. Details can be found in [16]. We compute Cg( ) and ( ) satisfying the COI ('Cone Of Influence') condition, and plot them in  The red regions enclosed by thick black lines denote their confidence levels higher than 95% against red noise. The light black curve denotes the 'cone of influence' (COI) outside which the edge effects might distort the picture. Black arrows represent phase angles of two records in the cross-wavelet analyses.
We can see from Figure 2 that the regions with strong cross-wavelet correlations (bright-color regions) usually appear intermittently, indicating that correlations between SSN and climate indices are not stable: they remain strong at a certain time, then become weak after that. The exception is for the millennial cyclicity. A high-power region at 1000 year persists strong nearly throughout the whole study interval for each climate index, shown as yellow to red regions of high power in each panel of this figure (but with relatively weak power for T A in (b)). This is the most evident for T Q . It means that SSN has some kinds of correlations to climate indices (temperature and δ 18 O) on this millennial cycle. In other words, the millennial components of SSN and climate indices are correlated for the whole study interval.
We adopt the global wavelet coherence C g (p) and angle strength S θ (p) at a given p to measure the cross-correlations between SSN and each climate index quantitatively. C g (p) is defined as (1) t 1 , t 2 , . . . , t n are discrete moments in the study interval, and C g (p) depicts the averaged correlation between SSN and climate signal at a given p. However, only high power of coherence is not sufficient to reveal the causal relationship between two variables. Their phase angles should be stable in addition to a high correlation if two variables are physically related. The phase angle strength S θ (p) is defined as where θ i (p) (I = 1, 2, 3, . . . , n) are the phase angles in the cross-wavelet transform at time t 1 , t 2 , t 3 , . . . , t n . S θ = 1 implies strictly stable phase angles, and S θ = 0 implies non-stable phase angles. Details can be found in [16]. We compute C g (p) and S θ (p) satisfying the COI ('Cone Of Influence') condition, and plot them in Figure 3. Both C g and S θ have a peak simultaneously at p = 1000 year (red dotted line) in cases of all four climate indices. For SSN and T Q , their wavelet coherence is strong (C g = 0.64) and stable (S θ = 0.96).
While for SSN and T A , their wavelet coherence is stable (S θ = 0.56), but not so strong (C g = 0.25) as that of SSN and T Q . For SSN and δ 18 O G , S θ = 0.81, C g = 0.37, at p = 990. 6 year. For SSN and δ 18 O D , S θ = 0.83, C g = 0.41 at this periodicity. These demonstrate that the wavelet coherences between SSN and these climate indices are evidently stable and generally strong for the 1000-year periodicity.

Millennial Oscillations and Solar Modulation
In order to show the millennial variations of both solar activity and climate change on Earth, we use a band-pass filter to extract millennial signals from data. The filter is centered at 1000 years with half-power points at 1000 ± 100 year. Figure 4 shows the output signals after the filter (thin lines). The millennial variations of SSN have a consistent trend with those of TQ and δ 18 OG, and positive correlations lie between them. As for TA and δ 18 OD, their changing trends are opposite to that of SSN with negative correlations. Figure 4 also provides a direct comparison between SSN and climate indices during the past 8640 years. As the main purpose of this paper is to probe long-term changes of solar/climate change, we show the 200-year running averages of the original data (thick lines). The corresponding millennial components are denoted as thin lines in this figure. Compared with solar variability, the variations of Earth's climate indices (TQ, TA, δ 18 OG, δ 18 OD) are smooth. Green vertical lines denote peaks (valleys) of SSN, TQ, and δ 18 OG that correspond to valleys (peaks) of TA, and δ 18 OD. The general trends of SSN, TQ, δ 18 OG are usually opposite to those of TA, δ 18 OD between vertical lines.
We compute the linear Pearson correlation coefficient (CC) between SSN and climate indices for their millennial components and show them in Table 1. The CCs for TQ and δ 18 OG are 0.76 and 0.58, demonstrating strong positive correlations between them. The CCs for TA and δ 18 OD are −0.74 and −0.67, implying clearly negative correlations with SSN. We also show the significance of the CC in Table 1, which is computed between two time series after taking into account the effect of their autocorrelations [22,23]. However, we should be careful in interpreting the statistical results. Considering the fact that δ 18 OG is more indicative of the local temperature [11], therefore the positive correlation between SSN and δ 18 OG indicates a positive correlation between SSN and Greenland temperature. While a smaller δ 18 OD value corresponds to a stronger East Asian Monsoon intensity [12], so the negative correlation between SSN and δ 18 OD indicates a positive correlation between SSN and Monsoon intensity in East Asia. In a summary, the millennial variation of SSN should be in-phase with variations of Lake Qinghai temperature, Greenland temperature, and East Asian Monsoon intensity, but anti-phase with the variation of Antarctica temperature. The strong and stable resonant

Millennial Oscillations and Solar Modulation
In order to show the millennial variations of both solar activity and climate change on Earth, we use a band-pass filter to extract millennial signals from data. The filter is centered at 1000 years with half-power points at 1000 ± 100 year. Figure 4 shows the output signals after the filter (thin lines). The millennial variations of SSN have a consistent trend with those of T Q and δ 18 O G , and positive correlations lie between them. As for T A and δ 18 O D , their changing trends are opposite to that of SSN with negative correlations.  We can present a coarse quantitative estimation for the Sun's modulation on the Earth's climate change on the millennial scale. Specifically, the 1000-year cycle of SSN varies within ±3.75; while the original SSN data (10-year resolution) changes from −8.0 to 77.9; thus, the millennial component accounts for 9% of its total variations for SSN. As to these four climate indices, the millennial variation of TQ is ±0.56 °C, accounting for 9% of its total variability; the millennial variation of TA is ±0.27 °C, accounting for 14% of its total variability; the millennial variation of δ 18 OG is ±0.16 (‰), accounting for 12% of its total variability; and the millennial variation of δ 18 OD is ±0.06 (‰), accounting for 5% of its total variability. As a summary, we can say that 9% of total solar variability had modulated 9%-14% of the Earth's temperature change and 5% of the East Asian Monsoon activity, on the millennial scale during our study interval as a rough estimate.

Summary and Discussion
In this paper, we analyzed the reconstructed datasets to probe periodicities and cross-correlations between solar variability (SSN) and climate indices on Earth (summer temperature at Lake Qinghai (TQ), air temperature at Vostok in Antarctica (TA), temperature index at Greenland (δ 18 OG), and Monsoon index at Dongge Cave in China (δ 18 OD)) during the past 8640 years (BC 6755-AD 1885) before the modern industrial era. We find that the millennial-scale cycle is the only common oscillation between SSN and these climate indices on Earth during our study interval. The millennial-cycle components of SSN and climate indices have both strong and stable correlations. When solar activity is strong (numerous SSN), it will give warming effect on Greenland and Lake Qinghai, a strengthening effect on the monsoon in East Asia, but a cooling effect on Vostok in Antarctica.  Table 1, which is computed between two time series after taking into account the effect of their autocorrelations [22,23]. However, we should be careful in interpreting the statistical results. Considering the fact that δ 18 O G is more indicative of the local temperature [11], therefore the positive correlation between SSN and δ 18 O G indicates a positive correlation between SSN and Greenland temperature. While a smaller δ 18 O D value corresponds to a stronger East Asian Monsoon intensity [12], so the negative correlation between SSN and δ 18 O D indicates a positive correlation between SSN and Monsoon intensity in East Asia. In a summary, the millennial variation of SSN should be in-phase with variations of Lake Qinghai temperature, Greenland temperature, and East Asian Monsoon intensity, but anti-phase with the variation of Antarctica temperature. The strong and stable resonant relationships between SSN and climate indices indicate that solar variability might have played a role in modulation on this millennial seesaw pattern of the Earth's climate change. We can present a coarse quantitative estimation for the Sun's modulation on the Earth's climate change on the millennial scale. Specifically, the 1000-year cycle of SSN varies within ±3.75; while the original SSN data (10-year resolution) changes from −8.0 to 77.9; thus, the millennial component accounts for 9% of its total variations for SSN. As to these four climate indices, the millennial variation of T Q is ±0.56 • C, accounting for 9% of its total variability; the millennial variation of T A is ±0.27 • C, accounting for 14% of its total variability; the millennial variation of δ 18 O G is ±0.16 (% ), accounting for 12% of its total variability; and the millennial variation of δ 18 O D is ±0.06 (% ), accounting for 5% of its total variability. As a summary, we can say that 9% of total solar variability had modulated 9-14% of the Earth's temperature change and 5% of the East Asian Monsoon activity, on the millennial scale during our study interval as a rough estimate.

Summary and Discussion
In this paper, we analyzed the reconstructed datasets to probe periodicities and cross-correlations between solar variability (SSN) and climate indices on Earth (summer temperature at Lake Qinghai (T Q ), air temperature at Vostok in Antarctica (T A ), temperature index at Greenland (δ 18 O G ), and Monsoon index at Dongge Cave in China (δ 18 O D )) during the past 8640 years (BC 6755-AD 1885) before the modern industrial era. We find that the millennial-scale cycle is the only common oscillation between SSN and these climate indices on Earth during our study interval. The millennial-cycle components of SSN and climate indices have both strong and stable correlations. When solar activity is strong (numerous SSN), it will give warming effect on Greenland and Lake Qinghai, a strengthening effect on the monsoon in East Asia, but a cooling effect on Vostok in Antarctica. Conversely, when solar activity is weak (few SSN), it will create a cooling effect on Greenland and Lake Qinghai, a weakening effect on the monsoon in East Asia, but a warming effect on Vostok in Antarctica. These millennial variations account for roughly 10% of the total changes for both solar variability and climate change.
We also find that the climatic response to solar forcing in Vostok of Antarctica is opposite to those in the northern hemisphere (Lake Qinghai, Dongge Cave, Greenland). Similarly, Eroglu et al. [24] discovered a seesaw relationship of the Holocene East Asian-Australian summer monsoon over the last 9000 years at millennial to sub-centennial time-scales, and concluded that solar activity acts as a driver in the seesaw relationship between the East Asian summer monsoon and the Indonesian-Australian summer monsoon by shifting the position of the Intertropical Convergence Zone. The opposite changing trend in the millennial temperatures of two hemispheres had been discovered for 20 years [25][26][27][28][29][30][31], which can be termed the "northern-southern hemispheric seesaw". More and more evidence indicates that this dynamic see-saw pattern seems to be a persistent feature of long-term climate change [27,28].
The complicated oceanic and/or atmospheric processes are thought to be possible explanations for this seesaw phenomenon, and the driving center of action for the millennial-scale seesaw could be initiated from the northern hemisphere [27,32]. In the case with emphasis on the northern hemisphere playing a leading role for millennial-scale oscillation or see-saw, an increase in the North Atlantic thermohaline circulation would warm the high-latitude northern hemisphere, but cool the southern hemisphere on a millennial scale [26,27]. The modeling based on the climate model demonstrated that the millennial variations of the northern temperature keep consistent with the global change, and only a small part of the very southern regions presents a converse trend with the global variation [29].
Therefore, we can give the following scenario: solar activity first generates positive feedback on the global Earth. For example, solar irradiance as well as geomagnetic storms produced by solar activities can have effects on the Earth's atmosphere [33][34][35]. The positive feedback in the ocean-atmosphere system can also amplify the response to solar irradiance variations [36,37]. Then this global change will lead to negative feedback in the local Antarctic through complex oceanic and/or atmospheric processes. However, it is important to note that our study here is only based on data from four places on Earth. For future steps, we need to find more evidences from different places all over the world to verify this 'north-south' seesaw pattern as well as more credible physical mechanisms accounting for the seesaw oscillations between climate indices in the northern and southern hemisphere. The next research goal must also be able to withstand worthy criticisms on the north-south seesaw mechanism originally leveled by Wunsch [38].