Possible link between Holocene East Asian monsoon and solar activity obtained from the EMD method

It is thought that East Asian monsoon (EAM) is linked and sensitive to solar activity. In this paper, we have decomposed the Dongge cave speleothem δ18O record (proxy for EAM), and114C and10Be (proxies for solar activity) time series into variations at different time scales with the empirical mode decomposition (EMD) method to reveal the possible link between the EAM variability and solar activity. There are some common cycles in the EAM and solar variability from centennial to millennial scales, indicating a possible link between EAM and solar activity at these time scales. The correlation between EAM and solar activity is much higher at millennial scales than at centennial scales, which means direct responses to the solar variation are more likely at time scales longer than a few hundred years. At ∼ 30, 60 and 600 yr time scales, the variation in EAM is amplified by the solar amplitude modulation at ∼ 100, 200 and 2200 yr time scales.


Introduction
The Holocene has long been regarded as an analog to the future of natural environment process, which began at the end of the Pleistocene (around 12 000 yr BP) and continues to the present.Detailed studies of Holocene climate change and the controlling mechanisms can provide a context for present and future global change, and increase our understanding of the underlying mechanisms of natural climate variability (Xue and Zhong, 2011).The East Asian monsoon (EAM) is an integral part of the global climatic system, which regulates global atmospheric circulation through heat and moisture transport from the warmest part of the tropical ocean, the west Pacific Warm Pool (WPWP), to higher latitudes (Wang et al., 2005), and plays a significant role in the socio-economic life of people in East Asia (Zong et al., 2006).Increasing evidence shows that there were abrupt changes and periodic fluctuations in East Asian monsoon on decadal to millennial scales during the Holocene (Sirocko et al., 1993;Hong et al., 2000;Gupta et al., 2003;Fleitmann et al., 2003).There are many factors affecting the East Asian monsoon, including orbitally induced insolation changes, changes in solar output, and changes in oceanic and atmospheric circulation (Wang et al., 2005).A number of existing high-resolution Holocene monsoon reconstructions from Indian and Asian monsoon domains, and some modelling studies, reveal an influence of solar variation on East Asian monsoon (Joussaume et al., 1999;Neff et al., 2001;Fleitmann et al., 2003;Gupta et al., 2003;Wang et al., 2005), which can be explained by a direct solar influence on the Intertropical Convergence Zone (ITCZ) that controls the monsoonal precipitation (Kodera, 2004).Asian monsoon is thought to be sensitive to relatively small changes in solar forcing (0.25 % change in solar output) (Neff et al., 2001;Fleitmann et al., 2003).The mechanism is not very clear (Cosford et al., 2008).It is suggested to be amplified by the Pacific Decadal Oscillation and El Niño-Southern Oscillation (ENSO) (Mann et al., 2005), and the North Atlantic Oscillation (Gimeno et al., 2003) or be due to the threshold effects (Dykoski et al., 2005).
Most of the previous studies explored the possible link between Asian monsoon and solar activity by direct comparisons among different proxy records (Hu et al., 2003;Gupta et al., 2003) or by the presence of some periodicities and coherence at certain bandwidths with respect to the 14 C and 10 Be data with spectral analysis (Dykoski et al., 2005;Xiao et al., 2006) and wavelet analysis (Witt and Schumann, 2005;Rousse et al., 2006;Cosford et al., 2008) (Brigham, 1988;Qian, 2002).Short-time Fourier transform is able to capture the time-dependence of frequency fluctuations from nonstationary data.However, in order to localize an event precisely in time, the window width must be narrow.Alternately, the frequency resolution requires longer time spans.This leads to conflicting requirements (Heisenberg-Gabor inequality) and restrains this method from many practical applications (Barnhart and Eichinger, 2011).The wavelet analysis has been applied to diagnose climate changes.It has good ability to make multi-resolution analysis in both time domain and frequency domain.However, the choice of wavelet basis functions limits the applicability of the technique, as the basic functions of wavelet transformation are fixed and do not necessarily match the shape of the considered data series in every instant in time (Torrence and Compo, 1998).
A relative new time-frequency analysis method, designated empirical mode decomposition (EMD), has been proposed by Huang et al. (1998).The method essentially decomposes the original signal into a number of intrinsic mode function (IMF) components based on the local characteristic time scale of the signal and a residue.The EMD method is self-adaptive because the IMFs, working as the basis functions, are determined by the signal itself rather than what is pre-determined.Therefore, the EMD method is efficient in data analysis (Huang et al., 1998;Huang and Wu, 2007;Wu et al., 2007).In this paper, we will apply the EMD method to reveal the possible link between EAM and solar activity.

