Natural periodicities and north – south hemispheres connection of fast temperature changes during the last glacial period : EPICA and NGRIP revisited

Introduction Conclusions References


Introduction
Some time ago, it has been 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) Antarctic, even if they are characterized by rather gradual warming and cooling (EPICAmembers, 2004).The claimed roughly periodic, non-sinusoidal character of these 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 dataset (Ditlevsen et al., 2007;Schulz, 2002).
DO events are widely described through models involving coupling between the Atlantic ocean thermohaline circulation and changes in atmospheric circulation and seaice cover (e.g.Rahmstorf, 2003).In particular, ice volume plays a relevant role in affecting the ocean-atmosphere system, as suggested by the fact that the Holocene interglacial period has not shown the extreme DO variability typical of the last ice age.
In this paper, we present a study of DO events in the two hemispheres, based on the analysis of the European Project for Ice Coring in Antarctica and North Greenland Ice-Core Project datasets.The Empirical Mode Decomposition technique is used to identify and characterize observed climate variability on different time scales.Introduction

Conclusions References
Tables Figures

The 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., 2010a, b;Capparelli et al., 2011;Vecchio et al., 2012).Both the δ 18 O time series are decomposed, by means of the EMD, into a finite number m of oscillating intrinsic mode functions (IMFs) as 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 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 Introduction

Conclusions References
Tables Figures

Back Close
Full satisfies the following 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 are 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 time scale, and a residue r m (t) which provides the mean trend, if any exists in the dataset.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 A(t) and frequency ω = dφ/dt.A typical average period T j can be estimated for all the IMFs.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 statistical significance of the IMFs with respect to a white noise has been verified through the test developed by Wu and Huang (2004).This is based on the constancy of the product between the mean square amplitude of each IMF and the corresponding Introduction

Conclusions References
Tables Figures

Back Close
Full 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 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 EPICA and NGRIP data, we obtain a set of m = 6 IMFs for EPICA (see Fig. 2), which we denote as C (S) j (t), and m = 7 IMFs, named C (N) j (t), for NGRIP (see Fig. 3).We remark that while r m (t) for the NGRIP dataset represents a real monotonic trend, for the EPICA sample it consists of an oscillating function with an offset representing the mean value of the δ 18 O proxy.Therefore, for EPICA we will consider r m (t), after subtracting its temporal mean r m (t) , as a genuine IMF.
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.The modes which are above the spread line can be considered significant at 99th percentile.For EPICA the modes j = 2-5 are significant, while the mode j = 1 lies just on the spread line.For NGRIP all the modes j = 1-6 result to be significant.
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 φ(t) = arctan[C * j (t)/C j (t)] and the instantaneous frequency ω j (t) and period T j (t) are given by ω j (t) = 2π/T j (t) = dφ j /d t.
The characteristic periods T j are estimated as the time averages of T j (t) and reported in Table 1.
More than one IMF is needed to reproduce the observed DO events.In fact, the main DO events are captured by the sum of the j = 1, 2, 3 modes, as shown in Fig. 5 where we report the superposition of the raw δ 18 O data for both the EPICA and NGRIP datasets, and the two reconstructions given by In order to investigate the variability at longer time scales with respect to those involved in DO events, we can consider the IMFs j = 4, 5, 6.Long term reconstructions, defined as, are reported as blue lines in Fig. 5. Reconstructions with long period IMFs j = 4-6 underline the long term trend, characterized by intervals in which DO events happen and intervals when these are not observed.
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 time-scales.A simple model of such a phenomenology is considered in the next section.Introduction

Conclusions References
Tables Figures

Back Close
Full

Potential analysis
Glacial-interglacial cycles are commonly modelled as transitions between different climate states (see e.g.Paillard, 2001).The same approach has been extended to DO events that are commonly modelled 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 EPICA records.To this purpose, the climate system is described in term 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 one-dimensional 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: 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 Kernel Density estimator (Silverman, 1998;Hall, 1992) where K is the kernel, a symmetric function that integrates to one, z i are the data points, and h > 0 is a smoothing parameter denoted as bandwidth.The corresponding variance is: For the following 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  (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 to calculate U(z) with the corresponding uncertainty from Eq. ( 11).In order to investigate the time evolution of the DO events, potentials have been separately calculated for N H (t), S H (t) and 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 panels a and b respectively.The best fit polynomials are of 2nd 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 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 4th 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.Introduction

Conclusions References
Tables Figures

Back Close
Full explored the possible connections between Greenland and Antarctica 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 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).
In order to investigate this issue, we study the cross-correlation between the EPICA and NGRIP reconstructions obtained from the EMD as illustrated in Sect.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 crosscorrelation 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. P N H S H (∆) 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 peak which is significantly higher than the others, with a maximum value of ≈ 0.66, is found in P N L S L (∆) at ∆ ≈ 3.6 kyr.This result supports the view according to which the Antarctic climate changes actually lead that of Greenland, but on a longer time-scale than previously reported.Introduction

