Objective identification of climate states from Greenland ice cores for the last glacial period

Introduction Conclusions References


Introduction
Measurements of δ 18 O for the last glacial period indicate large temperature changes occurring on timescales of centuries or less.They are most pronounced in Greenland ice cores where they are characterised by a rapid warming of up to 10 • C, followed by a slow cooling period.They are known as Dansgaard-Oeschger (DO) events and were numbered and dated by Dansgaard et al. (1993) using a timescale derived from a simple ice flow model and constrained by well established events at 11.5 kyr and 110 kyr BP (before present).There is evidence that the dramatic climate shifts at high latitudes in the Northern Hemisphere had global scale effects.For example, Blunier  2001) note general agreement of temperature changes between GISP2 and stalagmite records from the Hulu cave in China.Bond et al. (1993) used the abundance of planktic foraminifera in sediment records to show that sea surface temperatures in the NA are correlated with DO temperature cycles.They also note that there are longer periods of cooling, which each include several DO cycles.At the termination of these 10-15 kyr cycles there are signatures of large scale ice rafting events that correlate with a low temperature stadial.Known as Heinrich events, they are associated with collapse of the Laurentide ice sheet, thereby inundating the NA with freshwater, preventing the formation of North Atlantic Deep Water (NADW) and terminating the thermohaline circulation.The subsequent depletion of icebergs in the NA is thought to provoke a rapid reorganisation of the climate system and a temperature increase into a prominent interstadial followed by gradual cooling punctuated by a series of DO cycles.This does not explain the mechanism driving the shorter DO cycles, although they are considered to be linked with similar ice sheet calving affecting the salinity of the NA and the formation of NADW.
The apparent regular occurrence of DO events has lead some authors to suggest an external forcing or internal oscillation of the climate system.In particular the Fourier spectrum of the temperature record using the GISP2 time line has a significant peak at 1470 years (Grootes and Stuiver, 1997).However, the δ 18 O record from the Greenland Ice Core Project GRIP, completed the same year as GISP2, whilst also exhibiting DO cycles does not have this spectral peak.Neither does that from the North Greenland Ice Core Project (NGRIP) using the GICC05 timescale.Alley et al. (2001) suggest a simple model for the dynamics that includes weak periodic forcing and a stochastic resonance mechanism.Stocker and Johnsen (2003) link this to a bipolar thermal "see-saw", explaining the observed response in Antarctica and Braun et al. (2005) propose a resonance between two faster solar modes of variability to give a resulting weak signal at 1470 years.
It has been argued that Fourier methods are not suitable to detect recurrence patterns for such non-linear discrete events that are generated by a thresholding process.(2003) measures the deviation of the events from a perfect 1470 years signal and finds that their timing is well within the assumed dating error, although there are several 1470 year periods where no event occurs at all.Ditlevsen et al. (2007) assess the significance of the observed periodicity by defining null models and computing the Rayleigh measure and standard deviation of residuals for ensembles of surrogate data.

Rahmstorf
They find no evidence to reject the null model of randomly timed events.The exception is when using the GISP2 ice core record and omitting DO event 9 as numbered by Dansgaard et al. (1993).However, this does not change the result when using the more recent NGRIP ice core record.
The event detection algorithm of Rahmstorf (2003), which uses a criteria of 2‰ δ 18 O, does not include DO 9, potentially influencing the conclusion.Ditlevsen et al. (2005) defines a DO event as first up crossings of an upper level following a lower level, identifying several more events and concluding that they are randomly paced.
Given the sensitivity of results regarding the existence of periodicity to the definition of the events it is desirable to have a method that objectively detects their occurrence.
Previous work by Livina et al. (2010) used GRIP and NGRIP δ 18 O data to study the number of states in the climate for the last 60 kyrs using a polynomial fitting algorithm to windows of the data.They detect the two states corresponding to the stadial and interstadials of the last glacial and find that these merge to a single state around 25 kyrs BP.They detect more than two states for certain parts of the data although this seems to be dependent upon the size of the window used.Although their method allows for asymmetry between climate states this does not include the characteristic "saw tooth" features of Dansgaard-Oeschger cycles.
The aim of this paper is to study models that are able to capture the features of the cycles whilst also including the associated uncertainty in the identification of climate states.This uncertainty can then be included in subsequent analysis.We consider there to be hidden states corresponding to different climate regimes and identify these in a probabilistic sense.We apply a Bayesian statistical analysis to the three data sets NGRIP, GRIP and GISP2.We make no assumptions about the parameters (e. temperature and variance), instead these are computed from the posterior distribution using Bayes' rule.This is equivalent to a cluster analysis problem and will be discussed in Sect. 2. In Sect. 3 we analyse the results from the posterior simulation and find three distinct climate modes whilst in Sect. 4 we consider the periodic properties of these modes.