EMD method
According to the definition of the EMD method (Huang et al., 1998;Huang and Wu, 2007;Wu et al., 2007;Huang and Shen, 2005), the signal is decomposed into a finite number of IMFs, each of which represents a simple oscillatory mode and satisfies the following two conditions to ensure there are no riding waves (i.e. the wave motion formed by small-scale perturbation waves superposed on long-scale unperturbed waves): (1) in the whole dataset of the IMF, the number of extrema and the number of zero crossings must either be equal or differ at most by one; (2) at any point of the IMF, the mean value of the envelopes defined by the local maxima (relative maximum values that are greater than those values at the surrounding points), and local minima (relative minimum values that are smaller than those values at the surrounding points) is zero.
The decomposition method uses envelopes defined by the local maxima and the local minima of the original time series X(t), respectively.Firstly, all local extrema are identified, and then all local maxima are connected by a cubic spline to form the upper envelope and minima are connected to form the lower envelope.Their mean is designated as m 11 (t), and the difference between X(t) and m 11 (t) is the first component, h 11 (t), i.e. h 11 (t) = X(t) − m 11 (t) . (1) However, h 11 (t) is still not stationary.So, the above process is repeated.
Further iterations of this sifting process proceed until a standard deviation (SD) criterion is met.The SD is typically set between 0.2-0.3(Huang et al., 1998).
T is the time length.When SD criterion is satisfied after k iterations, then The first IMF component (c 1 = h 1k (t)) is separated from the rest of the data by Since r 1 still contains information of longer period components, it is treated as the new data and subjected to the same sifting process as described above.The process is repeated and r i is described as follows: The sifting process can be stopped by any of the following predetermined criteria: either when the component, c n , or the residue, r n , becomes so small that it is less than the predetermined value of substantial consequence, or when the residue, r n , becomes a monotonic function from which no more IMF can be extracted.Thus, we decompose X(t) into c i = IMF i (t) with i = 1, 2, 3 • • • , n and a residue, r n , which can be either the mean trend or a constant.
where ξ = 0.01 and the standard deviation of the mean-zero white noise is 10 % of the signal strength.
The signal was decomposed into four IMF components and a residue with the EMD method (Fig. 1).The frequencies of the four IMF components are 1/21, 1/40.8,1/99.7, and 1/280 Hz from high to low.The frequencies of IMF1, IMF2 and IMF3 are in good accord with the original three frequencies.The IMF4 and the residual may be mainly contributed by the Gaussian white noise.So, EMD method can decompose the original signal into a number of IMF components from high frequency to low frequency in order and adaptively.
Although the EMD method is powerful, one difficulty encountered when using EMD is the influence of the end point treatment.The envelopes are calculated using a cubic spline.However, splines are notoriously sensitive to end points.It is important to make sure that the end effects do not propagate into the interior solution (Lin and Wang, 2006).

