Research Article

Ocean circulation, ice shelf, and sea ice interactions explain Dansgaard–Oeschger cycles

See allHide authors and affiliations

PNAS November 20, 2018 115 (47) E11005-E11014; first published November 1, 2018; https://doi.org/10.1073/pnas.1802573115

  1. 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 15° C within a few decades. The cause of these transitions and their out-of-phase relationship with corresponding records from Antarctica remains unclear. Based on earlier hypotheses, we propose a model focusing on interactions between ice shelves, sea ice, and ocean currents to explain DO events in Greenland and their Antarctic counterparts. Our model reproduces the main features of the observations. Our study provides a potential explanation of DO events and could help assess more accurately the risk of abrupt climatic transitions in the future.

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.

Since their discovery in oxygen isotope (δ18O) profiles derived from Greenland ice cores (1⇓–3), substantial climatic fluctuations during the last glacial period, roughly 110–10 kyB.P., have attracted great interest in the paleoclimate community. These fluctuations are characterized by rapid, decadal-scale transitions from stadial cold to interstadial warm conditions, which are commonly termed Dansgaard–Oeschger (DO) events, followed by slow relaxations back to stadials within centuries to millennia (Fig. 1).

Fig. 1.
Fig. 1.

Variability of the last glacial interval as expressed by oxygen isotope ratios (δ18O). Blue indicates Greenland δ18O data obtained from the NGRIP (3) at a regular sampling rate of 20 y (13). Orange indicates Antarctic δ18O data from the WAIS ice core (12). As in ref. 12, the layer-counted NGRIP chronology GICC05 (14, 15) is rescaled by a factor of 1.0063, because the layer-counted WAIS divide deep ice core chronology (WD2014) (16), on which the WAIS δ18O record is shown, is synchronized to this rescaled chronology. The δ18O values are commonly interpreted as a proxy for atmospheric temperatures at the location of the ice core, with higher values indicating warmer temperatures. The training period for our model is the interval from 59 to 23 ky b2k, which roughly corresponds to Marine Isotope Stage 3 (MIS3). DO events are indicated by vertical magenta lines, Heinrich stadials are marked by grey shading, and MISs are indicated at the top of the figure. The thin vertical dashed lines indicate time steps in intervals of 20 ky. Inset shows the geographical locations of the NGRIP and WAIS sites and a sketch of the oceanic circulation, with warmer surface flow in red and colder bottom flow in blue.

The concentration of δ18O in ice cores is commonly interpreted as a qualitative proxy for the atmospheric site temperature. By combining δ18O with δ15N measurements from the same ice cores, in conjunction with firnification and heat diffusion models, increases of up to 16.5° C for single DO events have been indirectly estimated (cf. ref. 4 and references therein).

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 δ18O time series obtained from Antarctic ice cores show an antiphase relationship with the temporal evolution in Greenland, with gradual increases during Greenland stadials and gradual cooling during Greenland interstadials (10⇓–12) (compare with Fig. 1). A recent study based on a high-resolution Antarctic ice core estimates a delay of roughly 200 y between DO events in Greenland and the onset of cooling in Antarctica as well as between the return to stadial conditions in Greenland and the onset of warming in Antarctica. From the sign and magnitude of this delay, an oceanic north-to-south transmission of the climatic signal has been inferred (12).

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 δ18O time series from Greenland ice cores, (ii) the changing frequency of DO events across the last glacial, and (iii) the observed antiphase relationship between observed ice core δ18O data in the two hemispheres.

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 10° C (30). Furthermore, it has been shown that freshwater pulses induced by iceberg discharges can trigger DO-type oscillations via coherence resonance (31). The latter study also provides a possible explanation for the suggested relationship between DO events and Heinrich events (32⇓–34), which are characterized by massive iceberg discharges into the Labrador Sea. These discharges are evident as pronounced bands of ice-rafted debris in marine sediment cores (34). The Heinrich events themselves might have been triggered by warming subsurface waters in the northern North Atlantic during stadial conditions (35⇓⇓–38). Heinrich-type iceberg calving occurs during the cooler stadials, possibly acting as a feedback stabilizing the stadial conditions (39, 40). It should be noted here that—in addition to the Heinrich events, during which icebergs were mainly discharged into the Labrador Sea—there is empirical evidence for substantial iceberg discharges at several other locations around the northern North Atlantic (41) and in particular, into the Denmark Straight and Icelandic Sea (23, 42).

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 δ18O can, via its relationship to brine production (47), serve as an indirect proxy for the variability of sea ice cover. Data obtained from North Atlantic marine sediment cores indicate that strong sea ice cover changes are associated with the DO cycles: extensive sea ice cover occurs during stadials, while interstadials have substantially reduced sea ice cover (26, 48).

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 0≤C≤1, Greenland oxygen isotope ratios IG, near-surface atmospheric temperature over Greenland TG, and subsurface water temperature TNAW in the North Atlantic close to Greenland as well as the strength ψ of the AMOC and the Antarctic oxygen isotope ratios IA.

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, TG is low, TNAW is rising, the AMOC ψ is weak, and IA is high. During interstadials, in contrast, the ice cover is (at least partially) removed, TG is high, TNAW is dropping, the AMOC is strong, and IA is reduced.