Hidden state models
We assume that there are M unobserved states governing the dynamics of the climate.In one instance we consider hidden state models where each state corresponds to a different temperature regime where, for example, one state could correspond to the high temperature interstadial and the other to the cooler stadial.This model is only concerned with absolute values of temperature.However, considering the "saw tooth" nature of the DO cycles it may be more appropriate to allow for an asymmetry between warming and cooling phases.Therefore, we also consider models for which different states correspond to different forcing regimes.In this case it is the increments of temperature that are of primary interest.For example, there may be one state for the sudden temperature increase in a DO event and another for the slow cooling in the remainder of the cycle.The parameters of the model include the mean temperature (or mean temperature increase), and the variance within each state.These are not fixed but are inferred directly from the data using Bayes' rule.
We analyse models with up to four hidden states.The occurrence of each state and all of the associated parameters are estimated in a Bayesian framework from the posterior distribution.In this way no a priori thresholds need be defined.We also determine the most suitable model for the data by estimating their marginal likelihoods.This involves approximating the integral over the whole parameter space to obtain an estimate of the probability of the model given the data.In this way, models with too many parameters, that over fit the data, are penalised.The ratio of marginal likelihoods between two models is known as the Bayes' factor and quantifies the evidence Introduction

Conclusions References
Tables Figures

Back Close
Full The data sets we study are shown in Fig. 1.The NGRIP and GRIP chronology are derived from the ss09sea timescale whilst GISP2 is dated according to the Meese/Sowers timescale.Note the good agreement between the GRIP and NGRIP timescales isotope records, whereas GISP2 is similar only for the last 65 kyrs or less.
The discrepancy between GISP2 and the other two records is due to the differing dating methods.
For each data set we estimate the probability that it was generated by a given model.This will provide a reliable way to assess which model is most appropriate as a physical description of the climate, not being overly dependent on a single ice core.
For the increments model we assume that the time series can be modelled as a random walk dependent upon the state of the climate.This is intended to capture the differing distributions of the increments that occur during different climate regimes.Note that one could use a more general first order process, however we found that this did not significantly affect the identification of climate states.
Under the random walk model the data X i at time i evolves according to where i ∼ N (0,1) is the standard normal distribution, µ S i is the forcing and σ S i is the standard deviation in state S i .For N data points, the probability of the data given the unobserved state sequence is Analogously the probability of the data for the absolute value model is Full The following details are for the incremental model, where S j ∈ {0...M}; the absolute value model is simpler.Readers familiar with Bayesian inference can skip to Sect. 3.