Data source
Speleothems can have continuous deposition of calcium carbonate over long periods of time, and well-chosen speleothems are datable with high precision.Absolute ages can be determined by means of 230 Th dating by mass spectrometry (Dykoski et al., 2005;Cosford et al., 2008).Previous studies have shown that shifts in the oxygen isotope ratio (δ 18 O) of the stalagmite from the cave largely reflect changes in δ 18 O values of meteoric precipitation at the site,  which in turn relate to changes in the amount of precipitation and thus characterize the Asian monsoon strength (Dykoski et al., 2005).δ 18 O time series (from 9 ka BP to the present; the present is defined as the year 1950 AD) of stalagmite DA from Dongge cave, southern China, with an average temporal resolution of 4-5 yr are used as a proxy of the high-resolution absolute-dated Holocene Asian monsoon record (Wang et al., 2005).Chronology of the 962.5-mm-long stalagmite DA is established by 45 230 Th dates, all in stratigraphic order, with a typical age uncertainty of 50 yr (Wang et al., 2005).
Owing to the lack of complete direct observations of the Sun for the period before 1600 AD, cosmogenic radionuclides 14 C (recorded in tree rings) and 10 Be (preserved in polar ice cores) are used as indirect proxies of solar variability (Yiou et al., 1997;Stuiver et al., 1998;Wagner et al., 2001;Muscheler et al., 2007).During the periods of lower sunspot activity, solar wind intensity is reduced, which increases the influx of galactic cosmic rays.A higher influx of cosmic rays increases the production of 14 C and 10 Be in the atmosphere.The reverse is true during the periods of increased sunspot activity when less 14 C and 10 Be are produced.Therefore, changes in atmospheric 14 C and 10 Be can be related to changes in solar activity.In this paper, the tree-ring 14 C record (from 9000 yr BP to the present) is from the IntCal04 calibration curve with 1-standard deviation envelope and data with 1-standard deviation error bars in 14 C and calibrated age (Reimer et al., 2004) .The 10 Be record (from 304 to 9315 yr BP) is from Greenland Ice Core Project (GRIP) ice core (Yiou et al., 1997;Wagner et al., 2001;Muscheler et al., 2007) with 5 yr temporal resolutions, and we use the timescale published by Johnsen et al. (1997) (ss09 timescale) by counting annual layers in the GRIP core.
The datasets of δ 18 O, 14 C and 10 Be spaced at 10 yr intervals (Fig. 2) are decomposed into signals at different time scales for further comparisons with the EMD method to reveal the possible link between solar activity and EAM intensity on different time scales.

EMD analysis of δ 18 O, 14 C and 10 Be
The EMD method was utilized to extract the intrinsic cycles of the datasets of δ 18 O, 14 C and 10 Be (Fig. 3).The time series of δ 18 O was decomposed into seven IMF components and a residue, and the time series of 10 Be and 14 C were decomposed into six IMF components and a residue respectively.The cycle of each IMF component is calculated using zero-crossing method.
The cycles of the IMF components of 14 C centre on 43, 102, 215, 389, 750, 2291 yr.The cycles of 102, 215, 389 and 2291 yr are close to 104, 207, 350 and 2400 yr, which are strong in spectral analysis (Fig. 4b).These cycles are also revealed by previous studies (Stuiver and Braziunas, 1993;Ogurtsov et al., 2002).Although ∼ 750 and 2291 yr are not observed in the spectral analysis, ∼ 750 yr cycle is reported in another study (Clemens, 2005), and ∼ 2291 yr cycle is close to the Hallstatt cycle, which is revealed by Miller and Scott (1995) and Usoskin et al. (2009).
The cycles of the IMF components of 10 Be centre on 49, 103, 220, 480, 1005 and 2200 yr, which is similar to the results of Yin et al. (2007).Cycles of 49, 103, 220 and 1005 yr are also revealed in spectral analysis (Fig. 4c).
By the comparison of solar activity reconstructions from 10 Be GISP series and from 14 C series, Vonmoos et al. (2006) have shown that the agreement between them is good on centennial to millennial timescale.With the EMD method, the records of 14 C and 10 Be also show concurrent cycles on ∼ 100, 200, 1000 and 2200 yr.The correlation coefficients on these scales are higher than on multi-decadal scales (Table 2).Although the cycle of IMF1 of 14 C (44 yr) is close to that of 10 Be (49 yr), the correlation coefficient is very low (Table 2).Herein, good agreement between 14 C and 10 Be on centennial to millennial timescale implies that the variations at these cycles are more likely caused by the solar modulation than at multi-decadal scales.
Although 14 C and 10 Be are both used as proxy of solar activity, there are some discrepancies in their cycles.The cycle of IMF4 of 14 C is 389 yr, while that of 10 Be is 480 yr.The 389 yr cycle is close to a 360 yr tidal cycle (Keeling and Whorf, 2000), and observed in a high-resolution δ 18 O record of peat cellulose in the Jilin Province in China (Hong et al., 2000), and in slope sediments of the Great Bahama Bank (Roth and Reijmer, 2005).The worldwide appearance of this signal thus points to the involvement of oceanic and    (Roth and Reijmer, 2005).The ∼ 500 yr oscillation is thought to reflect a climatic process, involving thermohaline circulation and atmospheric processes that influenced atmospheric and oceanographic dynamics (Roth and Reijmer, 2005).Moreover, this cycle is possible to be a harmonic of the reported millennial-scale cycles or a subharmonic of the reported century-scale cycles.
Previous studies have pointed out that the short-term variations in 10 Be are mainly controlled by local climate (Muscheler et al., 2007;Usoskin et al., 2009).The EMD analysis shows that the amplitudes of 10 Be variation at short time scales are larger than those at long time scales (Fig. 3c).Moreover, the spectral analysis reveals there are many multidecadal cycles in 10 Be record (Fig. 4c).So, local climate may dominate the 10 Be data on short timescales as Usoskin et al. (2009) pointed out.Conversely, the amplitudes of shorttime variation in 14 C are smaller than those of long-term variation (Fig. 3b), which is suggested to be due to the amplification of carbon cycle at long time scales by Muscheler et al. (2007).