Fig. 2.
Fig. 2.

Schematic diagram of typical stadial and interstadial conditions during the last glacial interval together with the relevant observables included in the proposed model: δ18O in Greenland cores (IG), atmospheric temperature in Greenland (TG), subsurface water temperatures in the northern North Atlantic (TNAW), extent of ice cover close to Greenland (C), AMOC strength (ψ), and δ18O in Antarctic cores (IA); details are in the text. The blue arrows indicate the distinct atmospheric paths of δ18O from the evaporative source in the Atlantic Ocean to the ice core site in northern Greenland. Small and large fonts correspond to low and high values, respectively, of the observables. Inset shows the feedback mechanisms involved; here, solid arrows indicate a positive or enhancing influence, while dashed arrows indicate a negative or damping influence.

The coupled temporal evolution of TG, TNAW, ψ, and IA is given byṪG=−p0(TG−TG*)+10⋅p0TNAWΘ(TNAW)(1−C),[1]ṪNAW=−p0(TNAW−TNAW*)−10⋅p0TNAWΘ(TNAW)(1−C)−p0ψ,[2]ψ̇=ψ−ψ3−p0(TNAW−τ/2),[3]İA=−(IA−(IA*−ψ)),[4]where ψ̇=dψ/dt and so on.

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 TG, we, therefore, impose a fixed point TG*, around which the model dynamics takes place. The position of TG* in the model’s phase space is taken directly from paleoclimatic reconstructions, namely TG*=−51.5° C (4). Similarly, for Antarctic ice core δ18O, we impose a fixed point at IA*=−42‰ (12). To the best of our knowledge, however, no reconstructions exist for a fixed point TNAW* of North Atlantic subsurface water temperatures; its choice will be explained further below.

Fig. 2, Inset shows the feedback loop of the model proposed herein: subsurface water temperatures TNAW act positively on atmospheric temperatures TG via heat release from the upper ocean to the atmosphere. The sea–air heat exchange in the North Atlantic is proportional to TNAW, and it stops after TNAW falls below zero as indicated by the Heaviside function Θ, with Θ(x)=0 for x≤0 and Θ(x)=1 for x>0. The sea–air heat exchange is modulated by the ice cover C, because extensive ice cover insulates the atmosphere from the ocean (compare with Eq. 1).

The AMOC strength ψ is coupled to TNAW via a simplified version of the Stommel model (50), allowing for two alternative stable states of the AMOC (compare with Eqs. 2 and 3). We focus here on the negative thermal feedback between ψ and TNAW and neglect the positive salinity feedback of the Stommel model as in ref. 51; doing so allows one to investigate whether DO-type dynamics can be obtained without taking into account salinity changes (cf. ref. 52).

The parameter p0 sets the timescale at which the different variables are pulled toward their respective fixed points, and we set the timescale of the heat exchange between ocean and atmosphere to be an order of magnitude faster (compare with Eq. 1). The precise value of this prefactor does not impact the simulation results as long as it is much larger than p0.

The AMOC strength ψ determines atmospheric temperatures in Antarctica and hence, IA (compare with Eq. 4). For positive ψ, which corresponds to the strong AMOC mode, heat transport away from Antarctica is enhanced, and IA decreases accordingly. For negative ψ, which corresponds to the weak AMOC mode, IA increases, because less heat is transported away from the Southern Ocean.

