Ocean circulation, ice shelf, and sea ice interactions explain Dansgaard–Oeschger cycles
- aGeosciences Department and Laboratoire de Météorologie Dynamique (CNRS and Institute Pierre Simon Laplace, IPSL), École Normale Supérieure and Paris Sciences & Lettres (PSL) Research University, 75005 Paris, France;
- bResearch Domain IV–Transdisciplinary Concepts and Methods, Potsdam Institute for Climate Impact Research, 14473 Potsdam, Germany;
- cGrantham Institute, Imperial College, London SW7 2AZ, United Kingdom;
- dDepartment of Atmospheric and Oceanic Sciences, University of California, Los Angeles, CA 90095;
- eLamont-Doherty Earth Observatory, Columbia University, Palisades, NY 10964
See allHide authors and affiliations
Edited by Jean Jouzel, Laboratoire des Sciences du Climat et de L’Environ, Orme des Merisiers, France, and approved October 4, 2018 (received for review February 11, 2018)

Significance
Paleoclimatic proxy records from Greenland ice cores show that the last glacial interval was punctuated by abrupt climatic transitions called Dansgaard–Oeschger (DO) events. These events are characterized by temperature increases over Greenland of up to
Abstract
The last glacial interval experienced abrupt climatic changes called Dansgaard–Oeschger (DO) events. These events manifest themselves as rapid increases followed by slow decreases of oxygen isotope ratios in Greenland ice core records. Despite promising advances, a comprehensive theory of the DO cycles, with their repeated ups and downs of isotope ratios, is still lacking. Here, based on earlier hypotheses, we introduce a dynamical model that explains the DO variability by rapid retreat and slow regrowth of thick ice shelves and thin sea ice in conjunction with changing subsurface water temperatures due to insulation by the ice cover. Our model successfully reproduces observed features of the records, such as the sawtooth shape of the DO cycles, waiting times between DO events across the last glacial, and the shifted antiphase relationship between Greenland and Antarctic ice cores. Our results show that these features can be obtained via internal feedbacks alone. Warming subsurface waters could have also contributed to the triggering of Heinrich events. Our model thus offers a unified framework for explaining major features of multimillennial climate variability during glacial intervals.
- millennial climate variability
- nonlinear climate modeling
- north–south seesaw
- Bayesian parameter estimation
Since their discovery in oxygen isotope (
Variability of the last glacial interval as expressed by oxygen isotope ratios (
The concentration of
In addition to the North Greenland Ice Core Project (NGRIP) ice core, the general pattern of DO cycles has been observed in many Northern Hemisphere proxy archives, including other Greenland ice cores, but also, in terrestrial and marine sediments (5⇓–7). Most likely waiting times between subsequent DO events have been estimated to be of roughly 1,500 y (8), although there is no statistically significant evidence for a well-defined periodicity (9).
Corresponding
Since the DO events are outstanding examples of abrupt and dramatic climate transitions in the past, a better understanding of the underlying mechanisms is urgently needed to better assess the risk of abrupt climatic transitions in the future. Despite various promising advances, however, a comprehensive theory of the millennial-scale DO variability is still lacking. A physical theory that explains the DO cycles should reproduce (i) the sawtooth-shaped oscillations observed in
Most theories that have been proposed to explain the DO cycles focus on regime shifts in the Atlantic Meridional Overturning Circulation (AMOC) (17⇓–19), which transports water masses northward near the ocean’s surface and returns them southward in the deep ocean. A relationship between DO onsets and an AMOC strengthening has recently been inferred from a deep North Atlantic sediment core (20). There is, however, no general consensus whether existing paleooceanographic data suffice to support theories that propose AMOC shifts induced by freshwater pulses as triggers for DO events (21⇓–23). Furthermore, it is at least questionable how large-scale reorganizations of the oceanic circulation, which are likely to occur on much slower timescales (24), could cause the observed decadal-scale DO events on their own.
Variations in marginal ice sheets, ice shelves, and sea ice cover near Greenland and other North Atlantic basin coasts, possibly in concert with AMOC changes, have also been proposed to explain the observed DO cycles (23, 25⇓⇓⇓–29). Model results suggest that Nordic Sea sea ice retreat can increase winter temperatures by
Recently, a salt oscillator mechanism has been suggested as a possible explanation of the DO cycles obtained in a comprehensive climate model (43). This relaxation oscillator is triggered by Heinrich-type salinity disturbances, and the model results point to the formation of massive polynyas due to thermohaline convective instability, which in turn, leads to a rapid retreat of sea ice (44).
Alternatively, it has been suggested that warming subsurface waters in the northern North Atlantic can explain the DO events either by direct ice shelf melting from below (23, 41) or indirectly by destabilizing a proposed halocline (26, 45, 46). Empirical evidence for warming of subsurface waters before DO events is provided by planktonic and benthic foraminifera obtained from marine sediment cores (45). Benthic
A specific feature of the DO cycles, in particular for the longer ones, is that the cooling from interstadials to stadials takes place in two phases, with a rather slow decrease in the first phase followed by a considerably faster drop to stadial conditions in a second phase. A possible physical explanation for this behavior could be that, during the first phase, an ice shelf attached to Greenland grows at a relatively slow pace. After this ice shelf has reached a sufficient size, it can serve as an anchor for sea ice and thereby, provide favorable conditions for substantially faster sea ice expansion during the second phase (23, 41). For the shorter DO cycles, the ice shelf may not have been entirely removed during the previous transition from stadial to interstadial conditions; therefore, only the second phase of sea ice regrowth would be observed.
To our knowledge, studies focussing on sea ice or ice shelf variability to explain the DO events do not account for the antiphase coupling between Greenland and Antarctic temperatures (compare with Fig. 1 and ref. 12), and they do not account for the fact that, even in high northern latitudes, subsurface water temperatures are in phase with the temperature evolution observed in Antarctica (49). A possible explanation for these couplings is that reductions in North Atlantic subsurface water temperatures at the DO onset lead to a switch of the AMOC from its weak mode to its strong mode. Recent observational evidence of an approximately 200-y lag between Greenland and Antarctic temperatures, suggesting an oceanic north-to-south propagation of the climatic signal (12), supports this hypothesis. Note that, in such a setting, changes in AMOC strength would not be the cause but rather, a consequence of the DO events.
In this paper, we test the hypothesis that Greenland ice shelves and sea ice interact with the AMOC to produce the observed DO cycles and the shifted antiphase relationship between the two hemispheres. To do so, a conceptual model that comprises the key dynamical mechanisms involved is formulated in the next section.
Methods
Our model’s variables are the extent of sea ice cover C with
Typical stadial and interstadial conditions are sketched in Fig. 2, Upper and Lower, respectively, in accordance with available proxy evidence: During stadials, the extensive ice cover in the North Atlantic insulates the ocean from the atmosphere,
Schematic diagram of typical stadial and interstadial conditions during the last glacial interval together with the relevant observables included in the proposed model:
The coupled temporal evolution of
Such a conceptual model cannot be expected to obtain realistic temperatures and isotope ratios intrinsically, since it does not resolve energy fluxes locally. For the Greenland atmospheric temperatures
Fig. 2, Inset shows the feedback loop of the model proposed herein: subsurface water temperatures
The AMOC strength ψ is coupled to
The parameter
The AMOC strength ψ determines atmospheric temperatures in Antarctica and hence,
For the growth of the ice cover C, we impose a piecewise linear form: as long as the North Atlantic subsurface water temperatures
From the NGRIP
In accordance with the hypothesis to be tested here, the ice shelf growth speed
We infer the speed
When running our model thereafter for the entire glacial interval (see Fig. 4), we sample instead values for
The speed
After a time span with complete ice cover of the North Atlantic (i.e.,
The value of
The concentration of
In fact, such a cooling of source temperatures at the onset of DO events has been observed in the deuterium excess
For
Note that the specific choices of the fixed points for
The proposed model has the six above-mentioned variables
Optimal values for the four free parameters are determined via Approximate Bayesian Computations, a Bayesian approach to minimize the difference between observed and simulated summary statistics. For these summary statistics, the PDFs of the observed and simulated Greenland
Results
Our model is trained to reproduce the PDF of the NGRIP
Summary statistics of our DO cycle model. (A) Observed (blue) and simulated (red) PDFs of the Greenland ice core
Note that very similar values for
The dynamics simulated by our model is as follows. Starting in a situation with no ice cover, C first increases slowly at speed
Model simulation for the last glacial interval. (A) Benthic
Increasing
Regardless of the specific mechanism, the abrupt ice retreat causes the release of heat accumulated below the North Atlantic Ocean surface to the atmosphere, thereby rapidly increasing
Comparison between observed and simulated
Note that, for the simulation of the last glacial interval shown in Fig. 4, we have varied the positions of the fixed points
In agreement with observations, our model produces DO events only for glacial conditions and not for interglacial conditions (compare with Fig. 4) Furthermore, the frequency of DO events varied considerably across the last glacial; the blue line in Fig. 4B represents the observations. Our model quite accurately reproduces the temporal evolution of this frequency as shown by the red line and shading in Fig. 4B.
This result is due to two different mechanisms. First, the speed
We speculate that this prolonged AMOC slowdown around 20 ky b2k corresponds to the AMOC reduction inferred from more complex model simulations (61) and from
A direct comparison of the observed and simulated Greenland
Our model does not reproduce the high-frequency, large-amplitude oscillations that occurred shortly before some of the major DO events or the events that occurred toward the end of some interstadials. As noted above, these precursor and rebound events were excluded from the analysis beforehand, because they have been attributed to mechanisms other than those producing the main DO events (54).
The oceanic coupling mechanism between Greenland and Antarctica via the AMOC proposed herein provides a natural way to reproduce the observed antiphase relationship between Greenland and Antarctic temperatures, including a lagged response in Antarctica. Note that this mechanism differs from the minimal approach of coupling Greenland
Our results indicate that the direct coupling via the AMOC proposed herein is consistent with recent analyses based on high-resolution ice core records. The moderate warming in Antarctica during Greenland stadials and the moderate cooling during Greenland interstadials are reproduced in our simulations (compare with Fig. 4). The smoothing of the DO cycles in Antarctica is obtained in our simulations by the transmission of the signal from Greenland to Antarctica via the AMOC rather than by the heat reservoir of ref. 10. Furthermore, based on 1,000 simulated time series, we infer that the Antarctic maxima lag the Greenland DO events by
As noted above, our model is optimized only with respect to the PDF of the Greenland
Discussion
In addition to the statistical characteristics with respect to which we optimized the model, it reproduces other observed characteristics of DO cycles, including their sawtooth shape (Fig. 4D), the frequency of DO events across the last glacial interval (Fig. 4B), and the shifted antiphase relationship between
A key ingredient of our model is that heat transported northward by the AMOC accumulates below ice shelves and sea ice and eventually, removes the ice cover. So far, we only suggested this mechanism for a hypothesized Greenland ice shelf to trigger DO events; it could, however, apply—at least during some Greenland stadials—also for the Laurentide and Fennoscandian ice sheets or the ice shelves attached to either one of them. The ice sheets themselves could have been destabilized directly (60) or by the removal of the attached ice shelves. The resulting massive iceberg discharges would then cause the pronounced bands of ice-rafted debris in marine sediments that mark the Heinrich events (32⇓⇓⇓–36, 38, 41).
The proposed model thus qualitatively predicts the occurrence of Heinrich events as well and in particular, provides an explanation why Heinrich events occurred exclusively during Greenland stadial intervals. At least for the Laurentide ice sheet, a crucial factor in this context is the extensive sea ice cover of the Labrador Sea during Greenland stadials, which is documented by proxy data (65, 66). The presence of such an ice cover is necessary for subsurface warming to be relevant in destabilizing the Laurentide ice sheet, because the heat would otherwise be released to the atmosphere before reaching the Hudson Strait and grounding line of the Laurentide ice sheet.
The fact that Heinrich events occurred only during some Greenland stadials could be explained—along the lines of the hypothesis tested here—by the respective ice sheets and shelves having different geographical locations, with different strengths of the subsurface warming. Furthermore, the Laurentide and Fennoscandian ice sheets had, in all likelihood, longer characteristic timescales than the ice shelf proposed for eastern Greenland. In accordance with the binge–purge mechanism (67), a size threshold of the Laurentide and Fennoscandian ice sheets could still play a major role in the destabilization of these ice sheets and in their subsequent melting due to the subsurface warming.
It should be emphasized that, in such a setting, Heinrich events would not trigger DO events via the freshwater flux caused by iceberg calving. Instead, stadial conditions would lead to subsurface heat accumulation everywhere in the ice-covered North Atlantic. DO events are then caused as described above and are preceded by iceberg calving due to the destabilization of the Greenland ice shelf.
Moreover, it is important to note that, in our model setup, the changes in AMOC strength occur in response to the subsurface water temperature variations in the North Atlantic, which are induced by the ice shelf–sea ice mechanism that has been hypothesized heuristically in ref. 23 to explain the DO events. As already noted in the Introduction above, this stands in contrast to theories proposing AMOC variations as independent triggers for the DO events. Nevertheless, the AMOC does play a crucial role in helping our model to simulate DO events, since the northward oceanic transport of heat causes the subsurface heat accumulation underneath the ice cover during stadials, which eventually leads to the destabilization of the hypothesized ice shelf. Additional research, in particular on the relative dating between Greenland ice core records and oceanic tracers of the AMOC strength, will be needed to reach a consensus in this regard.
The main empirical support for the presence of an ice shelf to the east of Greenland comes from ice-rafted debris associated with each DO event: such debris is present in marine sediment cores from different parts of the adjacent seas (23, 39, 68, 69). These records suggest that the proposed ice shelf would have most likely been located to the east of Greenland, possibly narrowing the Denmark Strait during stadials.
Direct empirical evidence for the presence of this ice shelf during stadials is, however, so far not conclusive (42, 70⇓⇓–73). In principle, the two-step cooling after DO events could be explained by sea ice expansion alone if this expansion had occurred at two different speeds. The ice–albedo feedback (60) could have played a major role in this context: after initial slow sea ice expansion, the resulting cooling could have provided favorable conditions for subsequent significantly faster expansion. This alternative mechanism would also be consistent with our conceptual model. However, the abruptness of the switch from the slow to the fast cooling phase together with the proxy evidence of calving events for each Greenland stadial suggest that the ice shelf hypothesis is more likely.
The model proposed herein reproduces the key features associated with DO cycles during the last glacial interval, including the sawtooth shape of the DO oscillations in Greenland—with its two-step cooling from interstadials to stadials—as well as the correct waiting times between the DO events across the last glacial and the shifted antiphase relationship with Antarctica. No external forcing was included in our model to trigger the DO events, implying that these oscillations can be produced by internal feedbacks alone.
Since we have deliberately neglected salinity feedbacks in the formulation of our model, our results also show that these key features can, in principle, be explained by purely dynamic and thermodynamic arguments. The striking similarity between observed and simulated time series across the last glacial suggests that the ice shelf–sea ice mechanism, which has been proposed in earlier studies and is tested numerically with the model proposed herein, is a promising candidate for explaining the DO events as well as their connection to the Heinrich events. Our results emphasize, therefore, the need for additional investigations based on comprehensive ocean–atmosphere models with dynamical ice shelf and sea ice coupling.
Materials and Methods
Data.
We use a proxy record of oxygen isotope ratios (
In addition, in Fig. 1, we show, for comparison,
The NGRIP data used in this study are publicly available at www.iceandclimate.nbi.ku.dk/data/. The WAIS data shown in Fig. 1 were obtained from the supplementary information of ref. 12 (https://www.nature.com/articles/nature14401). The benthic
Parameter Estimation.
The quantity of interest in finding optimal parameter combinations π for a model given observed data D is the probability distribution
To avoid possible biases induced by Gaussian approximations of the likelihood function, we use instead so-called Approximate Bayesian Computations (78, 79) to find the parameters that optimize our model with respect to suitable summary statistics, namely the PDF of the NGRIP
i) Prescribe uniform prior distributions for the parameters with physically realistic ranges.
ii) Simulate time series with parameter combinations sampled from the uniform priors.
iii) Compute the PDF and the average stadial durations (SDs) for these simulated time series as summary statistics.
iv) Accept a parameter combination if the summary statistics of the simulation are within a prescribed tolerance of the observed statistics.
It can be shown that the joint PDF of the accepted parameter combinations approximates
Specifically, we accept a parameter combination
The joint PDF of the accepted parameter combinations as well as the optimal values of this joint PDF are shown in SI Appendix, Fig. S2. For our six-variable model, we have four parameters that need to take values from relatively narrow ranges to reproduce the statistics of the observed
For the summary statistics shown in Fig. 3 as well as for the simulations shown in Figs. 4 and 5, we used the parameter combination that maximizes the joint parameter PDF (blue lines in SI Appendix, Fig. S2).
Acknowledgments
N.B. acknowledges funding from the Alexander von Humboldt Foundation, the German Federal Ministry for Education and Research, and German Science Foundation Grant BO 4455/1-1. M.G. was partially supported by Multidisciplinary University Research Initiative of the Office of Naval Research Grant N00014-16-1-2073. This is Lamont–Doherty Earth Observatory (LDEO) Contribution #8260.
Footnotes
- ↵1To whom correspondence should be addressed. Email: boers{at}pik-potsdam.de.
Author contributions: N.B. designed research; N.B. performed research; N.B., M.G., and D.-D.R. analyzed data; and N.B. wrote the paper with contributions from M.G. and D.-D.R.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1802573115/-/DCSupplemental.
Published under the PNAS license.
References
- ↵
- Dansgaard W, et al.
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Rousseau DD, et al.
- ↵
- Schulz M
- ↵
- Ditlevsen PD,
- Andersen KK,
- Svensson A
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Henry LG, et al.
- ↵
- Piotrowski AM,
- Goldstein SL,
- Hemming SR,
- Fairbanks RG,
- Zylberberg DR
- ↵
- Pisias NG,
- Clark PU,
- Brook EJ
- ↵
- ↵
- Dijkstra HA,
- Ghil M
- ↵
- ↵
- ↵
- Ashkenazy Y,
- Losch M,
- Gildor H,
- Mirzayof D,
- Tziperman E
- ↵
- Menviel L,
- Timmermann A,
- Friedrich T,
- England MH
- ↵
- ↵
- ↵
- Timmermann A,
- Gildor H,
- Schulz M,
- Tziperman E
- ↵
- ↵
- ↵
- Hemming SR
- ↵
- ↵
- ↵
- Marcott SA, et al.
- ↵
- Bassis JN,
- Petersen SV,
- Mac Cathles L
- ↵
- Bond GC,
- Lotti R
- ↵
- ↵
- ↵
- Andrews JT,
- Dunhill G,
- Vogt C,
- Voelker AH
- ↵
- Peltier WR,
- Vettoretti G
- ↵
- Vettoretti G,
- Peltier WR
- ↵
- ↵
- Singh HA,
- Battisti DS,
- Bitz CM
- ↵
- ↵
- Hoff U,
- Rasmussen TL,
- Stein R,
- Ezat MM,
- Fahl K
- ↵
- ↵
- ↵
- Chen F,
- Ghil M
- ↵
- Djikstra H
- ↵
- ↵
- Capron E, et al.
- ↵
- Mitsui T,
- Crucifix M
- ↵
- ↵
- Steen-Larsen HC, et al.
- ↵
- Masson-Delmotte V, et al.
- ↵
- Parrenin F, et al.
- ↵
- Ghil M
- ↵
- Vettoretti G,
- Peltier WR
- ↵
- ↵
- Rypdal M
- ↵
- Boers N
- ↵
- Jennings AE, et al.
- ↵
- Simon Q,
- Hillaire-Marcel C,
- St-Onge G,
- Andrews JT
- ↵
- ↵
- Voelker AL, et al.
- ↵
- ↵
- ↵
- Funder S,
- Kjeldsen KK,
- Kjær KH,
- O Cofaigh C
- ↵
- Larsen NK, et al.
- ↵
- Bjørk AA, et al.
- ↵
- ↵
- ↵
- Busetto AG,
- Buhmann JM
- ↵
- Boers N, et al.
- ↵
- ↵
- Beaumont MA,
- Zhang W,
- Balding DJ
- ↵
- ↵
Citation Manager Formats
Article Classifications
- Physical Sciences
- Earth, Atmospheric, and Planetary Sciences