Relations between solar variability and East Asian monsoon at different time scales
Most of the cycles found in Sect. 3 have been reported from mid to high northern latitudes (Bond et al., 2001;Hu et al., 2003) to low-latitude regimes (Fleitmann et al., 2003), including the Asian monsoon region (Wang et al., 2005;Xiao et al., 2006;Thamban et al., 2007).The cycles of ∼ 118, 234, 1050 and 2000 yr closely match those of 14 C and 10 Be. ∼ 118 and 234 yr cycles are close to 90 yr (Gleissberg) and 207 yr (Suess) solar cycles.The century-type cycle of Gleissberg has a wide frequency band with a double structure consisting of 50-80 yr and 90-140 yr periodicities, and the Suess cycle is less complex showing a variation with a period of 170-260 yr (Ogurtsov et al., 2002).These two cycles have been reported from East Asian monsoon (Dykoski et al., 2005;Xiao et al.,2006;Cosford et al., 2008), Indian monsoon (Neff et al., 2001;Agnihotri et al., 2002;Fleitmann et al., 2003;Thamban et al., 2007) and American monsoon (Asmerom et al., 2007;Stríkis et al., 2011) based on a variety of climate proxies.The ∼ 1000 yr cycle is widely reported in stalagmites from Oman (Neff et al., 2001), sediments of the Arabian Sea (Thamban et al., 2007), peat cellulose from Hongyuan and Jinchuan (Xu et al., 2002) and lake sediments from Alaska (Hu et al., 2003).Its global nature due to the hemisphere-wide presence of this periodicity (Thamban et al., 2007) and its presence in 10 Be data support that it is related to the variations in solar activity.The ∼ 2000 yr cycle was recorded in sediments of the Arabian Sea (Thamban et al., 2007) and Yangtze River-derived mud in the East China Sea (Xiao et al., 2006).The cycle was close to the 2560 yr recorded in the Camp Century ice core δ 18 O profile (Dansgaard et al., 1984) and in the Okinawa Trough (Jian et al., 2000).The 14 C and 10 Be records both reveal ∼ 2200 yr cycle, which is also revealed by other studies (Stuiver and Braziunas, 1993;Lean, 2002) and is close to Hallstattzeit cycle.Hence, ∼ 2200 yr cycle may be linked to solar variability.Moreover, oceanic circulation is thought to play a significant role in transferring and modulating the variability at this time scale (Naidu and Malmgren, 1995;Thamban et al., 2007).
Although common cycles on ∼ 100, 200, 1000 and 2000 yr do not necessarily suggest a mechanistic link, they do indicate a possible link between the EAM and solar activity.As shown from Table 3, at millennial scales (∼ 1000 and 2000 yr), the correlations between EAM and solar activity are much higher than at centennial scales, which means solar activity has greater influences on EAM at millennial scales.
The ∼ 30 and 60 yr cycles are very common in not only the modern but also secular monsoon records, including East Asian monsoon and India monsoon (Minobe, 1997;Agnihotri et al., 2002;Clemens, 2005;Dykoski et al., 2005;Wang et al., 2005;Mazzarella and Scafetta, 2012).However, both 14 C and 10 Be records do not reveal these two cycles with the EMD method.It seems there is no direct relation between    solar variability and EAM variation at multi-decadal time scales.However, the atmospheric pressure changes associated with the North Atlantic Ocean (NAO), and changes in sea-surface temperatures and thermohaline circulation associated with the Atlantic Multidecadal Oscillation (AMO) exhibit a ∼ 60 yr cycle (Schlesinger and Ramakutty, 1994;Kerr, 2000;Gimeno et al., 2003;Knudsen et al., 2011).Moreover, the Pacific Decadal Oscillation (PDO) resembles the well-known ENSO system, but operates in the North Pacific and exhibits 25-35 yr and 50-70 yr cycles (Minobe, 1997).So, ∼ 30 and 60 yr cycles are thought to be linked to the Pacific Decadal Oscillation and the North Atlantic Oscillation (Schlesinger and Ramankuty, 1994;Gimeno et al., 2003;Knudsen et al., 2011;Mazzarella and Scafetta, 2012).
∼ 600 yr cycle is revealed in the Lianhua A1, Heshang HS4, Dongge DA and D4 speleothem δ 18 O records (Cosford et al., 2008).It is close to the prominent 649 yr cycle discovered by Damon and Sonett (1991) and is thought to be related to solar variability (Zhu et al., 2009).However, this cycle is not observed in 14 C and 10 Be records with the EMD method.Hence, it seems that solar variability has few or indirect influence on EAM at this time scale.
The residues of δ 18 O, 14 C and 10 Be have been compared to show their long-term trends (Fig. 5).The residue of δ 18 O increased gradually, which means the Asian monsoon weakened through the Holocene.It is thought to be caused by the gradual decrease in summer insolation (Overpeck et al., 1996).The long-term trend of δ 18 O is quite different from that of 14 C and 10 Be, which means that different driving forces control them.It is because the radionuclide production rates on long time scales (> 3000 yr) are probably caused by geomagnetic modulation and those on short time scales (< 3000 yr) are most likely of solar origin (Beer, 2000).There is a slight discrepancy between the long-term trend of 14 C and 10 Be.Whether this discrepancy is mainly due to changes in the 14 C or the 10 Be system cannot yet be decided without further information (Vonmoos et al., 2006;Muscheler et al., 2007).Both radionuclides exhibit the same solar signal despite the fact that different conditions influenced their concentrations before their deposition in the archives (Vonmoos et al., 2006).However, 10 Be has a better correlation with δ 18 O than 14 C (Table 3).Especially, at ∼ 200 yr time scale, correlation coefficient between 14 C and δ 18 O is very low (r = −0.0161),while that between 10 Be and δ 18 O is remarkable (r = 0.2385).It seems that there is almost no correlation between 14 C and δ 18 O at this time scale.Therefore, 10 Be is used as a proxy of solar activity for further comparison with δ 18 O in this article.
∼ 30, 60 and 600 yr cycles are absent in 10 Be and 14 C records with the EMD method.Although these cycles are observed in the spectral analysis in other studies (Neff et al., 2001;Fleitmann et al., 2003), they are too small to cause significant impact on the climate of the Earth.So, these three time scales are thought to be indirectly affected by solar variability through the interaction between atmospheric and ocean circulation that amplify this initial signal (Neff et al., 2001).With the EMD method, the amplitudes of these three IMFs of EAM change with time and, in fact, have characteristics of slowly varying wavelet packets (Fig. 4).So, there are possible amplitude modulations of EAM by some driving forces at these time scales.The envelope of an IMF is symmetric, so it is only necessary to examine the maximum values in order to identify the characteristics of this cycle (Lin and Wang, 2003).To reveal possible solar modulation of EAM, the envelopes of IMF1, IMF2 and IMF5 of δ 18 O have been compared with IMF2, IMF3 and IMF6 of 10 Be respectively (Fig. 6).The crests of the envelope of the IMF1 of δ 18 O largely coincide with the troughs of the IMF2 of 10 Be and vice versa (Fig. 6a), which means there is an inverse phase relation between them.Similarly, the envelopes of IMF2 and IMF5 of δ 18 O are also in inverse phase with the IMF3 and IMF6 of 10 Be in great part.So, ∼ 100, 200, 2000 yr solar variability are adaptations to ∼ 30, 60 and 600 yr EAM variation, and EAM variation is amplitude modulated by solar variation.
The ∼ 200 yr Suess cycle and 2200 yr Hallstattzeit cycle are thought to modulate the century-type cycle of Gleissberg, and the century-type cycle of Gleissberg is thought to modulate the 11 yr Schwabe cycle (Peristykh and Damon, 2003).In this article, ∼ 100, 200 and 2200 yr solar variations  are also thought to modulate EAM variations at ∼ 30, 60 and 600 yr, respectively.Moreover, solar variations at ∼ 100, 200, 2000 yr time scales are as large as at the other time scales in both spectral and EMD analysis.So, the changes in solar variation at these three time scales may amplify the changes in EAM at 30, 60 and 600 yr time scales.By this way, small changes in solar activity at 30, 60 and 600 yr may be able to trigger larger changes in EAM.Therefore, solar variability may be a plausible direct driver of EAM at these three time scales.