For the growth of the ice cover C, we impose a piecewise linear form: as long as the North Atlantic subsurface water temperatures TNAW stay below a threshold τ (i.e., as long as TNAW<τ), we haveĊ=s1, ifC<C¯,s2, otherwise.[5]Here, s1 denotes the speed at which the proposed ice shelf grows, and s2 denotes the speed at which sea ice expands after the ice shelf has formed. Note that 0≤C≤1, with C=0 representing no ice cover and C=1 representing complete ice cover in the North Atlantic.

From the NGRIP δ18O, it can be directly inferred that the switch from the slow to the fast cooling phase of the DO cycles occurs, on average, after roughly one-half of the δ18O increase at the prior DO event has been lost. In the proposed model setup, based on the hypothesis formulated by ref. 23, these slow and fast cooling phases correspond to the ice shelf and sea ice growth phases, respectively. The growth of the ice cover C thus switches from the lower speed s1 to the higher speed s2 at C¯=0.5. Varying C¯ within the interval [0.4,0.6] does not affect the results of our simulations.

In accordance with the hypothesis to be tested here, the ice shelf growth speed s1 is inversely proportional to the interstadials’ duration. When comparing the durations of interstadials throughout the last glacial with the temporal evolution of benthic δ18O (53)—as a proxy for the global ice volume and thus, the background global mean temperature—a clear nonlinear relationship arises, with longer interstadial durations for periods with lower global ice volume and shorter ones for higher ice volume (SI Appendix, Fig. S1). The benthic δ18O time series is low-pass filtered at (10 ky)−1 to remove any potential influence of the millennial-scale DO variability. Note that we excluded from the list of transitions to interstadials the so-called precursor and rebound events, which have been attributed to mechanisms other than those producing the main DO events (54).

We infer the speed s1 from the observed relationship between the interstadial durations and benthic δ18O using a 2D Gaussian kernel density estimator, which is illustrated in SI Appendix, Fig. S1C. For each simulated interstadial, we sample a value for s1 from this estimated joint probability density function (PDF). For the training period of our model from 59 to 23 ky b2k, we sample values for s1 from the joint PDF restricted to all benthic δ18O values for this time interval; these values are between the blue vertical lines in SI Appendix, Fig. S1C.

When running our model thereafter for the entire glacial interval (see Fig. 4), we sample instead values for s1 from the joint PDF given the benthic δ18O value at the respective time step. Note that, in doing so, we do not use the benthic δ18O as external forcing to trigger the DO events but only choose the growth speed of the proposed ice shelf in accordance with the background global mean temperature as represented by benthic δ18O. It has been previously shown that the interstadial durations exhibit a more robust relationship with the benthic δ18O than with the solar insolation curve (55). For benthic δ18O values higher than the maximum of the entire last glacial, we set the ice shelf growth speed s1 to zero. Physically, this setting corresponds to global background climatic conditions becoming so warm that the ice shelf will not grow any longer. Specifically, we use a threshold of −3.6 permil; varying this threshold between (i) the maximum of benthic δ18O for the last glacial and (ii) its minimum for the Holocene does not alter our simulation results.

The speed s2 at which sea ice expands is, according to the proposed hypothesis, related to the fast cooling phase from interstadials to stadials. In our simulations, each integration step corresponds to a time interval of 20 y in accordance with the temporal resolution of the NGRIP δ18O record used herein. In this record, the fast cooling occurs within at most a century, and we accordingly set s2 such that the ice cover goes from C=0.5 to C=1 within five simulation time steps.

After a time span with complete ice cover of the North Atlantic (i.e., C=1), the subsurface warming will destabilize the ice shelf and thereby, remove the entire ice cover (23); we thus set C=0 whenever TNAW≥τ. For the latter threshold, we chose τ=10° C to account for the fact that Greenland temperature increases during DO events are, on average, of this magnitude (4).