Conclusions References
Tables Figures

Back Close
Full

Conclusions
In the present paper we study both the EPICA 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.
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.
2. Due to the non-stationarity of the process, the dynamics of DO events is captured by three EMD modes, j = 1-3 for both the datasets.Therefore, the signal reconstructions obtained by summing the modes j = 1-3 reproduce the variability at the typical time scales of DO events.On the other hand, the reconstructions with the modes j = 4-6 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.
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 occurrence of a DO event can be described as an excitation of the system within the same climate state, rather than a transition between different states.Conversely, the longer time scale dynamics appear to be due to transitions between different climate states.Introduction

Conclusions References
Tables Figures

Back Close
Full changes lead that of Greenland, but on a longer time-scale than previously reported.
We 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 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.

Conclusions References
Tables Figures

Back Close
Full  Full    (t), remark that while r m (t) for eal monotonic trend, for the oscillating function with an e of the δ 18 O proxy.Therer r m (t), after subtracting its uine IMF.
test are shown in Figure 4. umed that the energy of the ise and it is thus assigned to h are above the spread line 99-th percentile.For EPICA nt, while the mode j = 1 lies RIP all the modes j = 1 − 6 where K is the kernel, a symmetric function that to one, z i are the data points, and h > 0 is a smoo rameter denoted as bandwidth.The corresponding 245 is: For the following application we use the Epan kernel, which is optimal in a minimum variance we chose h = 0.9 * s/n 0.2 , where s is the standard 250 of the data set.The confidence interval of f (z) . 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 B.P.).These fluctuations, known as Dansgaard-Oeschger (DO) events, are characterized by fast warmings and slow coolings, with temperature increasing up to 15 • C compared to glacial values.Less pronounced recurrent events, known as Antarctic Isotope Maxima (AIM), have been recognized over the same period in Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Dansgaard/Oeschger cycles and the Little Ice Age, Mechanisms of Global Climate Change at Millennial Time Scales, Geophysical Monograph 112 (AGU), 1999.1131 Capparelli, V., Vecchio, A., and Carbone, V.: Long-range persistence of temperature records induced by long-term climatic phenomena, Phys.Rev. E, 84, 046103, doi:10.1103/PhysRevE.84.046103, 2011.Discussion Paper | Discussion Paper | Discussion Paper | Fig. 1. δ 18 O data for the EPICA (upper panel) and NGRIP (lower panel) datasets in the period 20 ÷ 43.5 kyr B.P.

Fig. 5 Fig. 6 .
Fig. 5. δ18 O data (black solid lines), sums of the j = 1, 2, 3 IMFs (red dashed lines) and of the j = 4, 5, 6 modes (blue dash-dotted lines) for the EPICA (upper panel) and NGRIP (lower panel) datasets.An offset was applied to the IMF sums to allow visualization in the same plot.

2 Datasets and analysis technique
The European Project for Ice Coring in Antarctica (EPICA) dataset provides a record of 100 yr averages of the oxygen isotope δ 18 O from the EPICA Dronning Maud Land Ice Core (75.00• S, 0.07 • E, 2892 metres a.s.l.) covering the period 10÷51 kyr B.P. (EPICAmembers, 2006).Data for the Northern Hemisphere come from the North Greenland Ice-Core Project (NGRIP, 75.10 • N, 42.32 • W, 2917 m elevation) and consist in 50 yr averages of δ 18 O extending back to 123 kyr B.P. (North Greenland Ice Core Project

Table 1 .
Characteristic periods of both EPICA and NGRIP datasets obtained through the EMD.Errors are estimated as standard deviations of the instantaneous periods.The period of the j = 6 mode for EPICA refers to the residue r m (t) (see text).
The local extrema are interpolated by mean 105 splines to obtain the envelopes of the maxima an Then the difference h 1 (t) between δ 18 O and m 1 ( ter being the mean between the envelopes of the m minima, is computed.If this quantity satisfies the two conditions, (i) the number of extrema and zero 110 does not differ by more than 1 and (ii) the mean the envelopes obtained from the local maxima an is zero at any point, it represents an IMF.If the above are not fulfilled, the described procedure i on the h 1 (t) time series, and h 11 (t) = h 1 (t) − m 11 115 m 11 (t) is the mean of the new envelopes, is calcul process is iterated until, after s times, h 1s (t) fulfil properties.To avoid a loss of information about and frequency modulations, a criterion to stop the scribed procedure has been proposed(Huang et al.