Bayesian inference for hidden state models
To capture our uncertainty about the time spent in each climate regime an extra layer is introduced into the hierarchy that conditions S on an unknown parameter λ j , j ∈ {0...M}.This is the probability that S i is in state j , which is then distributed according to the multinomial distribution.In a Bayesian setting λ j , and the other parameters, are random variables with a prior distribution.By defining conjugate priors the model can be inferred using a Gibbs sampling algorithm (for an introduction to Bayesian inference using Markov Chain Monte Carlo (MCMC) methods see Gilks et al., 1995).The conjugate prior for λ j is given by a Dirichlet distribution with hyperparameter α.This respects the constraint j λ j = 1.For σ j one defines τ j = 1/σ 2 j and assigns a Gamma prior τ j ∼ Γ(a,b), whereas µ j is a priori normally distributed µ j ∼ N (0,σ 2 µ ).The hyperparameters a, b and σ µ are fixed.The posterior distribution can then be simulated by drawing sequentially from each of the conditionals.For µ j , we have where I[S i = j ] = 1 if S i = j or 0 otherwise.The posterior for τ j is the Gamma distribution Each unobserved state is sampled from the updated multinomial distribution with probability Tables Figures

Back Close
Full By cycling through the above distributions one obtains a dependent sample from the posterior.As is usual for MCMC algorithms we discard an initial "burn in" period to allow for convergence to the posterior distribution.
To estimate the probability of the model given the data we compute the marginal likelihood using the Laplace-Metropolis method as in Raftery (1996).For data X , parameter vector θ and model M, the marginal likelihood is As the hidden states are conditionally independent given θ the above sum can be performed analytically.Then the integral is estimated using samples from the posterior with the Laplace approximation where θ is the posterior mode and Ψ is the covariance matrix computed from the MCMC output.

Estimation of climate states
We first determine the model most supported by the data by estimating the marginal likelihoods using the method described above.The logarithm of the marginal likelihoods, shown in Table 1, were computed using 20 000 samples after a burn in period of 10 000.Note that the log marginal likelihoods are higher for GISP2 because of the smoothing required to produce 50 year running means.For all data sets there is strong evidence supporting the incremental model over the absolute value model.The log Bayes Factor, of the ratio between the two model types is typically of order 10 3 .This is decisive according to the guide scale in Kass and Raftery (1995).Table 1a implies that, for NGRIP and GISP2, the three state model is favoured.For GRIP the Bayes Factor 1216 Introduction

Conclusions References
Tables Figures

Back Close
Full is 1.2, which according to (Kass and Raftery, 1995) does not favour either model.This implies that, at least for NGRIP and GISP2, there is evidence to support three dynamic climate modes.The temperature decreases at the termination of the DO cycles have a significantly larger cooling effect than the gradual cooling state due to steady ice berg flux.We shall refer to these as termination events.Note that no more than three distinct states can be identified from any of the data sets implying that the increments model with three climate regimes is the best description of the dynamics.
Figure 2a shows the posterior probability of a DO event for the NGRIP data, i.e. for each time i the probability that S i = 1 (events with probability less than 0.5 omitted for clarity).The posterior distribution of hidden states is similar for the other data sets, although more events are detected in the GISP2 data (not shown).The canonical DO events of Dansgaard et al. (1993), indicated by a blue marker, all have posterior probability greater than 0.8, implying that the model is able to distinguish these events from the other states.There are, however several other events detected at the 0.8 probability level.Although not typically regarded as DO events in the literature, these fluctuations share similar statistics and may be important in the analysis of recurrence patterns.Many of the high probability regions cluster together and would previously have been classed as a single event.Here, we only regard two events as one if they occur in immediate succession.In this way we make little restriction on the minimum time period between events.
Figure 2b shows the equivalent for the termination events.As mentioned in the introduction, these may be associated with large ice rafting.The regions of high posterior probability do not clearly occur at the end of all DO cycles but, as discussed above there is statistical evidence that this third state should be included.
The mean waiting times between DO and termination events can be computed directly if we assume a random Poisson model where the mean waiting time T j is related to the probability λ j by λ j = 1−exp(−50/T j ).This is the probability of at least one event occurring within the 50 year window.Figure 3  waiting time of about 1400 years.However, the smaller values of GISP2 imply that it is attributing more data to the DO and termination events.The posterior distributions for the parameter of interest µ are shown in Fig. 4. Note that this parameter has units of δ 18 O‰ and can be considered to be proportional to the average temperature change within a 50 year period for each climate state.There is a clear distinction in this forcing between different climate states.µ 0 , corresponding to the relaxation state, has a mean of approximately −0.05 for NGRIP and GRIP and −0.03 for GISP2.The DO state is larger as expected and there is also good agreement for this between NGRIP and GRIP with a mean of approximately 1.90-2.00.Finally the termination state has mean −1.90, except for GISP2 with −0.60.In both cases GISP2 has smaller absolute magnitude forcings, most likely due to the smoothness of this time series.Note that the mean of 2.0‰ for the DO events is the same as the detection criteria of Rahmstorf (2003).Note that this is a result of our analysis and has not been a priori imposed.