The value of TNAW*, the fixed point of North Atlantic subsurface water temperatures, needs to be smaller than τ for the AMOC to contribute to an extra increase in TNAW, up to levels at which the proposed Greenland ice shelf is destabilized. For the model training phase, we set TNAW*=9° C. Varying the values of τ and TNAW* within a range of 3° C does not affect the main results of our model simulations. Furthermore, this fixed point should not be interpreted as the effective temperature of subsurface waters in the North Atlantic, which were most likely considerably colder. The effective values of TNAW resulting from our model simulations are substantially lower due to the release of heat to the atmosphere (see Fig. 4).

The concentration of δ18O in Greenland ice cores should not be regarded as a direct proxy of atmospheric temperatures at the site of the core. A simple model of Rayleigh distillation implies that δ18O depends on the temperature gradient between the evaporation source (i.e., the midlatitude Atlantic Ocean) and the site temperatures (56, 57). A substantial retreat of the North Atlantic ice cover C, which is supposed to occur for the DO events, reduces the source temperatures, because it shifts the source region northward (blue arrows in Fig. 2). This shift results in the effective cooling of source temperatures.

In fact, such a cooling of source temperatures at the onset of DO events has been observed in the deuterium excess d=δD−8δ18O (58). Therefore, both Greenland temperatures TG and the extent of the ice cover C impact the δ18O concentration in Greenland precipitation; accordingly, we model IG as a superposition of effects due to TG and C:IG(t+1)=IG*+p1TG(t)+p2(1−C(t))+p3ηt.[6]As for TG and IG above, we impose a base value IG*=−42.7‰, which is estimated directly from the NGRIP δ18O record. In the latter equation, ηt denotes Gaussian white noise with unit standard deviation σ.

For IA, the Antarctic δ18O, we assume a direct linear relationship to atmospheric temperatures, which can be inferred by comparing the West Antarctic Ice Sheet (WAIS) δ18O record with the Antarctic temperature stack (12, 59). Physically, the northward flow of the AMOC at shallow ocean levels implies that heat is always transported northward. In our model, melting by subsurface heat accumulation occurs, therefore, only in the North Atlantic.

Note that the specific choices of the fixed points for TG, IG, and IA do not impact the dynamical properties of our model: they are simply base values around which the variability is modeled to occur. When running our model for the entire glacial, we change the values of these fixed points in accordance with the slow variations in global ice volume, which serves as a proxy for the low-frequency global mean temperature evolution as given by the benthic δ18O (53).

The proposed model has the six above-mentioned variables (C,IG,TG,TNAW,ψ,IA) and four free parameters (p0,…,p3) that cannot be estimated directly from paleoclimatic reconstructions. Of these four parameters, only p0 is dynamically relevant (compare with Eqs. 1–5), while the remaining three (p1, p2, and p3) are merely used to estimate oxygen isotope ratios IG from Greenland atmospheric temperatures TG and ice cover C (compare with Eq. 6). The latter three parameters are thus only used to translate the model variables TG and C to the observed proxy variable δ18O that is actually measured in the NGRIP ice core (i.e., to establish the link between model variables and proxy observations).

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 δ18O time series as well as the observed and simulated average durations of stadials are taken into account (compare with SI Appendix). The training period for which the observed PDF and stadial durations are computed is the time interval from 59 to 23 ky b2k, which roughly corresponds to Marine Isotope Stage 3 (Fig. 1).

Results

Our model is trained to reproduce the PDF of the NGRIP δ18O time series for the training period from 59 to 23 ky b2k (Fig. 3A) as well the average stadial durations in this interval. The model also accurately reproduces the total waiting times between subsequent DO transitions (i.e., the stadial plus interstadial durations) (compare with Fig. 3B). The latter result can be explained from the fact that the correct average interstadial durations are obtained via sampling the growth speed s1 of the ice shelf, according to the observed relationship with the global ice volume, used as a proxy for global mean temperature variations at multimillennial timescales (SI Appendix, Fig. S1).

Fig. 3.
Fig. 3.

