Recent analysis of 38 globally distributed paleoclimatic records covering Marine Isotope Stage 3 (MIS 3) 60–26 ka demonstrated that the two leading empirical orthogonal functions (EOFs) explaining the data are the Greenland ice-core signal (“northern” signal) and the Antarctic ice-core signal (“southern” signal). Here singular spectral analysis (SSA) is used to show that millennial-scale variability of each of these two leading EOFs is characterized by two independent modes. The two modes of each EOF share similar relative distributions of variance, identical spectra, and, where each mode has spectral power, coherency spectra, which are significantly above the null hypothesis level at 95% confidence. The only difference between the modes of the northern and southern signals is that they are phase shifted. The phasing and long response time of the low-frequency mode, combined with its relationship to atmospheric CO2 and sea level, are consistent with coupled changes in the ocean, ice sheets, atmosphere, and carbon cycle, whereas the phasing and short response time of the high-frequency mode are consistent with an atmospheric transmission likely induced by changes in hemispheric sea ice distributions and attendant feedbacks.
Millennial-scale climate variability is a persistent feature in the spectrum of climate change (Mitchell 1976). Among the best-known expressions of this variability are the Dansgaard–Oeschger (D-O) Oscillations recorded in Greenland ice cores during Marine Isotope Stage 3 (MIS 3) 60–25 ka (Grootes et al. 1993). MIS 3 millennial-scale variability is also recorded in Antarctic ice cores, but with a fundamentally different structure (so-called A events) than seen in Greenland records (Barbante et al. 2006; Blunier and Brook 2001; Blunier et al. 1998). A key challenge in understanding the origin of these millennial-scale climate changes and their possible interhemispheric linkages involves distinguishing regional from global variability (Roe and Steig 2004; Stocker and Johnsen 2003; Wunsch 2003).
2. Data analysis
Previous statistical analyses of millennial-scale change focused on just the Greenland Ice Sheet Project 2 (GISP2) and Antarctic (Byrd) ice-core records (Roe and Steig 2004; Wunsch 2003). We base our results on analysis of a set of 38 published proxy records that significantly expand the global coverage (Clark et al. 2007; see supplementary Table 1 http://dx.doi.org/10.1175/2010JCLI3416.s1). Thirty-one of the 38 proxy records used in the original empirical orthogonal function (EOF) analysis (Clark et al. 2007) were based at least in part on the GISP2 time scale (Meese et al. 1997), either through synchronization of ice-core records with methane or through establishing tie points associated with particularly prominent climate events for that part of the record beyond radiocarbon dating control. Other records used in the analysis, however, were based on independently derived age models, some of which were known to differ from the GISP2 time scale, thus introducing some uncertainty in the EOF results. We have now transferred those age models based on the GISP2 time scale to the Greenland Ice Core Chronology 2005 (GICC05) that is based on annual layer counting back to 60 ka and is in good agreement with a number of independently dated records and reference horizons (Svensson et al. 2008). The changes in this new chronology are insignificant for the record less than 40 ka but are as much as 2500 yr older for the period from 40 to 60 ka. Our goal was to place all records that had been correlated to the GISP2 record into a common chronostratigraphic framework.
The revised EOF analysis remains largely unchanged, suggesting the result is robust to chronological errors; the two leading EOFs explain about 48% of the total data variability (global signals), with the remaining 52% explained by many more EOFs (regional signals). EOF1 (30% of the data variance) is essentially the Greenland ice-core signal of D-O variability (Fig. 1a), and is referred to as the “northern” signal based on its spatial distribution largely in the Northern Hemisphere, whereas EOF2 (18% of the data variance) is essentially the Antarctic ice-core signal of A events (Fig. 1b), and is referred to as the “southern” signal based on its spatial distribution largely in the Southern Hemisphere. Many of the climate records analyzed, however, are best described as a combination of the two primary EOF signals, indicating a complex pattern of variability (Clark et al. 2007).
We have further evaluated our results with respect to chronologic errors found in the MIS 3 datasets used in these analyses. In particular, the primary climate signals obtained from the EOF analysis remain unchanged when random errors in chronologies are imposed on the MIS 3 dataset. In these experiments we inserted random errors in the chronologies into each record used in the EOF analysis. The mean of the chronology errors inserted into any one record ranged from 250 to 450 yr. These mean errors are comparable to the mean differences between many of the chronologies assigned to GISP ice-core records or to chronologies adjusted to better fit a proxy to match the GISP ice core. The EOF analysis tends to average out these chronological errors and the extracted EOF signals have chronologies that are close to the chronologies of the EOF based on the unperturbed chronologies.
We note that in EOF analysis, the spatial patterns of the EOF are the eigenvectors of the data correlation–covariance matrix in which the vectors are constrained to be orthogonal (the vector dot products between EOF patterns is 0). In the case of the EOF time series, however, there is no constraint on the statistical relationships at the nonzero lag (i.e., if phase shifts are considered). Accordingly, we use singular spectral analysis (SSA; Ghil et al. 2002) to further examine the temporal relationships between the two primary signals identified by the EOF analysis (Clark et al. 2007) and evaluate whether they are responding to a common global forcing. SSA extracts a set of temporal modes of variability that describes some fraction of the original time series (see the supplementary information). Unlike many other time series techniques, however, SSA analysis does not require the assumption that the data are periodic or stationary.
We find that four SSA modes explain 87% of the variance in the northern signal (SSA1N–SSA4N) and four SSA modes explain 95% of the variance in the southern signal (SSA1S–SSA4S; Fig. 1), with a similar relative distribution of variance between the two sets of four modes (34.4%, 33.9%, 12.1%, and 6.8% in SSA1N–SSA4N; 41.1%, 37.8%, 11.8%, and 4.4% in SSA1S–SSA4S). The maximum in the frequency spectra of SSA2N and SSA1S is at 6.8 kyr and the maximum in the spectra of SSA1N and SSA2S is at 4.3 kyr (Figs. 2a,b; see supplementary information for specifics on spectral and cross-spectral analysis). We note that this apparent reversal in the primary two modes of each signal is likely insignificant because the relative distribution of variance in SSA1N and SSA2N is nearly indistinguishable (34.4% versus 33.9%). The spectra of the remaining two modes identified in the northern and southern signals are, within statistical error, nearly identical (the maximum is 2.6 kyr for the third SSA mode, and 1.6 kyr for the fourth SSA mode; Figs. 2c,d). In all cases, the cross-spectra results show that not only are the spectra for each pair essentially identical where each time series has significant variance, but also that where each mode has spectral power, the coherency spectra are all significantly above the null hypothesis level at 95% (Figs. 2e–h). The only difference between the modes of the northern and southern signals is that they have a significant phase shift (Figs. 2i–l). Further cross-spectral analysis, however, shows that the two intermediate modes have shared variance with the lowest- and highest-frequency modes, indicating that only the lowest- and highest-frequency modes are independent of each other. We thus conclude that the global signal of millennial-scale variability is composed of just two independent modes (SSA1S or SSA2N, and SSA4), with hemispheric differences within each mode reflecting the response time of the processes involved in transmission of the signals, which manifests itself in the phase shift shown in Fig. 2.
3. Discussion and conclusions
What might these modes of climate variability represent? In Fig. 3 we compare the two SSAN modes with the GISP2 δ18O record and the two SSAS modes with the Byrd δ18O record. The dominant low-frequency signals in these records are found in each of the first SSA modes, with a clear association between these modes and Heinrich events (Figs. 3a,b). Moreover, the low-frequency mode accounts for the majority (∼35%–40%) of the variance in both the northern and southern signals, indicating that the process responsible for this mode is the dominant one in millennial-scale global climate change (Roe and Steig 2004).
To gain insights into the origin of the low-frequency mode, we calculated the squared coherence and phase spectra of the original dataset with SSA2N. Only 6 of the 38 sites have no significant coherence, with the squared coherence among the remaining sites being high (>0.70) throughout the mid- to low latitudes (Fig. 4a). The calculated phases with respect to SSA2N show a complex spatial distribution, but the majority of sites fall into one of two patterns (Fig. 4b). Given the sign convention of SSA2N (high value = cold climate), most sites in the tropics and northern latitudes show an 180° out-of-phase relationship [dark blue and red sites, Fig. 4b, i.e., cold (warm) SSA2N coincides with cold (warm) proxy]. In contrast, most sites from the Southern Hemisphere are slightly phase shifted behind SSA2N with a phase relationship indicating a lagged bipolar seesaw like pattern between SSA2N and Southern Hemisphere proxies (e.g., warming in SSA2N leads cooling in Southern Hemisphere proxies; Fig. 4b).
Another indication of the origin of the low-frequency mode comes from looking at the phases between the northern and southern SSAs. If we assume that the northern modes are forcing the southern modes through a simple linear system, the response time can be calculated from the observed phase angle. At the dominant low-frequency mode this calculated response time is about 3.5 kyr (Table 1), which is considerably longer than the mixing time of the ocean, but is consistent with a significant role by ice sheets.
We note that the maximum in the spectra of the first mode (6.8 kyr) is similar to the maximum in the spectra of millennial-scale variability in atmospheric CO2 (Ahn and Brook 2008) and sea level (Clark et al. 2007; Siddall et al. 2003; Fig. 5), which are not included in our analysis. Cross-spectral analysis clarifies this relationship in showing that SSA2N is coherent and in phase with CO2 and SSA1S is coherent and in phase with sea level (Fig. 5).
The uniform importance of this low-frequency mode in the Northern and Southern Hemispheres, its phase relationship, its relatively long response time, and its relationship to CO2 and sea level are consistent with strongly coupled changes in the ocean, atmosphere, ice sheets, and carbon cycle as simulated in some models (Clark et al. 2007). Each cycle starts with an abrupt warming in the North Atlantic sector associated with resumption of the Atlantic meridional overturning circulation (AMOC) and reduced sea ice extent, with attendant climate responses reflected by the northern signal (SSA2N). Our finding that changes in CO2 are in phase with the northern signal is consistent with models that simulate coupled changes in the AMOC and the oceanic biological pump through the former’s control on the global nutrient inventory (Schmittner and Galbraith 2008).
The enhanced cross-equatorial heat transport in the Atlantic basin leads to a cooling of the Southern Ocean, which is amplified by a decrease in CO2 and increase in sea ice. The lagged seesaw response seen in our results from Southern Hemisphere sites (Fig. 4b) is likely induced by oceanic mechanisms that delay and attenuate the transmission of the signal into the Southern Ocean (Keeling and Visbeck 2005; Schmittner et al. 2003; Stocker and Johnsen 2003). The heat content anomaly associated with the cooler SSTs in the Southern Ocean is rapidly transmitted equatorward by the atmosphere and the shallow meridional circulation in the Pacific basin, where it cools equatorial SSTs (Liu et al. 2002). The cooler equatorial Pacific SSTs decrease temperatures over Northern Hemisphere ice sheets, inducing ice growth and calving of ice and associated freshwater flux to the North Atlantic (Clark et al. 2007), which then leads to a decrease in the AMOC. This induces the ocean and atmospheric teleconnections to the reverse sign, with warming in the Southern Ocean and equatorial Pacific and an increase in atmospheric CO2 increasing temperatures over the Northern Hemisphere ice sheets. The attendant ice sheet melting, however, sustains a freshwater flux to the North Atlantic, eventually causing a collapse of the AMOC that culminates in a Heinrich event and maximum expression of the atmospheric and oceanic responses.
We note that meridional shifts in the position of the ITCZ in response to changes in the AMOC and sea ice extent (Chiang and Bitz 2005; Zhang and Delworth 2005) may short circuit these teleconnections, with a significant decrease in phase between changes in the AMOC and equatorial Pacific SSTs. Although records from the equatorial Pacific are currently unable to distinguish between these two scenarios (Clark et al. 2007), our finding that the southern signal is in phase with sea level (Fig. 5) supports the oceanic transmission mechanism in controlling the mass balance of Northern Hemisphere ice sheets.
The dominant high-frequency peaks found in both ice-core records (D-O cycles in Greenland, smaller A events in Antarctica) are found in the high-frequency SSA4 mode (Figs. 3c,d), suggesting that the process responsible for this mode had a global effect. We note that although the spectral power for this mode (∼1.6 kyr) is similar to the commonly described “1470 yr” period for D-O cycles that was derived from the GISP2 time scale (Stuiver et al. 1997), it does not have a sharp spectral peak, which is consistent with a mode of variability associated with changes internal to the climate system rather than a periodic signal requiring external astronomical forcing (e.g., Rahmstorf 2003). However, the high-frequency SSA mode of the northern and southern signals accounts for only about 5% of the signal’s variance, which is considerably less than the amount of variance at these frequencies found in just the GISP2 and Byrd ice cores (Roe and Steig 2004). Moreover, while this mode is equally important in both the southern and northern signals defined by our EOF analysis, only a few of the Southern Hemisphere records are coherent with respect to SSA4N. Specifically, the SSA4N mode shows strong coherence with records from Greenland, the tropical and North Atlantic Ocean, the Indian Ocean, and the Santa Barbara basin, but only shows significant coherence at two Antarctica sites (Fig. 4c). It is important to emphasize, however, that at least part of the absence of this mode in Southern Hemisphere records is likely associated with chronology uncertainties that become increasingly important at this higher frequency.
The spatial pattern of coherence among the individual sites (Fig. 4c) indicates that the origin of this mode is in the Northern Hemisphere but, unlike the low-frequency SSA mode, it is either not uniformly transmitted to the Southern Hemisphere or it is modulated by regional processes. One insight into how this mode is transmitted is provided by its short response time (Table 1). In particular, these results suggest that over much of the frequency band where variance is concentrated in this mode, its phase is not significantly different from zero, which appears more consistent with interhemispheric propagation of the signal by an atmospheric rather than oceanic mechanism. In this case, modulation of the signal might be expected by regional oceanic and atmospheric processes (Wunsch 2003), particularly where transmission of the signal itself is weak, resulting in its nonuniform geographic distribution. One well-known example of atmospheric transmission of a climate mode with a complex geographic pattern of response is provided by the ENSO cycle. In any event, given that there is little proxy evidence for coherent changes in the AMOC during D-O cycles (Clement and Peterson 2008), we speculate that changes in sea ice (Gildor and Tziperman 2003) and attendant atmospheric responses (Chiang and Bitz 2005) provide a possible mechanism to transmit this signal globally.
Our results thus support global transmission of millennial-scale signals, with the highest likelihood of the signal being recorded by a climate proxy occurring in association with strongly coupled changes in the climate system that overwhelm modification by regional processes. Nevertheless, we find that a substantial fraction of the variability in any one proxy record is not accounted for, and our EOF analyses indicate that even where global signals are present, they frequently represent a combination of the leading EOFs (Clark et al. 2007). Moreover, although we have found that two dominant modes characterize global-scale millennial-scale variability, if an individual proxy is a mixture of these climate modes, then we would not expect the proxy to be exactly in-phase with either of the dominant climate signals. In the absence of a robust and independent chronology, these results thus suggest caution in correlating proxy records to one of the well-known MIS3 signals by so-called wiggle matching. Experiments matching individual proxy records to these climate modes suggest that mean errors in chronologies as large as 400 yr can occur.
We thank Carl Wunsch for helpful discussions, and Mark Siddall and two reviewers for reviews. The NSF Paleoclimate Program supported this research.
Corresponding author address: Nicklas G. Pisias, College of Oceanic and Atmospheric Sciences, Oregon State University, Corvallis, OR 97331. Email: firstname.lastname@example.org