Analysis of periodicity
As in Ditlevsen et al. (2007) we use the Rayleigh measure to quantify the amount of periodicity for discrete events.This is a measure for the average phase coherence for N events; for period T it is given by R For each period T we compute the Rayleigh measure and assess the significance of the result by comparing it to the distribution generated from a null ensemble of the original data.
The null ensemble is generated by sampling the data without replacement so that it is randomly ordered.Then hidden state sequences are generated by sampling from the posterior distribution conditioned upon the randomly ordered data and the other parameters (inferred for the original data).This means that the climate states in the ensemble will be equal in distribution to those in the actual data but their occurrence will be randomised.Introduction

Conclusions References
Tables Figures

Back Close
Full Figure 5 shows the Rayleigh measure over a range of periods T with the mean and 2σ values for the null ensemble.NGRIP and GRIP are within the 2σ value for all T , however, the value of GISP2 is significant for a period equal to 1450 years, in agreement with previous studies.

Conclusions
We propose statistical models able to identify different climate regimes from temperature proxy records with no threshold or minimum duration criteria.Using a Bayesian model comparison method we find that models based on the temperature dynamics, rather than absolute values, are better supported by the data.This leads one to consider how appropriate it is to define different climate regimes, for this region, based on absolute temperature differences, i.e. the stadial and interstadial states.Rather, we find it more useful to work in terms of the different forcing regimes, namely slow cooling due to gradual freshwater input to the North Atlantic, ice sheet collapse and then sudden temperature increase as the thermohaline circulation recovers.These states have distinct posterior distributions for the rate of temperature change within each climate state.We found that the temperature series is best characterised by this three state model with distinct forcing regimes.
Using this model we investigated the recurrence properties of the events.We found a significant signal at the 1450 year period for GISP2, agreeing with previous authors and supporting our method.This period is absent in the other data sets, in agreement with Ditlevsen et al. (2005).Due to the strong agreement between the climate regimes of GRIP and NGRIP, and their recurrence times, we assign a higher significance to the results from these data sets.We conclude that there is no lasting periodic signal for the rapid warming or cooling events of the last glacial period and suppose that the apparent periodicity seen in GISP2 is due to a small number of events that occur at the regular interval of 1450-1500 years.In future work we could incorporate an extra layer into these models, conditioning on dating uncertainty data.This would provide Introduction

Conclusions References
Tables Figures

Back Close
Full et al. (1998) observed an out of phase signal in Antarctic ice cores, whereas Wang Discussion Paper | Discussion Paper | Discussion Paper | et al. ( Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | supporting one model in favour of the other.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | displays the posterior waiting time distributions.There is good agreement between NGRIP and GRIP, both centred on a mean Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | a more robust estimate of the strength of any periodicities.Uncertainty estimates are now available for NGRIP 20 year mean averages.Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 1 .Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 1.Fifty year means for stable isotopic ratios (δ 18 O) from NGRIP, GRIP and GISP2 expressed in ‰ with respect to Vienna Standard Mean Ocean Water.Note that irregular spaced values in the GISP2 record have been interpolated to give a 50 year sample rate.

Table 1 .
Estimates of the marginal likelihood for each model with 1-3 hidden components for each data set.