Data availability and temporal resolution make it challenging to unravel the anatomy (duration and temporal phasing) of the Last Glacial abrupt climate changes. Here, we address these limitations by investigating the anatomy of abrupt changes using sub-decadal-scale records from Greenland ice cores. We highlight the absence of a systematic pattern in the anatomy of abrupt changes as recorded in different ice parameters. This diversity in the sequence of changes seen in ice-core data is also observed in climate parameters derived from numerical simulations which exhibit self-sustained abrupt variability arising from internal atmosphere-ice-ocean interactions. Our analysis of two ice cores shows that the diversity of abrupt warming transitions represents variability inherent to the climate system and not archive-specific noise. Our results hint that during these abrupt events, it may not be possible to infer statistically-robust leads and lags between the different components of the climate system because of their tight coupling.
Paleoclimatic records of the Last Glacial reveal a series of abrupt warming events occurring in the North Atlantic region, known as Dansgaard-Oeschger (D-O) events, with counterparts in lower latitudes1 and Antarctic climate archives2,3. Oxygen isotope (δ18O) profiles from Greenland ice cores provide master records of this climate variability4,5, illustrating fluctuations between Greenland Stadial (GS) phases with full glacial conditions and milder Greenland Interstadial (GI) phases (Fig. 1). The D-O climate variability is commonly linked to changes in the intensity of the Atlantic meridional overturning circulation (AMOC), resulting in heat transport changes from the low to the northern high latitudes6,7. However, no consensus exists yet to explain what triggers the abrupt warmings, characterized by Greenland surface temperature increases of 5–16 °C within a few decades to centuries8. Among the proposed paradigms, mechanisms involving changes in Nordic Seas sea-ice cover9, atmospheric circulation10, or the collapse of ice shelves11 have been investigated. Recent studies suggest that abrupt climate variability can result entirely from unforced12 or noise-induced oscillations of the coupled atmosphere-ice-ocean system that alter poleward energy transport (ref. 13 and 14 for reviews).
The mechanisms proposed to explain D-O event dynamics can be confronted with annual-to-decadal-scale observations of climatic changes across the globe over the GS–GI transitions. Indeed, such data sets provide a basis to map out the sequence of events, infer possible causal relations and evaluate hypothetical sets of governing mechanisms by comparing model output with the spatial expression and relative phasing of the observed changes, hereafter referred to as the “anatomy” of the changes. However, looking at the anatomy of abrupt events in paleoclimate data is challenging because it requires a high temporal resolution not attainable in most climatic archives, and because of relative dating uncertainties between paleoclimate records from different archives. Records of annual or close-to-annual resolution from Greenland ice cores overcome this challenge since they contain tracers recording conditions in different parts of the Earth System with each year’s precipitation, all in one archive. The δ18O value of Greenland ice is mainly affected by local surface temperature changes, past changes in precipitation seasonality, the temperature at the moisture source regions, and elevation changes15,16,17,18. Hence, although δ18O is not a direct temperature proxy, it can be used as a qualitative tracer of local Greenland surface temperature changes. The second-order parameter d-excess (d-excess = δD−8·δ18O) is commonly interpreted as a record of past changes in evaporation conditions or shifts in mid-latitude moisture sources17,19,20, whereas Ca2+ concentrations ([Ca2+]) in Greenland ice cores reflect both source strength and transport conditions from terrestrial sources, which are mainly the mid-latitude Asian deserts21,22. Finally, changes in Na+ concentrations ([Na+]) can be interpreted as qualitative indicators of the sea-ice cover extent in the North Atlantic at the stadial-interstadial scale23, whereas relative site accumulation rate changes can be estimated from the annual-layer thickness (denoted λ)24. Hence, ice-core multi-tracer studies are well suited to evaluate the precise phasing and duration of changes between different regions without relative dating uncertainty as all records come from the same core.
This approach was initially applied to characterize the sequence of events at the onsets of the Holocene, GI-1e (Bølling), and GI-8c25,26. For each of those transitions, a lead of a few years in changes in terrestrial aerosol concentrations, accumulation rate, and mid-latitude moisture sources relative to the changes in marine aerosols and the isotopic temperature was found. Such results suggest that the Greenland surface warming was preceded by changes in the conditions at the dust sources or changes to the transport to Greenland (e.g., rainfall-driven changes in aerosol washout). In parallel, the phasing between the high- and lower-latitude climate responses was investigated using ice-core gas-phase measurements: the δ15N of N2 as a tracer for Greenland surface temperature changes27,28 and the methane concentration (CH4) as a proxy for tropical climate change29,30. Although the first studies29,31 estimated a lag of a few decades of tropical CH4 emissions behind δ15N at the onset of the abrupt warmings, a more recent study32, focusing on the Bølling transition and using 5-yr-resolution δ15N and CH4 records, estimated that high- and low-latitude climate changes occurred essentially synchronously at that time, with Greenland surface temperature leading atmospheric CH4 emissions by 4.5+21−244.5−24+21 yrs, in agreement within errors with ref. 25.
Benefiting from the new NGRIP and NEEM high-resolution ice-core data sets, recent work extended the multi-tracer approach developed by ref. 25 and 26 to all transitions back to 60 ka b2k (thousand years before 2000 C.E.) and derived an average sequence of changes characteristic of the GI onsets by combining the estimated leads and lags for all studied transitions23. Based on the assumption that the relative timing differences between different tracers at all GI onsets are the result of the same underlying process, it was found that changes in both local precipitation and terrestrial dust aerosol concentrations led the change in sea-salt aerosol concentrations and δ18O of the ice by about a decade. Event-stacking-based approaches are often applied to extract the common signal from highly variable climatic records33,34,35,36. Although this is useful, it is also worth looking into the details of the sequence of changes over each event, especially considering the high diversity observed in the amplitude of the warming8, the shape and duration of GS and GI37,38 (Fig. 1), and the evolving climatic background state throughout the Glacial (orbital configuration, global ice volume, and atmospheric greenhouse gas concentrations). Taking this view, we observe that the results from ref. 23 illustrate a decadal-scale range in leads and lags from one event to the next when considering the onset of each individual transition. These differences can be interpreted as coming from different realizations of the same set of underlying mechanisms owing to noise processes in the archive and internal variability in the climate system, or alternatively as a suggestion that one common set of mechanisms or sequence of events may not adequately describe the processes of all rapid warming transitions.
The aim of this study is twofold. First, we investigate the anatomy of the D-O warming transitions down to 112 ka b2k using a multi-tracer approach relying on new and existing records from the Greenland NEEM (77.45°N, 51.08°W) and NGRIP (75°N, 42.3°W) ice cores. Having so many highly resolved ice-core records from two different locations over numerous D-O events provides the most comprehensive opportunity so far to assess the geographical representativeness of single ice-core records. Second, the anatomy of D-O warmings inferred from Greenland ice-core data is compared with new simulations from the coupled Community Climate System Model Version 4 (CCSM4) as the basis for discussing the processes involved in D-O warmings.
We use here new and existing water isotope measurements (δ18O, d-excess) at high resolution (5 cm) from the NGRIP ice core5 (Supplementary Data 1). The temporal resolution of the measurements corresponds to 1, 3, 4, 5 yr per sample at 10, 45, 80, and 105 ka b2k, respectively. We also include in our analysis, sections from the recent NEEM high-resolution water isotope records39 for which the 5 cm resolution corresponds to 1, 4, 7, 18 yr per sample at 10, 45, 80, and 105 ka b2k. We also present high-resolution NGRIP and NEEM [Ca2+] and [Na+] records annually interpolated and extended back to ~108 ka b2k (Methods, Supplementary Data 2). Finally, we use the NGRIP λ record back to 60 ka b2k obtained from the GICC05 annual-layer counting based on aerosol and visual stratigraphy records (Supplementary Data 3). We restrict our λ analyses to the last 60 ka as λ is modeled from the stable water isotope record below this age and, therefore, is not independent of δ18O. The GICC05 chronology is applied to NEEM by means of interpolation between reference horizons of mainly volcanic origin40. The NEEM annual-layer thicknesses are only available as averages between these unevenly spaced reference horizons, rendering the NEEM λ record unsuitable for this study. The NGRIP and NEEM data sets are reported on the GICC05 chronology back to 60 ka b2k and on the flow model-extended GICC05modelext chronology below this40,41. Age interpolation uncertainties limit the direct comparison of the absolute timing of changes between cores40.
We use a probabilistic characterization of the transitions to infer the timing, duration, and amplitude of the local and regional changes associated with each studied D-O warming. Following refs. 23,25, we determine the relative phasing of changes in the different data sets by fitting a ramp (i.e., a linear change in the raw or logarithmically-transformed data between two stable states) to each data series within a prescribed search interval across each GS–GI transition (Supplementary Figure 1, Supplementary Table 1, Supplementary Data 4). We describe the ramp by the temporal midpoint of the ramp, the duration of the transition, the data value before the transition, and the amplitude of the change. Our probabilistic model also accounts for additive noise with autocorrelation (Methods). Note that our method is conceptually similar to ref. 23 with only minor differences in the parameter priors, whereas the uncertainty estimation is different from that employed by ref. 25, which used the RAMPFIT method42. In the following, we only display results for transitions where the ramp-fitting technique provides an unequivocal solution, i.e., the timing and duration of the identified onset and end of the transitions do not change by more than a decade when the width of the search time window is varied (Methods, Supplementary Figure 3).