Natural periodicities and Northern Hemisphere – Southern Hemisphere connection of fast temperature changes during the last glacial period : EPICA and NGRIP revisited

We investigate both the European Project for Ice Coring in Antarctica Dronning Maud Land (EDML) and North Greenland Ice-Core Project (NGRIP) data sets to study both the time evolution of the so-called Dansgaard–Oeschger events and the dynamics at longer timescales during the last glacial period. Empirical mode decomposition (EMD) is used to extract the proper modes of both the data sets. It is shown that the time behavior at the typical timescales of Dansgaard–Oeschger events is captured through signal reconstructions obtained by summing five EMD modes for NGRIP and four EMD modes for EDML. The reconstructions obtained by summing the successive modes can be used to describe the climate evolution at longer timescales, characterized by intervals in which Dansgaard–Oeschger events happen and intervals when these are not observed. Using EMD signal reconstructions and a simple model based on the one-dimensional Langevin equation, it is argued that the occurrence of a Dansgaard–Oeschger event can be described as an excitation of the climate system within the same state, while the longer timescale behavior appears to be due to transitions between different climate states. Finally, on the basis of a cross-correlation analysis performed on EMD reconstructions, evidence that the Antarctic climate changes lead those of Greenland by a lag of ≈ 3.05 kyr is presented.