Summary statistics of our DO cycle model. (A) Observed (blue) and simulated (red) PDFs of the Greenland ice core δ18O time series. The blue and red shadings show the 2σ range of uncertainties for the observed and simulated PDFs, respectively. For the observed PDF, uncertainty estimates are derived by Monte Carlo sampling 200 times a subset of 1,000 data points of the total of 1,713 points. For the simulated PDF, uncertainties are derived from 1,000 simulated δ18O time series using the optimal combination of parameters (blue lines in SI Appendix, Fig. S2). (B) Observed (blue dashed–dotted line) vs. simulated (red dashed line) average stadial durations (SD); note that the two are practically indistinguishable. The red dotted lines indicate the 25th and 75th percentiles of the distribution of the simulated stadial durations shown in solid red. (C) Same as B but for the waiting times (WT) between subsequent DO transitions. Note that our model is trained with respect to the PDF shown in A and the average stadial durations shown in B. The waiting times shown in C are given by the sum of the stadial and interstadial durations as inferred from their relationship with the global ice volume (SI Appendix, Fig. S1), and they are not used for parameter optimization.

Note that very similar values for s1 would be obtained by including this parameter in the Bayesian optimization; still, to emphasize their physical meaning, we chose to sample the values for s1 directly from the joint PDF estimate with the global ice volume. The optimal parameter combination is (p0,p1,p2,p3)=(5.55,0.01,4.38,11.28). Here, p0 sets the timescale of the model’s dynamic evolution, and it is thus directly related to the simulated duration of stadials. The parameters p1, p2, and p3 are related to the PDF of IG and optimized to bring it close to the PDF of the NGRIP δ18O.

The dynamics simulated by our model is as follows. Starting in a situation with no ice cover, C first increases slowly at speed s1, consistent with a growing ice shelf attached to the Greenland ice sheet. After reaching a certain size, this ice shelf produces favorable conditions for extensive sea ice formation, which are further enhanced by the ice–albedo feedback (23, 60), implying a more rapid increase of C at speed s2 as shown by the two-step increase of C in Fig. 4C. An extensive ice cover, however, prevents the release of heat from the ocean to the atmosphere in the North Atlantic and leads to the relatively low TG values in Fig. 4E as well as to an accumulation of heat below the ice-covered sea surface (compare with Fig. 4F).

Fig. 4.
Fig. 4.

Model simulation for the last glacial interval. (A) Benthic δ18O (cyan) (53) used to determine the growth speed s1 of the proposed ice shelf (SI Appendix, Fig. S1). (B) Number of DO events in a 20-ky-wide moving window. The observed numbers are shown in blue, and the simulated numbers are in red, with shading indicating ±1 standard deviation σ as obtained from 1,000 independent simulations. (C) North Atlantic ice cover C (grey shading). (D) Greenland oxygen isotope ratios IG (light blue) together with a 200-y running mean (dark blue). (E) Greenland atmospheric temperatures TG (dark green), where white noise with standard deviation σ=1.5 permil has been added (light green variations around the dark green curve) to facilitate visual comparison with proxy reconstructions (4). (F) Subsurface North Atlantic water temperatures TNAW (red). (G) The strength ψ of the AMOC (black). (H) Antarctic oxygen isotope ratios IA (dark orange), where white noise with standard deviation σ=0.4 permil has been added (light orange) for comparison with the corresponding proxy values shown in Fig. 1. These illustrative time series are obtained by integrating our model with the optimal combination of parameters (blue lines in SI Appendix, Fig. S2) as determined using Approximate Bayesian Computations (compare with SI Appendix). The smooth thin curves in D–F and H indicate the slow multimillennial evolution of the fixed points of the respective variables, around which the dynamics takes place (compare with Eqs. 1, 2, 4, and 6). DO events are marked by vertical magenta lines. The time unit corresponds to 20-y intervals in accordance with the resolution of the published observational data (13).