Conclusions
EAM and solar activity share some common cycles on ∼ 100, 200, 1000 and 2000 yr time scales.This supports the idea that solar changes are partly responsible for changes in Holocene EAM intensity from centennial to millennial scales (Dykoski et al., 2005;Wang et al., 2005;Cosford et al., 2008).Moreover, EAM and solar activity have much higher correlation at millennial time scales than at centennial time scales.This result also supports the assumption that a direct response to the solar variation holds more likely for time scales longer than a few hundred years (Weber et al., 2004).
The 10 Be and 14 C records do not reveal cycles at ∼30, 60 and 600 yr, and, hence, it seems that there is no direct relation between EAM and solar activity at these time scales.However, EAM intensity at these time scales is amplitude modulated by solar activities at ∼ 100, 220, and 2200 yr, which implies that solar activity is a plausible driver of EAM at these time scales.So, EAM variation at ∼ 30, 60 and 600 yr time scales may be directly amplified by solar activity.It may be a possible mechanism to explain why small changes in solar output at these time scales can cause large changes in EAM.

Fig. 4 .
Fig. 4. Power spectrum using REDFIT3.6 (Schulz and Mudelsee, 2002) for δ 18 O (a), 14 C (b) and 10 Be (c) records.Default parameters are used except that the parameter of n50 representing the Welch's overlap segment averaging is 4 and Welch spectrum window is chosen (refer to Schulz and Mudelsee, 2002 for details). Fig.5

Fig. 6 .
Fig. 6.Comparisons between the envelopes of IMFs of δ 18 O and IMFs of 10 Be.(a) IMF1 envelope of δ 18 O and IMF2 of 10 Be.(b) IMF2 envelope of 1δ 18 O and IMF3 of 10 Be.(c) IMF5 envelope of δ 18 O and IMF6 of 10 Be.

Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union. H. Y. Liu et al.: Possible link between Holocene East Asian monsoon and solar activity contained
. Spectral analysis is a global analysis that gives the frequency components in the signal.It lacks the time localization of the spectral components.This type of analysis works best when the input signal is linear and stationary

Table 1 .
Cycles of δ 18 O, 14 C and 10 Be with the EMD method.

Table 2 .
Correlation coefficients(r) of IMFs between 14 C and 10 Be.