Introduction
Some time ago, it was suggested that the climate of the Holocene period is much more variable than believed, being perhaps part of a millennial-scale pattern (Denton and Karlén, 1973;O'Brien et al., 1995;Bond et al., 1997).Oxygen isotope δ 18 O data from Greenland ice cores (North Greenland Ice Core Project members, 2004) reveal that the sub-Milankovitch climate variability presents abrupt temperature fluctuations which dominated the climate in Greenland between 11 and 74 thousand years before present (kyr BP).These fluctuations, known as Dansgaard-Oeschger (DO) events, are characterized by fast warmings, lasting few decades, with temperature increasing up to 15 • C compared to glacial values, followed by plateau phases and slow coolings of several centuries.
The origin of DO events is usually attributed to mechanisms involving the coupling between the Atlantic Ocean thermohaline circulation (THC), changes in atmospheric circulation and sea-ice cover (e.g., Rahmstorf, 2002).THC bistability (Broecker et al., 1985), internal oscillations in the THC volume transport (Broecker et al., 1990), and latitude shifts of oceanic convection (Rahmstorf, 1994) are some of the most popular ideas developed in this context.Another type of abrupt climate change found in paleoclimatic records is represented by the so-called Heinrich events (Heinrich, 1988).These are cooling events which are recorded in North Atlantic Ocean sediments and are characterized by variable spacings of several thousand years.Heinrich events are associated with massive iceberg releases from the Laurentide ice sheet into the North Atlantic (e.g., Bond et al., 1992;Broecker, 1994).Some relations between DO and Heinrich events have been highlighted in previous works (e.g., Bond et al., 1992;Rahmstorf, 2002).
Recurrent events known as Antarctic Isotope Maxima (AIM), less pronounced than DO, have been recognized also in Antarctic records, even if they are characterized by much more gradual warming and cooling with respect to DO (EPICA-members, 2004).Evidence has been presented, in previous works, in support of both synchronous evolution between north and south (e.g., Bender et al., 1994) and the existence of systematic lags between Antarctica and Greenland warming trends, with the first leading the second by 1-2.5 kyr (Blunier et al., 1998).This findings have been discussed in several works in the framework of the interhemispheric coupling.According to the classical bipolar seesaw hypothesis (Broecker, 1998;Stocker, 1998), an antiphase relation should exist between north and south.However, it was shown more recently that thermal bipolar seesaw models, obtained by introducing the effect of a thermal reservoir, are able to reproduce much of the variability observed in ice core records (Stocker et al., 2003;Barker et al., 2011).
The claimed roughly periodic, non-sinusoidal character of DO events, with a basic period of about T 1450-1500 yr, has been questioned on the basis of the fact that a null random hypothesis for these events cannot be fully rejected (Ditlevsen et al., 2007;Schulz, 2002;Peavoy and Franzke, 2010).This is obviously due to the fact that the criteria used for a definition of DO events cannot be firmly established.The usual numbered DO events were firstly determined visually (Dansgaard et al., 1993), but when different definitions are used, some events are excluded from the list or additional events are included, thus changing the basic periodicity (Ditlevsen et al., 2007;Schulz, 2002;Rahmstorf, 2003;Alley et al., 2001;Ditlevsen et al., 2005;Bond et al., 1999).For example, Schulz (2002) showed that the basic cycle is determined solely by DO events 5, 6, and 7, even if a fundamental period of about 1470 years seems to control the timing of the onset of the events.On the other hand, different periodicities in the range 1-2 kyr have been found from deep-sea sediment cores in the sub-polar North Atlantic (Bond et al., 1999).Of course, the main problem comes from the non-stationarity of the oxygen isotope data sets (Ditlevsen et al., 2007;Schulz, 2002), and advanced signal processing tools turned out to be useful to characterize the observed climatic variability (e.g., Solé et al., 2007b).
In this paper, we present a study of climate variability in Antarctica and Greenland, based on the analysis of two data sets coming from the European Project for Ice Coring in Antarctica (EPICA) and North Greenland Ice-Core Project (NGRIP).The empirical mode decomposition technique (EMD) is used to recognize the dominant variability patterns and characterize the climate transitions which can be associated with DO events and longer timescale changes.
A cross-correlation analysis is also performed on EMD reconstructions with the aim of identifying the existence of correlation lags at different timescales between the two signals.The paper is organized as follows.Section 2 is devoted to the description of the data sets, the EMD technique and the results of its application on the time series under consideration.The potential analysis and the characterization of climate transitions are illustrated in Sect.3. The crosscorrelation analysis and the results about north-south asynchrony are described in Sect. 4. Discussion and conclusions are given in Sect. 5.

Data sets
The EPICA data set provides a record of the oxygen isotope δ 18 O from the EPICA Dronning Maud Land (EDML) Ice Core (75.00 • S, 0.07 • E; 2892 m a.s.l.) covering the period 0-150 kyr BP (EPICA-members, 2006(EPICA-members, , 2010)).Data for the Northern Hemisphere come from the NGRIP project (75.10 • N, 42.32 • W; 2917 m elevation) and extend back to 123 kyr BP (North Greenland Ice Core Project members, 2004).We consider here the interval 20-120 kyr BP for both the data sets, since this is the interval in which significant temperature changes, which are the focus of the present work, are observed (see Fig. 1).Both the data sets are synchronized using the AICC2012 age scale (Veres et al., 2013;Bazin et al., 2013).

Empirical mode decomposition
To identify the main periodicities and their amplitudes, we applied the empirical mode decomposition (EMD), a technique designed to investigate non-stationary data (Huang et al., 1998).The EMD has been extensively and successfully used in different fields (Cummings et al., 2004;Terradas et al., 2004;Franzke, 2009;Vecchio et al., 2010aVecchio et al., , b, 2012;;Capparelli et al., 2011).In the context of paleoclimatic studies, the EMD was used by Solé et al. (2007a) to investigate the oscillation patterns in GRIP, Vostok and EPICA time series of temperature proxies.
In the present work, the EDML and NGRIP δ 18 O time series are decomposed, by means of the EMD, into a finite number m of oscillating intrinsic mode functions (IMFs) as (1) The functional form of the IMFs, C j (t), is not given a priori, but obtained from the data through the algorithm developed by Huang et al. (1998).As a first step, the local extrema, minima and maxima, are identified in the raw data δ 18 O.The local extrema are interpolated by means of cubic splines to obtain the envelopes of the maxima and minima.Then the To avoid a loss of information about amplitude and frequency modulations, a criterion to stop the above described procedure has been proposed (Huang et al., 1998) by defining calculated between two consecutive iterations.The iterations are stopped when σ is smaller than a threshold σ th , which is typically fixed at 0.3 (Huang et al., 1998) where P denotes the Cauchy principal value and C * j (t) is the complex conjugate of C j (t), the instantaneous phase can be calculated as φ(t) = arctan[C * j (t)/C j (t)] and the instantaneous frequency ω j (t) and period T j (t) are given by 190 ω j (t) = 2π/T j (t) = dφ j /dt.The characteristic period T j of each IMF can be estimated as the time average of T j (t).
Since the decomposition is local, complete, and orthogonal, the EMD can be used as a filter by reconstructing the signal through partial sums in Eq. (1) (Huang et al., 1998;195 Cummings et al., 2004;Vecchio et al., 2010a).
The use of empirical EMD functions, characterized by time-dependent amplitude and phase, allows to overcome some limitations of Fourier analysis.The latter requires linear systems and periodic or stationary data, and its appli-200 cation to non-linear, non-stationary data can produce misleading results for the following reasons.(1) The Fourier uniform harmonic components do not carry local information: many components are needed to build up a solution that corresponds to non-stationary data thus resulting in a 205 energy spreading over a wide frequency range.As a consequence the energy-frequency distribution of non-linear and non-stationary data is not accurate.(2) Fourier spectral analysis uses linear superposition of sinusoidal functions, thus several components are mixed together in order to reproduce 210 deformed waveforms, local variations, or the fictitious periodic boundary conditions imposed by the analysis.(3) Sinusoidal functions are usually far from being eigenfunctions of the phenomenon under study for non-linear/non-stationary data.In these situations, when the data are far to be peri-215 odic, linear and stationary, an empirical decomposition such as the EMD provides a better description of the analysed phenomenon.
The statistical significance of the IMFs with respect to a white noise can be verified through the test developed by 220 Wu and Huang (2004) that is based on the constancy of the product between the mean square amplitude of each IMF and the corresponding average period when the EMD is applied to a white noise series.This allows to derive, for each IMF, the analytical mean square amplitude spread function for dif- difference h 1 (t) between δ 18 O and m 1 (t), the latter being the mean between the envelopes of the maxima and minima, is computed.If this quantity satisfies two conditions -(i) the number of extrema and zero crossings does not differ by more than 1 and (ii) the mean value of the envelopes obtained from the local maxima and minima is zero at any point -it represents an IMF.If the conditions above are not fulfilled, the described procedure is repeated on the h 1 (t) time series, and h 11 (t) = h 1 (t)−m 11 (t), where m 11 (t) is the mean of the new envelopes, is calculated.This process is iterated until, after s times, h 1s (t) fulfils the IMF properties.To avoid a loss of information about amplitude and frequency modulations, a criterion to stop the above described procedure has been proposed (Huang et al., 1998) by defining calculated between two consecutive iterations.The iterations are stopped when σ is smaller than a threshold σ th , which is typically fixed at 0.3 (Huang et al., 1998).When the first IMF C 1 (t) is obtained, the "first residue" r 1 (t) = δ 18 O -C 1 (t) is calculated and processed in the same way as previously explained.A new IMF C 2 (t) and a second residue r 2 (t) are then obtained.This procedure is carried on until C s or r s is almost constant or when the residue r s (t) is monotonic.
As a result of this method, the original signal is decomposed into m empirical modes, ordered in increasing characteristic timescale, and a residue r m (t) which provides the mean trend, if any exists, in the data set.Each C j (t) represents a zero mean oscillation C j (t) = A(t) sin φ(t) (being φ(t) the instantaneous phase) experiencing modulation both in amplitude and frequency.Performing the Hilbert transform of each IMF where P denotes the Cauchy principal value and C * j (t) is the complex conjugate of C j (t), the instantaneous phase can be calculated as and the instantaneous frequency ω j (t) and period T j (t) are given by ω j (t) = 2π/T j (t) = dφ j /dt.The characteristic period T j of each IMF can be estimated as the time average of T j (t).
Since the decomposition is local, complete, and orthogonal, the EMD can be used as a filter by reconstructing the signal through partial sums in Eq. ( 1) (Huang et al., 1998;Cummings et al., 2004;Vecchio et al., 2010a).
The use of empirical EMD functions, characterized by time-dependent amplitude and phase, allows overcoming some limitations of Fourier analysis.The latter requires linear systems and periodic or stationary data, and its application to nonlinear, non-stationary data can produce misleading results for the following reasons.(1) The Fourier uniform harmonic components do not carry local information: many components are needed to build up a solution that corresponds to non-stationary data -thus resulting in energy spreading over a wide frequency range.As a consequence the energy-frequency distribution of nonlinear and non-stationary data is not accurate.(2) Fourier spectral analysis uses linear superposition of sinusoidal functions; thus several components are mixed together in order to reproduce deformed waveforms, local variations, or the fictitious periodic boundary conditions imposed by the analysis.(3) Sinusoidal functions are usually far from being eigenfunctions of the phenomenon under study for nonlinear/non-stationary data.In these situations, when the data are far to be periodic, linear and stationary, an empirical decomposition such as the EMD provides a better description of the analyzed phenomenon.
The statistical significance of the IMFs with respect to white noise can be verified through the test developed by Wu and Huang (2004) that is based on the constancy of the product between the mean square amplitude of each IMF and the corresponding average period when the EMD is applied to a white noise series.This allows the derivation, for each IMF, of the analytical mean square amplitude spread function for different confidence levels.Therefore, by comparing the mean square amplitude of the EMD modes obtained from the data under analysis to the theoretical spread function, the IMFs containing information at a given confidence level can be discerned from purely noisy IMFs.

Empirical mode decomposition results
Applying the EMD procedure to the EDML and NGRIP data, we obtain a set of m = 15 IMFs for EDML (see Fig. 2), which we denote as C (S) j (t), and m = 17 IMFs, named C (N) j (t), for NGRIP (see Fig. 3).Some IMFs for both data sets, e.g., C (S) 11 (t), show a "mode mixing" behavior consisting of signals containing different frequencies or a signal of a similar timescale residing in different IMFs (Huang et al., 1998).Mode mixing, a consequence of signal intermittency, is troublesome in signal processing where the main purpose is the signal cleanliness.However, this effect is not crucial in our analysis since in the following we will discuss signals obtained through partial sums by adding IMFs ranging in a wide range of timescales.
The results of the significance test are shown in Fig. 4. In performing the test, it was assumed that the energy of the IMFs j = 0 comes only from noise, and it is thus assigned to the 99 % line (Wu and Huang, 2004).The modes which are above the spread line can be considered significant at 99th percentile.The significant modes are j = 5-14 and j = 4-16 for the EDML and NGRIP data respectively.
The characteristic periods T j are reported in Table 1.
The dynamics of the DO events is reconstructed by the sum of the j = 6-10 NGRIP IMFs (which have characteristic periods between 0.7 kyr and 3.3 kyr).For EDML, using the same characteristic period range as NGRIP (0.7-3.3 kyr), the modes j = 5-8 are selected.Therefore the two "short" timescale reconstructions N H (t) and S H (t), for NGRIP and EDML respectively, are defined as In order to investigate the variability at longer timescales with respect to those involved in DO events, the modes j = 11-16 and j = 9-14 are used for NGRIP and EDML respectively.The corresponding two "long" timescale reconstructions N L (t) and S L (t) are thus given by In Fig. 5 the δ 18 O original data (black lines), the short timescale (red lines) and the long timescale (blue lines) reconstructions are shown for the EDML (upper panel) and NGRIP (lower panel) data sets.
The results of the EMD decomposition suggest that the evolution of cooling-warming cycles over the investigated period can be described in terms of two dynamical processes occurring on different timescales.A simple model of such a phenomenology is considered in the next section.

Potential analysis
Glacial-interglacial cycles are commonly modeled as transitions between different climate states (see, e.g., Paillard, 2001).The same approach has been extended to DO events that are commonly modeled as a two-state, cold/stadial and warm/interstadial, system.In the following we apply the method developed by Livina et al. (2010) to identify the number of climate states present in both NGRIP and EDML records.To this purpose, the climate system is described in terms of a nonlinear system with many dynamical states, and we assume that transitions among states are triggered by a stochastic forcing.A very simple model for this is the onedimensional Langevin equation (Livina et al., 2010): where z represents, in our case, the oxygen isotope δ 18 O, U (z) is the potential, σ is the noise level and W is a Wiener process.A stationary solution of Eq. ( 8) is found when U (z) is a polynomial of even order and positive leading coefficient.The order of the polynomial fixes the number of available states can be evaluated from the data by performing a polynomial fit of the probability density function (pdf) calculated as 305 representing a stationary solution of the Fokker-Planck equation Equation ( 9) establishes a one to one correspondence between the potential and the the pdf of the system so that where p emp (z) is the empirical pdf extracted from data.To estimate p emp (z) we used the well known Kernel Density Estimator method (Silverman, 1998;Hall, 1992), described in Appendix A. This procedure allows to calculate U (z) and 315 the corresponding uncertainty from Equation (11).In order to investigate the time evolution of the DO events, potentials have been calculated for N H (t), S H (t), N L (t), S L (t), and reported in Fig. 6.The potentials associated with N H (t), S H (t) and the corresponding best fits are shown in 320 panels a and b respectively.The best fit polynomials are of 2 nd order, thus corresponding to a single well potential.Therefore, the occurrence of a DO event could not be due to a transition of the system to a different dynamical state but states: for example, a second-order polynomial corresponds to a system with single-well potential and one state, while a fourth-order polynomial, corresponding to a double-well potential, identifies a system with two states.The number of states can be evaluated from the data by performing a polynomial fit of the probability density function (pdf) calculated as representing a stationary solution of the Fokker-Planck equation Equation ( 9) establishes a one-to-one correspondence between the potential and the the pdf of the system so that where p emp (z) is the empirical pdf extracted from data.To estimate p emp (z) we used the well-known kernel density estimator method (Silverman, 1998;Hall, 1992), described in Appendix A. This procedure allows the calculation of U (z) and the corresponding uncertainty from Eq. ( 11).In order to investigate the time evolution of the DO events, potentials have been calculated for N H (t), S H (t), N L (t), and S L (t), and reported in Fig. 6.The potentials associated with N H (t), S H (t) and the corresponding best fits are shown in panels a and b respectively.The best fit polynomials are of second order, thus corresponding to a single-well potential.to an excitation of the system within the same state.On the other hand, the U (z) calculated for N L (t), S L (t) are well fitted by 4 th order polynomials, indicating the occurrence of a double-well potential, thus suggesting that the high activity, when the DO events are observed, and low activity periods correspond to different states of the climatic system.
It is worth to briefly discuss the differences of our approach with respect to Livina et al. (2010).We perform the potential analysis on the EMD reconstructions of the EDML and NGRIP data-sets, using the time range 20-120 kyr BP.Livina et al. (2010) used only Greenland data-sets, in particular they considered the GRIP and NGRIP δ 18 O series in the interval 0-60 kyr BP.Moreover we calculate the potential shape from the full time range of the EMD reconstructions, while they used sliding windows of varying length through each data-set to detect the number of climate states as a func-340 tion of time.

The North-South asyncrony
An important aspect of the polar climate dynamics is the understanding of how Earth's hemispheres have been coupled during past climate changes.Bender et al. (1994) explored 345 the possible connections between Greenland and Antarctica climate, suggesting that partial deglaciation and changes in ocean circulation are the main mechanisms which transfer Therefore, the occurrence of a DO event could not be due to a transition of the system to a different dynamical state but to an excitation of the system within the same state.On the other hand, the U (z) calculated for N L (t) and S L (t) is well fitted by fourth-order polynomials, indicating the occurrence of a double-well potential, thus suggesting that the high activity, when the DO events are observed, and low activity periods correspond to different states of the climatic system.
It is worth briefly discussing the differences of our approach with respect to Livina et al. (2010).We perform the potential analysis on the EMD reconstructions of the EDML and NGRIP data sets, using the time range 20-120 kyr BP.Livina et al. (2010) used only Greenland data sets; in par-ticular they considered the GRIP and NGRIP δ 18 O series in the interval 0-60 kyr BP.Moreover we calculate the potential shape from the full time range of the EMD reconstructions, while they used sliding windows of varying length through each data set to detect the number of climate states as a function of time.

The north-south asynchrony
An important aspect of the polar climate dynamics is the understanding of how earth's hemispheres have been coupled during past climate changes.Bender et al. (1994)  warmings in Northern Hemisphere climate to Antarctica.Different analyses showed that southern changes are not a 350 direct response to abrupt changes in North Atlantic thermohaline circulation, as is assumed in the conventional picture of a hemispheric temperature seesaw (Morgan et al., 2002).
The study of the global concentration of methane recorded in ice cores allowed to infer that Antarctic climate changes lead 355 that of Greenland by 1 − 2.5 kyr over the period 47-23 kyr before present (Blunier et al., 1998).More recently, Barker et al. (2011) observed a a near zero-phase anticorrelation between the Greenland temperature anomaly and the rate of change of Antarctic temperature.Based on their correlation 360 analysis and on the thermal bipolar seesaw model (Stocker et al., 2003), the same authors constructed, using the Antarctic record, an 800 kyr synthetic record of Greenland climate which reproduces much of the variability fo the last 100 kyr.
In order to investigate this issue, we study the cross-365 correlation between the EDML and NGRIP reconstructions obtained from the EMD as illustrated in Sect.2.3.The crosscorrelation coefficient P yz (∆) between two time samples y(t k ) and z((t k ) is defined as where brackets denote time averages and ∆ the time lag.We calculated the cross-correlation coefficient both for the short time-scale reconstructions (Eqs.( 4) and ( 5)) and for the long time-scale ones (Eqs.( 6) and ( 7)).The two coefficients P N H S H (∆) and P N L S L (∆) are shown in Fig. 7.
375 displays oscillations with many peaks of similar amplitude at both negative and positive lags.This behaviour is typically obtained when oscillating signals having nearly the same frequency are compared.In this case, it is not possible to identify the leading and the following process.On the other hand,  climate, suggesting that partial deglaciation and changes in ocean circulation are the main mechanisms which transfer warmings in Northern Hemisphere climate to Antarctica.Different analyses showed that southern changes are not a direct response to abrupt changes in North Atlantic thermohaline circulation, as is assumed in the conventional picture of a hemispheric temperature seesaw (Morgan et al., 2002).
The study of the global concentration of methane recorded in ice cores allowed Blunier et al. (1998) to infer that Antarctic climate changes lead those of Greenland by 1-2.5 kyr over the period 47-23 kyr before present.More recently, Barker et al. (2011) observed a near-zero-phase anticorrelation between the Greenland temperature anomaly and the rate of change of Antarctic temperature.Based on their correlation analysis and on the thermal bipolar seesaw model (Stocker et al., 2003), the same authors constructed, using the Antarctic record, an 800 kyr synthetic record of Greenland climate which reproduces much of the variability for the last 100 kyr.
In order to investigate this issue, we study the crosscorrelation between the EDML and NGRIP reconstructions obtained from the EMD as illustrated in Sect.2.3.The crosscorrelation coefficient P yz ( ) between two time samples warmings in Northern Hemisphere climate to Antarctica.Different analyses showed that southern changes are not a direct response to abrupt changes in North Atlantic thermohaline circulation, as is assumed in the conventional picture of a hemispheric temperature seesaw (Morgan et al., 2002).The study of the global concentration of methane recorded in ice cores allowed to infer that Antarctic climate changes lead that of Greenland by 1 − 2.5 kyr over the period 47-23 kyr before present (Blunier et al., 1998).More recently, Barker et al. (2011) observed a a near zero-phase anticorrelation between the Greenland temperature anomaly and the rate of change of Antarctic temperature.Based on their correlation analysis and on the thermal bipolar seesaw model (Stocker et al., 2003), the same authors constructed, using the Antarctic record, an 800 kyr synthetic record of Greenland climate which reproduces much of the variability fo the last 100 kyr.
In order to investigate this issue, we study the crosscorrelation between the EDML and NGRIP reconstructions obtained from the EMD as illustrated in Sect.2.3.The crosscorrelation coefficient P yz (∆) between two time samples Fig. 5. δ 18 O data (black lines), short time scale reconstructions N H (t) and S H (t) (red lines), and long time scale reconstructions N L (t) and S L (t) (blue lines) for the EDML (upper panel) and NGRIP (lower panel) datasets.An offset corresponding to the temporal mean of the δ 18 O original data was applied to the EMD reconstructions in order to allow visualization in the same plot.y(t k ) and z((t k ) is defined as 370 where brackets denote time averages and ∆ the time lag.We calculated the cross-correlation coefficient both for the short time-scale reconstructions (Eqs.( 4) and ( 5)) and for the long time-scale ones (Eqs.( 6) and ( 7)).The two coefficients P N H S H (∆) and P N L S L (∆) are shown in Fig. 7.
375 displays oscillations with many peaks of similar amplitude at both negative and positive lags.This behaviour is typically obtained when oscillating signals having nearly the same frequency are compared.In this case, it is not possible to identify the leading and the following process.On the other hand, a shape characterized by a single significant peak, with a maximum value of ≈ 0.73 and ∆ = 3.05 ± 0.19 kyr, is found in P N L S L (∆).The uncertainty on the correlation peak position was estimated by means of the following procedure.The y(t k ) and z((t k ) is defined as where brackets denote time averages and the time lag.We calculated the cross-correlation coefficient both for the short timescale reconstructions (Eqs.4 and 5) and for the long timescale ones (Eqs.6 and 7).The two coefficients P N H S H ( ) and P N L S L ( ) are shown in Fig. 7. P N H S H ( ) displays oscillations with many peaks of similar amplitude at both negative and positive lags.This behavior is typically obtained when oscillating signals having nearly the same frequency are compared.In this case, it is not possible to identify the leading and the following process.On the other hand, a shape characterized by a single significant peak, with a maximum value of ≈ 0.73 and = 3.05 ± 0.19 kyr, is found in P N L S L ( ).The uncertainty in the correlation peak position was estimated by means of the following procedure.The errors on the age scale available in the EDML and NGRIP 385 data were used in a Montecarlo algorithm to estimate the time errors on N L (t) and S L (t).More specifically, we calculated 10 3 realizations of the long time-scale reconstructions varying randomly the age scale position of each data point within the error windows.Then, we calculated the corresponding 390 10 3 cross-correlations between the EDML and NGRIP long time-scale reconstructions and the peak positions for each of them.Following this procedure we obtained the above mentioned estimate of 0.19 kyr for the error on the position of the long time-scale cross correlation peak.

395
The result obtained through the cross-correlation analysis supports the view according to which the Antarctic climate changes actually lead that of Greenland, but on a longer timescale than previously reported.

400
In the present paper we study both the EDML and NGRIP datasets, in order to investigate their periodicities and the possible existence of synchronization of abrupt climate changes observed in the North and South hemispheres.Our results can be summarized as follows.405 1. Proper EMD modes, significant with respect to a random null hypothesis, are present in both datasets, thus confirming that natural cycles of abrupt climate changes during the last glacial period exist and their occurrence cannot be due to random fluctuations in time.errors on the age scale available in the EDML and NGRIP data were used in a Monte Carlo algorithm to estimate the time errors on N L (t) and S L (t).More specifically, we calculated 10 3 realizations of the long timescale reconstructions varying randomly the age-scale position of each data point within the error windows.Then, we calculated the corresponding 10 3 cross-correlations between the EDML and NGRIP long timescale reconstructions and the peak positions for each of them.Following this procedure we obtained the above-mentioned estimate of 0.19 kyr for the error on the position of the long timescale cross-correlation peak.
The result obtained through the cross-correlation analysis supports the view according to which the Antarctic climate changes actually lead those of Greenland, but on a longer timescale than previously reported.

Conclusions and discussion
In the present paper we study both the EDML and NGRIP data sets, in order to investigate their periodicities and the possible existence of synchronization of abrupt climate changes observed in the Northern and Southern hemispheres.Our results can be summarized as follows.
1. Proper EMD modes, significant with respect to a random null hypothesis, are present in both data sets, thus confirming that natural cycles of abrupt climate changes during the last glacial period exist and their occurrence cannot be due to random fluctuations in time.
2. Due to the non-stationarity of the process, the variability at the typical timescales of DO events is reproduced through signal reconstructions obtained by summing more than one EMD mode, more precisely j = 6-10 for NGRIP and j = 5-8 for EDML.On the other hand, the reconstructions obtained summing the successive modes can be used to describe the climate evolution at longer timescales, characterized by intervals in which DO events happen and intervals when these are not observed.
3. By comparing the EMD signal reconstructions to the results of a simple model based on the one-dimensional Langevin equation, evidence is found that the occurrences of DO events can be described as excitations of the system within the same climate state, rather than transitions between different states.Conversely, the longer timescale dynamics appear to be due to transitions between different climate states.for NGRIP and j = 5−8 for EDML.On the other hand,

415
the reconstructions obtained summing the successive modes can be used to describe the climate evolution at longer time scales, characterized by intervals in which DO events happen and intervals when these are not observed.420 3.By comparing the EMD signal reconstructions to the results of a simple model based on the one-dimensional Langevin equations, evidence is found that the occurrence of DO events can be described as excitations of the system within the same climate state, rather 425 than transitions between different states.Conversely, the longer time scale dynamics appear to be due to transitions between different climate states.
4. On the base of a cross correlation analysis between the NGRIP and EDML EMD reconstructions, performed to 430 investigate the North-South asynchrony, it is found that the clearest correlation occurs between the long-scale reconstructions at a lag ∆ = 3.05 ± 0 − 19 kyr, which supports the view according to which the Antarctic climate changes lead that of Greenland, but on a longer 435 time-scale than previously reported.
The methodological novelty of the present work is represented by the fact that we use EMD reconstructions to investigate the climate dynamics at different time scales and to highlight, through a potential analysis of the EMD re-440 constructions, some characteristics of the climate transitions.The main new results are: 1) the presence of two different climate transitions occurring at different time scales; 2) the finding of a significant correlation between long time scale Antarctic and Greenland signals, with the first leading the 445 second, at a lag of ≈ 3 kyr which was not found in past studies of paleoclimatic records.
The EMD filtering procedure applied on EDML and NGRIP datasets and the correlation analysis based on it are able to identify time scale dependent dynamical features of 450 the climate evolution which have not been underlined in previous works.How the interhemispheric coupling mechanisms, such as the THC, can be involved in the behavior highlighted by our study is a question which needs to be investigated in future works.We suggest that the results of 455 our correlation analysis results, and in particular the correlation found in association with the long time scale dynamics, could be explained in the framework of seesaw models.But, since the correlation lag (≈ 3 kyr) obtained from our analysis is quite different from the characteristic thermal time 460 scale (about 1 -1.5 kyr) of previous bipolar seesaw models (Stocker et al., 2003;Barker et al., 2011), it would be necessary to build up a thermal seesaw model starting from our EMD filtered long time scale series to properly deal with this problem.

465
We also plan to extend the study presented in this paper to other paleoclimate records, in order to investigate the role of other physical processes, such ocean-atmosphere interaction and solar irradiance variations, in the climate system evolution.Moreover, the results presented in this paper could 470 be useful for the theoretical modelling of the climate evolution in polar regions, in order to characterise the processes involved at different time scales and the coupling between northern and southern hemispheres.

Kernel Density Estimator
The empirical pdf p emp (z) in Equation ( 11) is estimated through the Kernel Density Estimator technique (Silverman, 1998;Hall, 1992).The Kernel Density Estimator is defined 4. On the basis of a cross-correlation analysis between the NGRIP and EDML EMD reconstructions, performed to investigate the north-south asynchrony, it is found that the clearest correlation occurs between the long-scale reconstructions at a lag = 3.05 ± 0.19 kyr, which supports the view according to which the Antarctic climate changes lead those of Greenland, but on a longer timescale than previously reported.
The methodological novelty of the present work is represented by the fact that we use EMD reconstructions to investigate the climate dynamics at different timescales and to highlight, through a potential analysis of the EMD reconstructions, some characteristics of the climate transitions.The main new results are (1) the presence of two different climate transitions occurring at different timescales as well as (2) the finding of a significant correlation between long timescale Antarctic and Greenland signals, with the first leading the second, at a lag of ≈ 3 kyr, which was not found in past studies of paleoclimatic records.
The EMD filtering procedure applied on EDML and NGRIP data sets and the correlation analysis based on it are able to identify timescale-dependent dynamical features of the climate evolution which have not been underlined in previous works.How the interhemispheric coupling mechanisms, such as the THC, can be involved in the behavior highlighted by our study is a question which needs to be investigated in future works.We suggest that the results of our correlation analysis results, and in particular the correlation found in association with the long timescale dynamics, could be explained in the framework of seesaw models.But, since the correlation lag (≈ 3 kyr) obtained from our analysis is quite different from the characteristic thermal timescale (about 1-1.5 kyr) of previous bipolar seesaw models (Stocker et al., 2003;Barker et al., 2011), it would be necessary to build up a thermal seesaw model starting from our EMD filtered long timescale series to properly deal with this problem.
We also plan to extend the study presented in this paper to other paleoclimate records, in order to investigate the role of other physical processes, such as ocean-atmosphere interaction and solar irradiance variations, in the climate system evolution.Moreover, the results presented in this paper could be useful for the theoretical modeling of the climate evolution in polar regions, in order to characterize the processes involved at different timescales and the coupling between the Northern and Southern hemispheres.The empirical pdf p emp (z) in Eq. ( 11) is estimated through the kernel density estimator technique (Silverman, 1998;Hall, 1992).Given a set of data points z i (i = 1, ..., n), the kernel density estimator of the pdf is defined as where K is the kernel, a symmetric function that integrates to one, and h > 0 is a smoothing parameter denoted as bandwidth.The corresponding variance is For our application we use the Epanechnikov kernel, which is optimal in a minimum variance sense, and we chose h = 0.9 • s/n 0.2 , where s is the standard deviation of the data set.The confidence interval of f (z) has been evaluated through a bootstrap.Being the bootstrap t statistic (Hall, 1992), the symmetric confidence interval with coverage probability 1 where u * α/2 is the bootstrap estimate of the quantile defined by P (t * (z) ≤ u * α/2 ) = α/2 and u * 1−α/2 is the bootstrap estimate of the quantile defined by P (t * (z) ≤ u * 1−α/2 ) = 1 − α/2 (Hall, 1992).In order to remove asymptotic biases which are not correctly taken into account by the bootstrap procedure, we perform an under-smoothing by choosing a smaller bandwidth h * = 0.9 * s/n 0.25 < h (Hall, 1992).
This procedure allows the calculation of p emp (z), U (z) from Eq. ( 11) and the corresponding uncertainty by performing error propagation in the same equation.

Figure 1 .
Figure 1.δ 18 O data for the EDML (upper panel) and NGRIP (lower panel) data sets in the period 20-120 kyr BP

Figure
Figure 2. IMFs C (S) j (t) and residue r (S) m (t) for the EDML data set.

Figure
Figure 3. IMFs C (N ) j (t) and residue r (N ) m (t) for the NGRIP data set.

Fig. 4 .
Fig. 4. Normalized IMF square amplitude E j vs. the period T j for the EMD significance test applied to the EDML (upper panel) and NGRIP (lower panel) IMFs.The dashed lines represent the 99th percentile (see the text for details).

Fig. 5 .
Fig.5.δ 18 O data (black lines), short time scale reconstructions N H (t) and S H (t) (red lines), and long time scale reconstructions N L (t) and S L (t) (blue lines) for the EDML (upper panel) and NGRIP (lower panel) datasets.An offset corresponding to the temporal mean of the δ 18 O original data was applied to the EMD reconstructions in order to allow visualization in the same plot.
380a shape characterized by a single significant peak, with a maximum value of ≈ 0.73 and ∆ = 3.05 ± 0.19 kyr, is found in P N L S L (∆).The uncertainty on the correlation peak position was estimated by means of the following procedure.The

Figure 4 .
Figure 4. Normalized IMF square amplitude E j vs. the period T j for the EMD significance test applied to the EDML (upper panel) and NGRIP (lower panel) IMFs.The dashed lines represent the 99th percentile (see the text for details).

Fig. 4 .
Fig. 4. Normalized IMF square amplitude E j vs. the period T j for the EMD significance test applied to the EDML (upper panel) and NGRIP (lower panel) IMFs.The dashed lines represent the 99th percentile (see the text for details). 380

Figure 5 .
Figure5.δ 18 O data (black lines), short timescale reconstructions N H (t) and S H (t) (red lines), and long timescale reconstructions N L (t) and S L (t) (blue lines) for the EDML (upper panel) and NGRIP (lower panel) data sets.An offset corresponding to the temporal mean of the δ 18 O original data was applied to the EMD reconstructions in order to allow visualization in the same plot.

Figure 7 .
Figure 7. Cross-correlation coefficients between EDML and NGRIP short timescale (P N H S H ( ), top panel) and long timescale (P N L S L ( ), bottom panel) EMD reconstructions as functions of the time lag .

Table 1 .
Characteristic periods of the significant IMFs obtained for the EDML and NGRIP data sets through the EMD.Errors are estimated as standard deviations of the instantaneous periods.The period T 16 calculated for the j = 16 NGRIP mode is not reported as it is not sufficiently shorter than the time series length, although the mode is significant.