Increasing TNAW, in turn, makes the subsurface waters lighter and prevents them from sinking to deeper ocean levels; thereby, the AMOC is weakened as seen in Fig. 4G. As soon as TNAW reaches a threshold τ, the ice cover is quickly removed (i.e., C is set to zero in the model (compare with Fig. 4C). Specifically, the rapid ice retreat could be either caused by direct melting of the proposed ice shelf from below (23) or caused by a destabilization of the water column and the removal of a halocline that was present during stadial conditions (26, 44).

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 TG, which is strikingly apparent in Fig. 4E. This rapid increase marks the DO transition from stadial to interstadial conditions in Greenland, and it is indicated by the vertical magenta lines in Fig. 4; these lines appeared already in Fig. 1 and will be used again in Fig. 5. The resulting cooling of North Atlantic subsurface water temperatures TNAW leads to a switch of the AMOC back to its strong mode and hence, to a relative cooling of Antarctic air temperatures and thus, IA (compare with Fig. 4H).

Fig. 5.
Fig. 5.

Comparison between observed and simulated δ18O time series for Greenland and Antarctica. (A) Observed time series of δ18O obtained from the NGRIP ice core with 20-y time steps from ref. 13 on the GICC05 chronology rescaled by a factor of 1.0063 to synchronize it with the WD2014 chronology of the WAIS δ18O record (cf. ref. 12). (B) Observed time series of δ18O obtained from the WAIS ice core on the WD2014 chronology (12). (C) Simulated time series IG of Greenland δ18O. (D) Simulated time series IA of Antarctic δ18O. To facilitate the visual comparison with the observed record shown in B, white noise with a standard deviation σ=0.4 has been added to the simulated time series of IA. The vertical magenta lines indicate the simulated DO events as in Figs. 1 and 4, while the vertical cyan lines mark the relative maxima of the simulated time series of IA. Based on 1,000 simulated time series, we infer that the Antarctic maxima lag the Greenland DO events by 315±33 y, which is just slightly larger than the value of 218±92 y reported in ref. 12. As in Fig. 4, each time step of the two simulated time series corresponds to a 20-y interval, and the same optimal parameter combinations are used. The time range corresponds to the temporal coverage of the WAIS record.

Note that, for the simulation of the last glacial interval shown in Fig. 4, we have varied the positions of the fixed points IG*, TG*, TNAW*, and IA* according to the multimillennial global temperature variability by settingf*(t)=a−bδ18Obenthicmax{δ18Obenthic}−min{δ18Obenthic}.Here, f* denotes the fixed point in question, and a and b are adjusted for each variable such that the observed differences between glacial and interglacial conditions are met (thin lines in Fig. 4 D–F and H).

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 s1 of the ice shelf growth, which is inversely proportional to the duration of interstadials, is sampled as a function of the global ice volume taken as a proxy of the background climate. Technically, this relationship is quantified by the 2D Gaussian kernel estimate of the joint density between s1 and benthic δ18O as described above. Second, during the last glacial maximum, around 20 ky b2k, the fixed point TNAW* of the subsurface North Atlantic water temperatures is so low that, even with the additional warming contributed by the prolonged AMOC slowdown, TNAW does not cross the threshold τ to destabilize the ice shelf.

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 Pa231/Th230 reconstructions (62). It should be noted, however, that, in the latter proxy evidence, the AMOC reduction starts roughly 2,000–3000 y later than in our simulations. Furthermore, we note that our model produces, on average, two to three transitions in the time interval between the last glacial maximum and the transition to the Holocene, although there is only one major event in the observed record (compare Fig. 5 A and C). This explains the slight offset between observed and simulated DO event frequencies between 35 and 10 ky b2k (Fig. 4B).

A direct comparison of the observed and simulated Greenland δ18O time series in Fig. 5 A and B, respectively, shows that they are in very good qualitative agreement. In particular, our model reproduces the two distinct steps in the cooling from interstadials to stadials. After a DO event, δ18O and atmospheric temperatures in Greenland decrease slowly during the first phase and more rapidly back to stadial conditions in the second phase (23, 41). The inclusion of a slowly growing ice shelf in the first phase in addition to rapidly expanding sea ice in the second phase is essential for our model’s reproduction of this feature of the observations. Setting C¯=0 in Eq. 5 would correspond to a model that only includes sea ice expansion at the fast speed s2, while setting C¯=1 would correspond to a model with only ice shelf growth at slow speed s1. Both settings would lead to time series with only one phase of cooling, which is not consistent with the observations.

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 δ18O variability directly to a conceptual heat reservoir representing the Southern Ocean (10), although the two approaches do not contradict each other. The latter approach provided an explanation for the smoothing of the signal in Antarctica via a conceptual heat reservoir, but it did not aim to reproduce the temporal lag, which was discovered only recently (12).

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 315±33 y in our model, which is reasonably close to the value of 218±92 y estimated from the corresponding observations (12).

As noted above, our model is optimized only with respect to the PDF of the Greenland δ18O record and the average duration of stadials, which were both evaluated during the training period of 59–23 ky b2k. The similarity between the observed (Fig. 4B) and simulated (Fig. 4D) Antarctic δ18O is, therefore, unexpected and provides additional support for the proposed coupling mechanisms.

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 δ18O in Greenland and Antarctica (Fig. 5 C and D). Additional evidence for the hypotheses tested here is provided by recent results showing significant variance and autocorrelation increases before DO events in the high-frequency band of the NGRIP δ18O time series (63, 64). A destabilized, fluctuating ice cover before a DO event would very likely lead to high-frequency changes in the difference between source and site temperatures. Since this difference is directly proportional to the δ18O values obtained from ice cores, a fluctuating ice cover could, therefore, mechanistically explain the statistical precursor signals.

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 (δ18O) from the NGRIP. The δ18O ratios obtained from ice cores are commonly interpreted as a proxy for atmospheric temperature variability at the site of the ice core (2, 3, 74, 75). The layer-counted chronology of this record is the Greenland Ice Core Chronology 2005 (GICC05) (15), which starts at roughly 60 ky b2k. We use the most recently published version of this record, with values of δ18O given at a temporal resolution of 20 y (13).

In addition, in Fig. 1, we show, for comparison, δ18O data from the WAIS ice core (12). Furthermore, we used the stack of benthic δ18O reported in ref. 53 to estimate the multimillennial variability of global ice volume across the last glacial interval. This time series has a temporal resolution of 1 ky; it is first low-pass filtered at (10 ky)−1 to exclude any impact of the millennial-scale DO variability and thereafter, linearly interpolated to the GICC05 chronology.

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 δ18O stack data can be downloaded from lorraine-lisiecki.com/LR04stack.txt.

Parameter Estimation.

The quantity of interest in finding optimal parameter combinations π for a model given observed data D is the probability distribution P(π∣D) of π given the observations D. Bayes’ Theorem implies that this distribution can, in principle, be obtained from P(D∣π), the likelihood function of the model, given the observations. The common Maximum Likelihood Estimation method is based on the fact that the parameters that maximize the likelihood function are also the parameters that maximize P(π∣D). The precise form of the likelihood function is, however, very hard and often impossible to guess, and existing approaches, therefore, commonly resort to Gaussian approximations (76, 77).

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 δ18O time series and the average stadial duration (SD). This optimization is implemented by the following rejection algorithm.

  • 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 P(π∣D) if the chosen summary statistics are sufficiently informative (80, 81). We thus take the parameter combination that maximizes the joint parameter PDF as the best estimate of the optimal parameter combination.

Specifically, we accept a parameter combination (p0,…,p3) if the following two conditions are fulfilled:supI|Pobs(IN)−Psim(IN)|<0.05,and|SD¯obs−SD¯sim|SD¯obs<0.01;here, Pobs denotes a Gaussian kernel density estimate of the PDF of the observed NGRIP δ18O time series, and Psim is the corresponding estimate of the PDF of the simulated IN time series using the parameters (p0,…,p3). Furthermore, SD¯ denotes the average stadial duration. The training period for which the observed summary statistics are computed is the interval from 59 to 23 ky b2k.

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 δ18O time series.

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.

Back to top
Article Alerts
Email Article
Citation Tools
Request Permissions
Ocean circulation, ice shelf, and sea ice interactions explain Dansgaard–Oeschger cycles
Niklas Boers, Michael Ghil, Denis-Didier Rousseau
Proceedings of the National Academy of Sciences Nov 2018, 115 (47) E11005-E11014; DOI: 10.1073/pnas.1802573115
del.icio.us logo Digg logo Reddit logo Twitter logo CiteULike logo Facebook logo Google logo Mendeley logo

Article Classifications

Proceedings of the National Academy of Sciences: 115 (47)
Table of Contents

Submit

Jump to section

You May Also be Interested in

Sunset.
Resolving the red sky paradox
A study suggests that red dwarfs are not visible in the sky due to factors that inhibit complex life to develop on or around them.
Image credit: Pixabay/pixel2013.

Similar Articles