Past Patterns of millennial variability over the last 500 ka

. Millennial variability is a robust feature of many paleoclimate records, at least throughout the last several glacial cycles. Here we use the mean signal from Antarctic climate events 1 to 4 to probe the EPICA Dome C temperature proxy reconstruction through the last 500 ka for similar millennial-scale events. We ﬁnd that clusters of millennial events occurred in a regular fashion over half of the time during this with a mean recurrence interval of 21 kyr. We ﬁnd that there is no consistent link between ice-rafted debris deposition and millennial variability. Instead we speculate that changes in the zonality of atmospheric circulation over the North Atlantic form a viable alternative to freshwater release from icebergs as a trigger for millennial variability. We suggest that millennial changes in the zonality of atmospheric circulation over the North Atlantic are linked to precession via sea-ice feedbacks and that this relationship is modiﬁed by the presence of the large, Northern Hemisphere ice sheets during glacial periods.


Introduction
(Background) Ice-core proxy records of high latitude Northern Hemisphere temperature reveal a distinctive pattern of repeated rapid warming events of 8-15 • C during Marine Isotope Stage 3 (MIS 3), known as Dansgaard-Oeschger (D-O) events (Dansgaard et al., , 1993Oeschger et al., 1984). These rapid warmings are interspersed with cold periods such that MIS 3 is a period of substantial millennialscale climate variability. This variability is found throughout much of the Northern Hemisphere in marine sediments and also continental records (Shackleton et al., 2000;Wang et al., 2001;Voelker et al., 2002;Rohling et al., 2003;Denton Correspondence to: M. Siddall mark.siddall@bristol.ac.uk et al., 2005). D-O events in MIS 3 appear in groups with a longer warm period followed by several shorter warm periods, interspersed with cold periods (Bond and Lotti, 1995). These D-O groups, sometimes called Bond cycles, end in a cold culmination, expressed as a so-called Heinrich event of massive ice-rafted debris (IRD) deposition in the North Atlantic between about 40 and 50 • N (see overview in Hemming, 2004). Such indications of substantial iceberg release have motivated modelling studies which consider freshwater pulses into the North Atlantic as a trigger for rapid cooling during D-O variability (Stocker and Wright, 1992). Blunier et al. (1998) and Blunier and Brook (2001) used variations in the concentration of atmospheric methane (a globally well-mixed gas) in air bubbles enclosed within the ice, to synchronise ice-core records from Antarctica and Greenland. This work showed that D-O events correspond to slower, smaller changes in Antarctica. Longer D-O stadials in Greenland correspond to the largest variations in Antarctic temperature (Antarctic events A1 to A4) (Stocker and Johnsen, 2003). High-resolution temperature reconstructions from Antarctic ice cores have confirmed that shorter D-O events are also associated with smaller Antarctic warming events and that the magnitude of the Antarctic warming is proportional to the duration of D-O stadial events (EPICA Community members, 2006). This paper is concerned with the larger Antarctic climate events, corresponding to the longer D-O stadials during MIS 3.
The detailed phase response between Antarctic and Greenland records during MIS 3 has been described by Stocker and Johnsen (2003) -in response to rapid warmings in Greenland, Antarctica starts to gradually cool and in response to rapid coolings in Greenland, Antarctica starts to gradually warm. The physical mechanism responsible for this "bipolar seesaw" is linked to an active meridional overturning circulation (MOC) in the Atlantic, which transports heat northward at all latitudes, heating the north and cooling the south (Crowley, 1992 drawn to the north from the south along the axis of the Atlantic, causing rapid cooling in Greenland and gradual warming in Antarctica. In this way it is suggested that cooling during D-O stadials, plausibly caused by iceberg release also affects Antarctic temperature (Hemming, 2004). Here we consider the link between IRD deposition in the north Atlantic and the larger Antarctic climate events over the last 500 kyr.
Although millennial-scale variability has been predominantly found in high-resolution records of MIS 3, some evidence for it exists also in earlier intervals (Oppo et al., 1998;McManus et al., 1999). Larrasoaña et al. (2003) found that the occurrence of sub-Milankovitch (or sub-orbital) variability became significant from as early as ∼0.95 Ma. A study of millennial events recorded in the Vostok ice core has suggested that broadly similar north-south phase relationships to those of MIS 3, as described in Stocker and Johnsen (2003), existed over the last four climatic cycles (Delmotte et al., 2004). Siddall et al. (2006Siddall et al. ( , 2007 described very similar millennial-scale relationships between temperature reconstructions in the Vostok ice core, methane concentrations and a simple model of the bipolar seesaw during MIS 8, compared to MIS 3. Cave speleothem records have revealed millennial-scale variability within the last several glacial cycles, which is similar in appearance to that of MIS 3 ( Bar-Matthews et al., 2003;Yuan et al., 2004).
Several authors have suggested that millennial-scale climate variability is a function of ice volume, relating periods of more intense variability to periods with increases ice volume (Oppo et al., 1998;McManus et al., 1999;Siddall et al., 2006Siddall et al., , 2007. Other work has found millennial variability during periods with reduced ice volume (Martrat et al., 2007;Desprat et al., 2009). Most recently Ganopolski et al. (2009) have linked episodes of millennial variability to the complex phase relationships between Antarctic temperature and atmospheric CO 2 concentrations during the terminations, requiring millennial variability during periods of reduced ice volume.
Ice core records from Antarctica extend down to 800 ka BP (EPICA Community Members 2004), although adequate resolution to consider millennial changes only exists for the last 500 kyr. New, multi-century resolution sea-level data from the Red Sea extend down to 500 ka BP (Rohling et al., 2009). Combined with IRD data from the north Atlantic (Oppo et al., 1998;McManus et al., 1999), these data allow a systematic analysis of millennial climate variability over the last 500 kyr (Fig. 1). Here we take a simple approach to analysing millennial variability over the last 500 kyr using these records. Although 800 kyr of data is available from the EPICA Dome C (EDC) ice core, we consider only the last 500 kyr for several reasons: (i) Red Sea sea-level and north Atlantic IRD data only extend to the last 500 kyr and we compare our results to these and; (ii) the resolution of the EDC deuterium record deteriorates with depth and millennial comparison is not possible much beyond 500 kyr ago.

Method
Our method is simple but effective. Essentially we take a representative A-event (in this case the mean of A-events 1 to 4) from the MIS 3 section from the EDC ice core record as recorded in the deuterium temperature proxy, and look for similar events in the rest of the EDC record. We further validate our results by carrying out a series of Fast Fourier Transforms with overlapping windows (known as a spectrogram analysis) in addition. Our approach consists of five steps: 1. Removing the long-term trend from the EDC record We first treat the last 500 kyr of EDC deuterium data on the EDC3 time scale to remove the 6 kyr trend (Fig. 1a). We achieve this by carrying out least-square linear regressions on 6 kyr windows centered on each data point in turn and removing the temperature trend. We will call the detrended, evenly sampled EDC deuterium record AA detrended (Fig. 1b). The window length of 6 kyr is chosen to effectively remove the trend during the glacial terminations. Sensitivity tests varying the length of this 6 kyr window by a factor of two (between 3 ky and 12 kyr) show no impact on our results.

Defining the probe signal
We define an average A-event and use this event to probe the rest of the signal to look for similar events. To define an average A-event we took 4 kyr windows centered around 37.9, 46.2, 52 and 58.4 kyr ago representing A-events 1 to 4, respectively (Fig. 2). The four windows were overlain around their centres and the probe signal constructed by binning the combined data using the mean over 0.1 kyr intervals. Although, in the example given here we have used an average A-event, we could have used any of the four events individually because our result is not sensitive to which of the MIS 3 A-events we use. We will call the probe signal A mean .

Probing AA detrended with A mean
We compare AA detrended with A mean by calculating the standard difference between AA detrended and A mean for overlapping windows separated by 0.1 kyr The standard difference is then normalised to its greatest value over the last 500 kyr. We call this result AA probed (Fig. 1c).
4. Calculate the standard deviation of AA probed over 8 kyr windows Finally, because AA probed is sensitive to both the presence of, and the phase of, millennial variability similar to that of A2, it exhibits both low (in phase) and high (out of phase) values during periods of millennial variability. Contrastingly, AA probed shows little variation during periods without millennial variability. Therefore we also calculate the standard deviation from the mean of AA probed over 8 kyr windows overlapped by 0.1 kyr.  (Oppo et al., 1998;McManus et al., 1999). Light grey bars represent periods of sustained ice sheet growth and the absence of millennial variability and dark grey bars represent periods with ice volume greater than 80 m sea-level equivalent. This result is normalised by the maximum standard deviation over the 500 kyr study period. The output of this calculation (Fig. 1d) is simple to interpret because where the standard deviation of AA probed is high, the record contains millennial signals similar to A mean .

Validation using a spectrogram analysis
Because of the novelty of our approach, we have also carried out a more routine spectral analysis of AA smoothed using FFTs over 30 kyr windows with a 29.999 kyr overlap standard. The resulting spectrogram is shown in Fig. 3a, where the millennial band (3 to 5 kyr period) is highlighted. By taking the mean of the spectrogram within this millennial window of 3 to 5 kyr period, we produce a time series of the mean spectral power in this band. This is compared to precession forcing in Fi. 3b and the running standard deviation of AA probed in Fig. 3c. We note that the two approaches in Fig. 3c give a very similar result, confirming that our result is not dependent on our novel approach.

Results
The results of our analysis are given in Fig. 1, where they are also compared to the Red Sea sea-level reconstruction and the IRD record of ODP 980 from the north Atlantic. The results are repeated in Fig. 3, where they are compared to our spectrogram results.
The most obvious result is that millennial variability similar to that of MIS 3 is repeated throughout the last 500 kyr in the EDC deuterium record (high values in Figs. 1d and 3 and Table 1). Periods of increased variability occur with striking regularity throughout the last 500 kyr. Far from being unusual events, they intermittently occupy approximately half of the record of the last 500 kyr.
In order to compare our results with reconstructions of ice volume, we plot the Red Sea sea-level curve (Fig. 1e). Clustering of millennial-scale variability is found to occur during two distinct phases within the glacial cycles, namely during periods of intermediate ice volume and extensive IRD deposition and during glacial terminations (including the early parts of interglacials) (Fig. 1e, f and Table 1). Conversely, millennial-scale variability also appears to be absent during two distinct phases within the glacial cycles, namely during periods of extensive ice-sheet growth (marked with light grey bars in Fig. 1 and referenced to periods of sea-level fall in Fig. 1e), and during glacial maxima with more than 80 m of sea-level equivalent ice volume (marked with dark grey bars in Fig. 1 and noted in Table 1).
Over the last 500 kyr, we observe twenty-four clusters of increased/decreased millennial-scale variability. These clusters are typically separated according to a "cycle" of ∼21 kyr, a value similar to that of precession. However, the relationship between precession and millennial variability is not simple. Precession is generally in phase with the presence of millennial variability during periods of intermediate ice volume (40 to 80 m sea-level equivalent) and in antiphase during periods of reduced ice volume (40 m or less sea-level equivalent) (Figs. 1d, e, f, 3 and Table 1). Between periods when clusters of millennial variability are in and out of phase with precession, there are transition intervals when precession leads the millennial variability.

Discussion
Previous work has pointed out the links between millennialscale variability and intermediate ice volume of around 40 to 80 m ice volume equivalent (e.g. Oppo et al., 1998;Mc-Manus et al., 1999;Siddall et al., 2007). Other work has suggested that the link between millennial-scale variability and ice volume may be less important (Martrat et al., 2007;Desprat et al., 2009). Yet other work has focused on the role of millennial-scale variability in glacial terminations (Ganopolski et al., 2009).
Our signal analysis approach finds that there are two distinct regimes during which millennial-scale variability develops. One is associated with intermediate ice volume and episodes of IRD deposition (in line with the conclusions of Oppo et al., 1998;McManus et al., 1999, andSiddall et al., 2007), and the other is associated with glacial terminations (in line with the conclusions of Ganopolski et al., 2009). The The natural log of the mean power in the spectrogram with a period between 3 and 5 kyr. Precession forcing from Berger and Loutre (1991) is shown for comparison. (C) The standard deviation of AA probed over 8 kyr windows overlapped by 0.1 kyr and the natural log of the mean power in the spectrogram with a period between 4 and 5 kyr. Light grey bars represent periods of sustained ice sheet growth and the absence of millennial variability and dark grey bars represent periods with ice volume greater than 80 m sea-level equivalent (see Fig. 1e).
Millennial variability: H=high, L=low Phase with respect to precession: O=out of phase, T=transition, I=in phase Sea level: H=high, L=low, D=deglaciation, F=falling IRD: Y=yes, N=no similarity between both classes of millennial-scale events in the EPICA Dome C deuterium records suggests that both classes reflect the same processes.
If millennial-scale variability exists for approximately half of the time during the last 500 kyr, then there is an equally important reverse question: "What characterises periods which have reduced millennial variability or for which millennial variability is absent?" In Fig. 1 these periods are characterized by either ice volumes greater than ∼80 m icevolume equivalent or periods of sustained ice sheet growth. Many authors have shown evidence for millennial-scale sealevel fluctuations associated with millennial-scale climate variability during MIS 3 (Chappell 2002;Cutler et al., 2003;Siddall et al., 2003Siddall et al., , 2008Rohling et al., 2008;Arz et al., 2007). We suggest that sustained ice-sheet growth requires the absence of millennial variability to allow the continuous growth of ice sheets, rather than the growth and collapse of ice sheets during periods such as MIS 3. The lack of millennial variability at the glacial maxima may be evidence that the presence of the large ice sheets issues in a relatively stable climate during those periods.
Given that only one of our classes of millennial-scale event is linked with IRD deposition, we infer that massive iceberg calving may not be the ultimate cause of millennial variability. We note that IRD events only last several centuries (less than D-O stadials, which last 1-2 kyr) and peak at the conclusion of D-O stadials (Hemming, 2004). This suggests that IRD events are caused by the D-O stadials, rather than being the cause of D-O stadials. Hence, IRD deposition (iceberg calving) does not appear to be the driving mechanism, but merely a feature of the class of millennial-scale variability that occurs during periods of intermediate ice volume.
The above implies that another causative mechanism is required to explain D-O variability. Because the recurrence interval of clusters of millennial variability is close to precession any discussion of processes explaining millennial vari-ability needs to consider precession forcing. It has been suggested that changes in the zonality of atmospheric circulation over the North Atlantic might instead be the cause of millennial variability by switching the latitude of the polar front between different modes (Seager et al., 2002;Seager and Battisti, 2007). As noted above, the mean recurrenceperiod of clusters of millennial-scale variability in Fig. 1c,d is 21 kyr, similar to precession. This implies that changes in the zonality of atmospheric circulation over the North Atlantic may be linked to precessional cycles. In turn this implies zonal teleconnections in the Northern Hemisphere. Such links are indeed suggested by Wang et al. (2008) who found both strong precessional and millennial variability in Chinese speleothem records.
If precession forcing sets the atmospheric preconditions for millennial variability, what mechanisms might link precession to individual millennial events? Several possible links have been suggested linking the atmospheric circulation in the Northern Hemisphere to millennial changes in the MOC. Changes in sea ice extent have been suggested to generate seesaw-type temperature variations in the MOC (Kaspi et al., 2004) and are linked to variations in atmospheric circulation in the Northern Hemisphere (Li et al., 2005). Sea-ice feedbacks with the MOC might therefore be the mechanism behind individual millennial events, given the correct preconditions with respect to precession forcing. A second explanation might be given by the direct link between the magnitude of the wind stress in the North Atlantic and the stability of the MOC found by Ashkenazy and Tziperman (2007). Such an explanation would remove the requirement of seaice feedbacks to force changes in the MOC. Finally Friedrich et al. (2009) have used the LOVECLIM model to find a feedback between sea-ice anomalies and the atmospheric circulation. This feedback in turn drives atmospheric anomalies over Hudson Bay which form a freshwater pulse into the Labrador Sea, which then causes the MOC to collapse. In order to consider the possible link between precession forcing and the occurrence of millennial variability, we compare our result with precession forcing in Fig. 1d. Clusters of millennial variability are in phase with precession during periods of intermediate ice volume but alter phase during glacial maxima (when there is little or no millennial variability). Clusters of millennial variability then remain in antiphase during the glacial termination and periods of reduced ice volume. During the transition from reduced to intermediate ice volume the phase realigns again to put the reoccurrence of millennial variability back in phase with precession forcing. Marked changes to the atmospheric circulation due to the presence of large ice sheets and extended sea ice at the Last Glacial Maximum (LGM) have been found in atmospheric models (e.g. Li and Battisti, 2008). Otto-Bliesner et al. (2006) also found that atmospheric circulation was sensitive to the detailed configuration of the north American ice sheets at the LGM.
In summary, we speculate that the existence of ice sheet topography and albedo effects modifies the influence of precession forcing on millennial variability and that this may explain the existence of two classes of millennial variability. One class of millennial variability is linked to precession forcing, the other class requires precession forcing modified by the existence of ice sheets of intermediate size. Because precession forcing affects how much insolation an entire hemisphere receives in a certain season, precession may be an important driver of changes in seasonal sea ice. We therefore suggest that sea-ice feedbacks are the most likely candidate for linking millennial variability to precession forcing.

Conclusions
It is clear from our analysis that millennial-scale variability is a much more common feature of Quaternary climate than has hitherto been realised, taking up a significant proportion of the variability during periods of intermediate icesheet growth and glacial terminations.
Millennial variability occurs during distinct periods in the glacial cycles, namely periods with ice volume equivalent to 40-80 m of sea-level lowering, and glacial terminations.
Different periods give rise to the absence of millennial variability, namely periods of extensive ice-sheet growth, and periods with ice volume equivalent to more than 80 m of sealevel lowering.
The lack of a unique relationship between ice volume and millennial-scale variability calls for a causative mechanism other than iceberg release into the north Atlantic. We suggest that changes in the zonality of atmospheric circulation over the North Atlantic may be a viable alternative (Seager et al., 2002;Seager and Battisti, 2007). Because we find the mean recurrence-period of the millennial-scale variability clusters to be ∼21 kyr, we suggest that changes in the zonality of atmospheric circulation over the North Atlantic may be linked to precession. Because the phase relationship between precession forcing and clusters of millennial variability changes over the glacial cycle, we suggest that the existence of ice sheets modifies the influence of precession forcing on millennial variability. This suggestions is supported by Chinese speleothem records, which show strong precessional and millennial variability over the last two glacial cycles, indicating important changes in the zonality of atmospheric circulation on these time scales (Wang et al., 2008).
If millennial variability requires preconditions with respect to precessional forcing, it is unlikely to occur due to natural forcing for some thousands of years from now.