A mechanism for reconciling the synchronisation of Heinrich events and Dansgaard-Oeschger cycles

Domain: earth_science
Published: 2024-04-05
Editorial outcome: major_revision

Abstract
The evolution of the northern hemispheric climate during the last glacial period was beset by quasi-episodic iceberg discharge events from the Laurentide ice sheet, known as Heinrich events (HEs). The paleo record places most HEs into the cold stadial of the Dansgaard-Oeschger cycle. However, not every Dansgaard-Oeschger cycle is associated with a HE, revealing a complex interplay between the two modes of glacial variability. Here, using a coupled ice sheet-solid earth model, we introduce a mechanism that explains the synchronicity of HEs and Dansgaard-Oeschger cycles. Unlike earlier studies, our mechanism does not require a trigger during the stadial. Instead, the atmospheric warming signal during the interstadial of the Dansgaard-Oeschger cycle causes enhanced ice stream thickening that leads to the HE during the late interstadial. We demonstrate that this mechanism reproduces the key HE characteristics and provides an explanation for synchronous HEs from different regions of the Laurentide ice sheet.

Introduction
Heinrich events are associated with massive iceberg discharge events from the Laurentide ice sheet - the ice mass covering the North American continent during past glacial periods. Subsequent melting of the icebergs left behind a distinct layer of ice-rafted debris on the North Atlantic ocean floor that ultimately led to the identification of Heinrich events1,2. During the Marine Isotope Stage 3 period (MIS3; ca. 60-25 ka BP), four Heinrich events have been identified in the paleo records2. Remarkably, most of them occurred during the cold phase of an Dansgaard-Oeschger cycle3, a period that is more commonly associated with ice mass gain rather than ice mass loss. Dansgaard-Oeschger cycles are characterised by a period of rapid, decadal warming of up to 14 °C in the high northern latitudes4, followed by a more gradual cooling spanning several centuries. Temperature reconstructions from ice cores indicate a dominant recurrence interval of ~1,500 years for Dansgaard-Oeschger cycles5. To this day, the underlying mechanisms interlinking Dansgaard-Oeschger cycles and Heinrich events remain poorly understood.

Over the last decades many trigger mechanisms for Heinrich events have been put forward in search of a better understanding of the governing processes as well as reproducing the characteristics of Heinrich events as derived from the paleo record6–10. Early theories focused on internally-driven oscillations6,7,11–13, suggesting that Heinrich events follow a two-stage pattern referred to as binge-purge cycle6. During the binge phase, the ice sheet becomes thicker causing an increase in basal temperatures until it crosses a critical thickness and temperature threshold, resulting in the onset of rapid basal sliding. The ensuing purge phase is characterised by elevated ice discharge that is maintained until the internal heat supply decreases and basal temperatures fall below the pressure-melting point. The validity of this mechanism has since been challenged because it struggles to explain the observed timing of Heinrich events within the Dansgaard-Oeschger cycle8,9. The occurrence of most Heinrich events during cold stadials, also referred to as phase locking7,14, suggests a common external climatic trigger. Motivated by observational evidence indicating an ocean subsurface warming signal prior to several Heinrich events15,16, more recent theories revolve around the ocean as key climate component for controlling the timing of Heinrich events8–10. Some studies suggest that the observed subsurface warming of the ocean led to the demise of an ice shelf wedged-in between Greenland and the Laurentide ice sheet8,9. The sudden disintegration of the ice shelf destabilised the Hudson ice stream, leading to an acceleration of the ice located upstream, very similar to the recently reported behaviour of regions in the Antarctic17. While this may explain the timing of Heinrich events within the Dansgaard-Oeschger cycle, it is incompatible with proxy observations, which provide no firm evidence of the existance of such a persistent ice shelf. This shortcoming inspired the development of a mechanism that links Heinrich events to the intrusion of warm subsurface ocean waters, as observed at marine outlet glaciers in Greenland and Antarctica10. The warm ocean water initiates widespread ice-sheet retreat through a positive feedback mechanism referred to as marine ice cliff instability. Subsequent glacio-isostatic uplift, as a result of the reduced ice mass, cuts off the path of the warm ocean waters to the ice stream, terminating the Heinrich event and leading to a regrowth of the ice sheet. The use of simple coupled models, however, has revealed that two weakly coupled autonomous oscillators, in this case the ice stream and the ocean, with different natural frequencies can synchronise with ice-ocean interactions alone, if the external forcing that both oscillators are exposed to is strong enough to overcome the natural frequency of the two oscillators14,18,19.

All previous theories have in common that they have been specifically tailored towards the glaciological and climatological setting of the Hudson ice stream7–10,19. Some studies7,10 have hinted, but without demonstrating it, that their mechanism could explain the observational evidence of synchronous Heinrich events from Hudson ice stream and other regions of the Laurentide ice sheet20. Naturally, synchronous Heinrich events from mechanisms based on the permanently marine-terminating Hudson ice stream are restricted to ice streams with such characteristics. However, the proxy record cannot dismiss the idea of synchronous Heinrich events from land-terminating ice streams that only become marine-terminating during the actual Heinrich event21. Support for this is provided by mapped geomorphological features from land-terminating ice streams of the Laurentide ice sheet that indicate surging22,23, even though the timing of these surges remains poorly constrained. A modern-day analogue for this hypothesis, albeit on a much smaller spatial scale than Heinrich events, comes from observations of synchronous surges from land-terminating mountain glaciers24. In the following, we propose a mechanism for Heinrich events that can trigger a pan-ice sheet response and is applicable to land and marine-terminating ice streams exhibiting distinctly different climatological settings. As a result, we are able to explain all of the dominant observed characteristics from the paleo record.

Results
Heinrich events in our simulations are preconditioned by the binge-purge mechanism6. However, the final trigger for the Heinrich events that determines the exact onset of the purge phase is provided by the Dansgaard-Oeschger cycle forcing. Conceptually, in our mechanism the rapid transition from stadial to interstadial results in a surface temperature warming that is accompanied by an elevated surface mass balance. Changes in surface temperature and surface mass balance instantaneously affect the geometric and thermal evolution of the ice stream. The warmer surface temperatures lead to increased ice deformation and enhanced frictional heat generation, while the higher surface mass balance flux increases ice stream thickening. Together these processes ensure that a critical ice thickness and basal temperature threshold, required for the onset of the internal instability, is reached toward the end of the interstadial phase (Fig. 1a, b). Because a Heinrich event is only triggered when the size of the region permitting rapid basal sliding crosses a certain threshold, the onset of the purge phase is typically delayed into the late interstadial of the Dansgaard-Oeschger cycle (Fig. 1c). Our canonical simulated Heinrich event lasts between 1000 and 1200 years and adds between 2.7 and 4 m to the global sea level with peak discharge rates reaching up to 45 mSv. In accordance with the paleo record, not every Dansgaard-Oeschger cycle is associated with a Heinrich event3, but, owing to the ice stream’s natural frequency, it takes between 3-5 Dansgaard-Oeschger cycles for the next Heinrich event to occur.

To demonstrate our mechanism we use simulations with a coupled ice sheet-solid earth model. Under constant MIS3 climate forcing, Heinrich events from the Hudson ice stream do not exhibit an inherent phase locking signal to Dansgaard-Oeschger cycles (henceforth referred to as reference simulation; see Methods: Initial state and Dansgaard-Oeschger cycle forcing). However, if we prescribe a 1500 year long Dansgaard-Oeschger cycle in addition to the constant MIS3 climate forcing, Heinrich events cluster in the late interstadial phase during a Dansgaard-Oeschger cycle (Fig. 2). This is illustrated by comparing an ensemble of simulations with prescribed and spatially homogeneous Dansgaard-Oeschger cycle forcings with the reference simulation (Supplementary Fig. 1). For the construction of our ensemble simulations mimicking a Dansgaard-Oeschger cycle forcing, we focus on two key parameters: (i) the surface temperature and (ii) the surface mass balance, which have previously been shown to affect the timing of Heinrich events25,26. The prescribed magnitudes range from 2 to 14° C for surface temperature (ΔTemp, see Methods: Dansgaard-Oeschger cycle forcing) and 0.009–0.1 m yr−1 for the surface mass balance (ΔSMB) and are motivated by reported values from ice core reconstructions4 as well as model simulations27. When we compare the empirical Heinrich event distributions of the forced simulations with the unforced reference simulation, we find that, based on the Kuiper’s test for two distinct empirical distributions (see Methods: Coupled simulations), 17 out of the 24 forced simulations show a Heinrich event distribution that is statistically significantly (R2 > 0.9) different from the reference simulation and exhibit phase locking (Supplementary Fig. 2). Only the smaller temperature and surface mass balance simulations do not show phase locking. The strength of the observed phase locking becomes stronger for larger surface mass balance perturbations which is in accordance with an earlier study26. Additional sensitivity simulations reveal that for the stronger perturbation scenarios (ΔTemp ≥ 10°C and ΔSMB ≥ 0.075 m yr−1) only one of the two atmospheric forcings is required to achieve phase locking (Supplementary Fig. 3). Hereby, perturbations to only the surface mass balance result in Heinrich events occuring deeper into the stadial phase of the Dansgaard-Oeschger cycle. This underlines that the timing of internally-driven Heinrich events is modulated by the prescribed atmospheric perturbations.

Since the paleo record indicates that Heinrich stadials tend to be longer than non-Heinrich stadials4,28, we performed a subset of the ensemble with a shorter 1000-year-long Dansgaard-Oeschger cycle, a 2000-year-long Dansgaard-Oeschger cycle, and experiments in which the Dansgaard-Oeschger cycle varies stochastically between 1000 and 2000 years in 100 year increments. Our simulations show that longer Dansgaard-Oeschger cycles result in enhanced phase locking and push the onset of the Heinrich events deeper into the stadial of the Dansgaard-Oeschger cycle (~135 years later, Supplementary Fig. 4). In contrast, the shorter 1000 year long Dansgaard-Oeschger cycle shows weaker phase locking (Supplementary Fig. 4). The response to an enhanced noise level in the Dansgaard-Oeschger cycle length is that the phase locking becomes weakened (Supplementary Fig. 5). Despite the Heinrich events being more spread out, the vast majority of the events still occurs during the cooling phase of the Dansgaard-Oeschger cycle (Supplementary Fig. 5). Many events show even an onset much deeper into the stadial phase which is more in accordance with reconstructions. Regardless of the length of the Dansgaard-Oeschger cycle, phase locking is persistent in almost all our simulations. This is corroborated by additional simulations with a more realistic, spatially-varying forcing (Supplementary Fig. 6). These simulations also display a distinct phase locking signal (Supplementary Fig. 7).

The main characteristic of Heinrich events from the paleo record that previously proposed mechanisms7–10 have failed to demonstrate is the phenomenom of synchronous Heinrich events from the Hudson ice stream and other regions of the Laurentide ice sheet. Paleo ice stream reconstructions indicate that other ice streams may have been active during a number of the recorded Heinrich events20, such as Mackenzie ice stream, which is located in the north-west of the Laurentide ice sheet. Here, we show that our proposed Heinrich mechanism also leads to phase locking for Mackenzie ice stream, albeit weaker than for Hudson ice stream, resulting in synchronous Heinrich events from both ice streams. Mackenzie ice stream is located in a distinctly different glaciological and climatological setting in comparison to Hudson ice stream, but also exhibits Heinrich-like events in our simulations. Therefore, it presents an ideal case to test the generalisation of our mechanism. The surface temperatures are ~8–10°C warmer than in the Hudson area, resulting in shorter Heinrich event intervals. The timing of the events is primarily controlled by surface temperature rather than surface mass balance26. During the binge phase, Mackenzie ice stream terminates on land and only becomes marine-terminating for short periods during the purge phase. Due to the absence of a well-confined subglacial trough, the periodicity of Heinrich events for Mackenzie ice stream shows an enhanced variability. These circumstances result in a noisier distribution of Heinrich events with weaker phase locking throughout our experiments in comparison to Hudson ice stream (Supplementary Fig. 8). However, phase locking becomes stronger for larger surface temperature perturbations. Similar to the Hudson ice stream, longer Dansgaard-Oeschger cycles also result in stronger phase locking and push the occurence of Heinrich events into the stadial phase of the Dansgaard-Oeschger cycle (Supplementary Fig. 9).

To illustrate that our proposed mechanism supports synchronous Heinrich events from Hudson and Mackenzie ice streams, we use our simulation ensemble with spatially homogeneous Dansgaard-Oeschger perturbations and compute the likelihood of a Heinrich event from Mackenzie ice stream within a time window of 300 and 600 years of a Heinrich event from Hudson ice stream. To improve the signal-to-noise-ratio, we average the likelihood of synchronous Heinrich events across the surface mass balance ranges. This step is justified because the periodicity of Heinrich events from Mackenzie ice stream has shown to be more sensitive to surface temperature changes than surface mass balance changes26. The results show that the likelihood of synchronous Heinrich events substantially increases for higher surface temperature perturbations from ~8% to ~17% and ~8% to ~22% for the 300 year time window and for the 600 year window, respectively (Fig. 3). This corresponds to a total of four more synchronous Heinrich events for the 300 year window and six more synchronous Heinrich events for the 600 year window. An even stronger signal for synchronous Heinrich events (up to 50%) is observed for the longer Dansgaard-Oeschger cycle experiments, supporting the result that longer Dansgaard-Oeschger cycles enhance phase locking. As a consequence, every fifth Heinrich event for the 1500-year-long Dansgaard-Oeschger cycle and every second Heinrich event for the 2000-year-long Dansgaard-Oeschger cycle from the Hudson ice stream is accompanied by a Heinrich event from the Mackenzie ice stream. This matches well with reconstructions from ice-rafted debris that show contributions from source regions other than the Hudson area for some, but not all recorded Heinrich events20,29.

Here, we have introduced a Heinrich event mechanism in which internally driven ice-sheet oscillations are paced by an atmospheric perturbation signal caused by Dansgaard-Oeschger cycles. Unlike earlier theories, our mechanism does not need a trigger event during the stadial phase of the Dansgaard-Oeschger cycle. This mechanism successfully reproduces all of the main characteristics of Heinrich events as evidenced by the paleo record. We demonstrate that our mechanism not only lead to phase locking of Heinrich events in the late interstadial phase of Dansgaard-Oeschger cycles, but can also produce synchronous Heinrich events from the Hudson and Mackenzie ice streams. This simulated behaviour is robust under a wide range of forcing scenarios. A distinct advantage of our mechanism is that it is not restricted to marine-terminating ice streams, but applies equally to land-terminating ice streams that only become marine-terminating during the actual Heinrich event. Moreover, our findings indicate that the length of activation and stagnation cycles for present-day ice streams, which are important contributors to sea-level rise, may be controlled by atmospheric perturbations.

Methods
In our simulations, we employ the coupled ice sheet-solid Earth model mPISM-VILMA, consisting of the modified Parallel Ice Sheet Model (mPISM30,31) and the global VIscoelastic Lithosphere and MAntle model (VILMA32). The ice dynamics in mPISM are based on the superposition of the shallow ice approximation (SIA) and shallow shelf approximation (SSA), making it suitable to model slow ice flow, fast ice flow as well as grounding-line motion. The model has previously been shown to successfully reproduce the dynamics of internally-driven Heinrich events26,31. We use a northern hemispheric setup of mPISM with a horizontal mesh resolution of 10 km. The three-dimensional temperature distribution of the ice sheet is computed based on the enthalpy method33. In mPISM, we modified the enthalpy advection to use the c-grid that is also used in PISM for mass advection.

The surge characteristics are strongly affected by the basal sliding parameterisation. In all our simulations, we apply a non-linear Weertmann-type friction law that relates the basal shear stress τb to the basal sliding velocity ub of the form:

τb = (τc / u0^q) * ub^q : (1)

Here q is the basal sliding coefficient, which is set to 0.25 in all simulations, τc is the yield stress, and u0 is the threshold velocity at which the magnitude of the basal shear stress equals τc. We set u0 to 70 m yr−1 for all simulations. The yield stress τc is computed via the Mohr-Coulomb criterion

τc = tan(ϕ) * Ntill : (2)

Here, ϕ represents the till friction angle that is parameterised based on bedrock elevation. The friction angle varies linearly from 15° to 30° for bedrock elevations between -300 m and 400 m and remains constant at the upper or lower limit otherwise. The effective pressure (Ntill) is determined from the ice overburden pressure and the till saturation. Subglacial till water is simulated based on a non-conserving undrained plastic bed model34. This means that meltwater is produced and consumed locally, and horizontal transport of meltwater is neglected. To enhance basal sliding for the ice-sheet surges, we spread 50% of the basal heating effect from the sliding in each cell over the four neighbour cells. This results in a faster surge propagation and shorter surge periods that is more in accordance with reconstructions.

The formation of icebergs is approximated using two separate parameterisations. First, we apply so-called Eigencalving35, where the detachment of icebergs is based on the spreading rates of the ice shelf. In addition, we also apply a thickness calving criterion which removes ice at the calving front that is thinner than 100 m at a maximum rate of one grid cell per time step.

The evolution of the solid Earth is modelled by VILMA using a 1D Earth structure. This implicitly assumes that the structure of the solid Earth varies in depth but is otherwise spatially homogeneous. The coupling of mPISM and VILMA is performed every 100 model years and includes the exchange of the ice thickness distribution and glacial isostatic adjustment fields between the two models. As a global model, VILMA requires a global ice thickness distribution as input. Since we focus on the evolution of the Laurentide ice sheet, we keep the ice thickness distribution of all other northern and southern hemispheric ice sheets constant through time.

The initial state for all our mPISM-VILMA simulations was taken at 36 kyr before present from an accelerated Marine Isotope State 3 (MIS3) simulation with the Max Planck Institute for Meteorology Earth System Model (MPI-ESM) coupled to mPISM and VILMA. The state represents MIS3 climate conditions and is characterised by extensive ice sheets such as the Laurentide, Fennoscandian, and Russian ice sheets in the northern hemisphere and the Antarctic ice sheet in the southern hemisphere.

All simulations are forced with a time-constant MIS3 background climate. The climate forcing includes ocean temperature and ocean salinity to compute basal melt rates under ice shelves. The atmospheric forcing, consisting of ice surface temperature and surface mass balance fields, is provided through an energy balance model36. The surface temperature hereby corresponds to the temperature of the lowest layer in a snow model, restricting ice surface temperature to values below freezing. Applying constant MIS3 forcing ensures that the model drift and corresponding ice volume differences within our simulations remain low (<3%). To superimpose the spatially homogeneous rapid warming and gradual cooling behaviour of Dansgaard-Oeschger cycles as derived from the paleo record, we constructed a forcing function that mimics this behaviour of the following form:

forcing = (1/d) * arctan(d * sin(2π * (t - ti) / T) / (1 - d * cos(2π * (t - ti) / T))) : (3)

Here d represents the distortion factor and is set to 0.75 in all our simulations to reproduce the rapid warming-gradual cooling cycle. T is the length of the Dansgaard-Oeschger cycle. Depending on the simulation, T takes values of 1000, 1500, or 2000 years. The variables ti and t represent the current and the initial time of the model simulation. To ensure the function values of Eq. (3) map between [−1, 1], we scale the weights according to:

forcing_scaled = forcing / max(forcing) : (4)

It is important to note that the weights of the forcing function sum to zero when integrated over an entire Dansgaard-Oeschger cycle. This means that the background climate is not altered during the simulations and only on time scales shorter than the length of the Dansgaard-Oeschger cycle does the forcing function result in a transient perturbation of the climate forcing. The final forcing functions for the surface temperature and the surface mass balance are described by:

ΔSMB = ΔSMB * forcing_scaled : (5)

ΔTemp = ΔTemp * forcing_scaled : (6)

where Δ(SMB, Temp) are the prescribed anomaly magnitudes of the respective simulation.

All presented perturbation simulations were integrated for a total of 57,000 years, equivalent to 38 Dansgaard-Oeschger cycles. This resulted in a minimum of 9 Heinrich events for Hudson ice stream and 11 Heinrich events for Mackenzie ice stream. As Mackenzie ice stream is located in a climate close to a threshold under which Heinrich events cannot be maintained anymore26. A few simulations exist in which Mackenzie ice stream turns into a persistent ice stream. However, since Mackenzie ice stream predominatly serves the purpose of showing the possibility of synchronous Heinrich events, we deem this less critical. To permit a more robust statistical analysis, the reference simulation was integrated for a total of 123,000 years, corresponding to 82 Dansgaard-Oeschger cycles. To determine whether Heinrich event distributions from the perturbation simulations are significantly different from the Heinrich event distribution of the reference simulation, we used the Kuiper’s test. This test is closely related to the more widely-known Kolmogorov–Smirnov test, but has the additional property that it remains invariant under cyclic transformation, making it well-suited to the application of cyclic variations such as our Dansgaard-Oeschger forcing.

To evaluate the robustness of our results with the idealised forcing, we ran additional experiments with a more realistic and spatially-varying surface temperature and surface mass balance Dansgaard-Oeschger cycle forcing. The forcing is constructed from a MIS3 simulation with a comprehensive Earth System Model that includes interactive ice sheet and solid earth components (see “Methods” section: Initial state). As this model does not exhibit the characteristic Dansgaard-Oeschger cycles of MIS3, we approximate the expected Dansgaard-Oeschger cycle atmospheric forcing as the difference between years with strong Atlantic Meriodional Overturning Circulation (AMOC; >14 Sv) and years with weak AMOC (<12 Sv). Overall, the AMOC varied between ~8–18 Sv. Because we still anticipate the forcing magnitudes to be lower than from a full Dansgaard-Oeschger cycle, we perform simulations where we multiply the forcing magnitudes by constants of 1, 2, 3 and 4. The constructed 2D forcing shows similar perturbation magnitudes to the range covered in the spatially homogeneous simulations (Supplementary Fig. 6). We therefore observe a phase locking behaviour of Hudson ice stream that closely resembles the results of the simulations with the more idealised forcing (Supplementary Fig. 7).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study suggests a new hypothesis for Heinrich Events—mass iceberg evacuations from the\nLaurentide ice sheet. The cause of Heinrich Events have long been debated because they occur during\nthe coldest portions of the glacial cycle. The original hypothesis for Heinrich Events, proposed by Doug\nMacAyeal involves a two-stage process called the binge-purge cycle where the ice is first frozen to its\nbed, but then reaches the pressure melting point at the base. The transition to the melting point at\nthe base results in a surging behavior of the ice sheet. Although this hypothesis resulted in a time\nscale for Heinrich Events roughly comparable to that observed for Heinrich Events, the predicted\nHeinrich Events were not synchronized with the DO cycles. Subsequent hypotheses have largely\nrevolved around the fact that paleo-proxies indicate that ocean warming associated with DO events\nprecedes Heinrich Events and these subsequent studies have assumed that the warming ocean can act\nas a trigger for Heinrich Events.\nThe authors show here that the original binge-purge hypothesis can be resuscitated if temperature\nand snowfall are taken to be a function of DO events. In the authors model, this results in a Heinrich\nevent cycle that is phase locked to the DO events. The phase locking proposed here seems similar to\nthe phase locking found by Mann et al., (2021), although that study assumed an ocean trigger.\nMoreover, the novelty of this is the syncing of the binge-purge cycle with the DO events. But the\nclassic binge-purge hypothesis proposed seems like it is still at the heart of the mechanism proposed.\nOverall, I think the manuscript is generally well written and put together. I include a couple of\nsuggestions and questions below.\n1. Some questions about physics and mechanisms. What causes iceberg calving in the current model?\nTo be more precise, during the purge phase of the DO cycles the ice sheet bed reaches the pressure\nmelting point. This causes the ice stream to surge. But this could cause the ice stream to thin and\nadvance, perhaps forming an ice shelf? Presumably the same mechanism would work on land-\nterminating ice streams, but would not produce an armada of icebergs? Fundamentally, what causes\nthe ice to break and the deluge of calving events? The surging part of the model is clear to me, but\nsurging glaciers observed today don’t typically create armadas of icebergs even when they terminate\nin the ocean. Figure 2 looks like a small ice shelf is forming, but what controls the extent of the ice\nshelf?\n2. A related question is how do the authors treat sliding, till deformation and the thermodynamics of\nthe bed? Sliding, heat flow and the dynamics at the ice-till interface seems crucial to the model and it\nwould be helpful to see this described in the methods section with a little bit more detail.\n3. I don’t think it is entirely true that other mechanisms proposed aren’t able to reproduce behavior of\nthe non-Hudson ice streams. The ocean forcing hypothesis suggests widespread ocean warming which\nwould trigger widespread retreat. The retreat of Greenland and other glaciers driven by ocean\ntemperatures is analogous to what is observed currently. Where it seems the ocean forced model will\nresult in significant differences is that the model proposed by the authors will apply to both terrestrial\nand marine margins. This is hinted at in the paper, but not elaborated on. I would encourage the\nauthors to take readers through the evidence that land-terminating ice streams also surged at the\nexact same time as marine-terminating ice streams. I think this would be powerful evidence that it\nisn’t an ocean signal that is driving these events and would provide strong support for the authors\nmodified binge-purge model.\n4. My next comment relates to how the authors have driven their model. The DO events observed in\nthe record are not periodic. There have been some hypothesis that DO events are stochastic poisson\nprocesses. What I would like to know is whether the phase locking is a consequence of the assumption\nthat the DO events are periodic and occur quasi-regularly. What happens if the onset of DO events are\nset according to the observed records (or stochastically) with long gaps between some of the DO\nevents. Do you still avoid Heinrich Events during the long gaps? Does the phase locking only occur in a\nsmall portion of the record? Or does it occur across large portions of the record?\n5. This is a minor point, but why not use the temperature and mass balance perturbation associated\nwith DO events inferred from the Greenland ice sheet ice cores? This seems like it would be especially\nuseful in a simulation driven with the timing and magnitude of observed DO events. I realize that this\nwould not necessarily apply to the Laurentide, but it would provide a magnitude.\nMann, L. E., Robel, A. A., & Meyer, C. R. (2021). Synchronization of\nHeinrich and Dansgaard-Oeschger events through ice-ocean interactions. Paleoceanography and\nPaleoclimatology", "reviewer_2": "This manuscript on “A mechanism for reconciling the synchronization of Heinrich events and\nDansgaard Oeschger cycles” discussed a new modeling-based theory for millennial-scale glacial\nclimate variability and the interactions between two different modes of variability. The central\nargument is that atmospheric processes at the surface of the Laurentide ice sheet allow for\nsynchronization between internally generated ice sheet variability (Heinrich events) and ocean\nvariability (D-O events). This synchronization potentially helps resolve several outstanding questions\nin the literature including the occurrence of Heinrich events during D-O stadials and the activation of\nmany Laurentide ice streams during some Heinrich events.\nThe particular idea here (synchronization of Heinrich and DO events through surface temperature\nand SMB) is novel, though it has unacknowledged similarities to previous work on Heinrich-DO\nsynchronization. The manuscript does a good job explaining why these mechanisms explain certain\nenigmas within the paleoclimate records. It does so by building on well established models and ideas\nfrom climate and ice sheet dynamics.\nMajor concerns:\n1. This manuscript argues that certain aspects of the results improve on prior theories in ways which\nare not clearly supported by the model results.\n(a) In the abstract and later, a distinction is made between a \"trigger mechanism\" and\n\"preoconditioning\". It is not clear to me how these are different. It could be argued that the\natmopsheric warming \"triggers\" the Heinrich event by causing the basal thermal state to warm. To\nme, this difference seems to be more rhetorical than substantive. A related point is that in systems of\nsynchronized oscillators causality isn’t straightforward. The oscillations in the two systems (ice sheet\nand ocean in this case) would happen regardless of the interaction between them, but both systems\ninteract with each other to produce a situation in which the timing of the events is consistent. One\ndoes not \"cause\" the other, but their timing is related.\n(b) In the abstract and elsewhere, it is argued that the mechanism put forward in this study \"provides\na simple explanation for the observational evidence of synchronous Heinrich events from different ice\nsheets\". The MacKenzie and Hudson Strait ice streams are in the same ice sheet, so its no clear that\nthis is in fact demonstrated in this study. Given that the model configuration in question seemingly\nincludes the Fennoscandian ice sheet, it would seem that a similar analysis could be done to determine\nwhether the atmospheric mechanism also produces synchronization with Fennoscandian ice streams.\nHowever, if such evidence is lacking, then it should not be claimed. Perhaps more accurate would be to\nsay that it \"provides a simple explanation for the observational evidence of synchronous discharge\nevents from different ice streams within the Laurentide Ice Sheet\".\n2. The central idea in this study is the most similar to a recent study by Mann et al. (2021 in\nPaleoclimatology and Paleoceangoraphy) which is also able to explain the synchronization of Heinrich\nand DO events, but instead using ice-ocean interactions to generate synchronization between a flow\nline glacier model and a simple ocean model. Clearly this manuscript proposes a different\nsynchronization method and uses more realistic models, but that prior study provides many potential\npoints of comparison (including more comprehensive parameter sweeps). One thing that is not clear is\nwhether the Mann et al. ocean-based mechanism for synchronization, if included in a more complex\nmodel like this one, could also explain all the aspects of the paleoclimate record that the authors focus\non in this study. The paper would benefit from a more extensive citation and comparison to the\nprevious study by Mann et al.\n3. One surprising result is that surface temperature perturbations from DO events have any signficant\nimpact on ice sheet sliding on time scales of a few thousand years. The diffusion time scale of heat\nthrough thick ice is tens of thousands of years. This is why in the original papers on the B-P theory for\nHeinrich events, MacAyeal discounts the role of surface temperature changes in driving Heinrich\nevents because they would take too long to diffuse through the ice sheet. It would seem that a more\nrobust explanation is needed here, particularly with some plots showing the influence of adding the\natmospheric perturbation on the basal heat budget of ice streams.\nMinor comments:\nLine 87: In what sense do you mean \"statistically significant different\"? What metrics characterizing\nthe Heinrich event are you measuring and how is that difference measured?\nLine 90: Why does it seem that the SMB perturbation has more of an effect than the temperature\nperturbation?\nLine 93: Increased surface temperature will increase snow accumulation at high elevation but also\nincrease surface melting at low elevations? Is the latter not important at the edges of the ice sheet?\nLine 95: What is the mechanism explaining why \"The onset of the Heinrich event is, however, delayed\ninto the early phase of the transition into the stadial of the Dansgaard-Oeschger cycle\"\nLine 99: Similarly, can you explain why it takes 3-5 D-O cycles before the next Heinrich events occurs?\nIs it simply that H events are paced internally by ice streams? Do the Hudson Strait and MacKenzie\nice streams have the same natural periodicity for discharge events?\nLine 140: Can you say with significance that an increase from 8% to 17% or to 22% would not just\nhappen by chance? How many more Heinrich events is this? Is it more than one?\nLine 158: What makes this an advantage? Why do you need this mechanism to apply to land-\nterminating glaciers? It wouldn't help explain anything about the marine sediment record which is only\ninformed by marine-terminating glacier discharge.\nL163: related to point #3 above, it would be helpful to explain here how internal temperature and\nSMB are modeled in mPISM-VILMA since those processes are central to the mechanism in this study.", "reviewer_3": "Review of the manuscript \"A mechanism for reconciling the synchronisation of Heinrich events and\nDansgaard-Oeschger cycles\" by Dr Schannwell and colleagues.\nThis paper presents a novel mechanism that is introduced to explain the synchronous occurrence of\nHeinrich events and the stadial phase of Dansgaard-Oeschger cycles during the last glacial period.\nHeinrich events, quasi-episodic iceberg discharge events from the Laurentide ice sheet into the North\nAtlantic, have been a significant feature of the northern hemispheric climate during the last ice age\ncycle and is a difficult process to study because of the large range of temporal and and spatial scales\nthat are involved and the small number of Heinrich events observed. Using simulations with a coupled\nice sheet-solid earth model, the researchers present a novel mechanism that eliminates the need for a\ntrigger mechanism during the stadial phase. Instead, it is proposed that the interstadial phase of the\nDansgaard-Oeschger cycle, characterised by atmospheric warming, preconditions the ice sheet for a\nHeinrich event, causing it to occur during the subsequent stadial phase.\nThis paper re-visits the binge-purge hypothesis and presents a fresh re-interpretation on the\nmechanisms behind Heinrich events. While I personally resonate with the subsurface ocean warming\nmarine ice sheet instability hypothesis which occurs during stadials, it is important to provide other\npossible theories that may help explain a dominant feature of ice age climate. I feel the paper is well\nwritten and the study rests on some solid previous work by the group. I would like to see some of my\nconcerns addressed which I list below.\nGeneral Structure of the paper:\nI like Figure 2 and I think this nice graphic should be the highlight of the paper. I would suggest\nmoving this to Figure 1. It should not be too hard to rework the text on pages 4 and 5. Start with Pg 5\nL92. I have more specific comments on the figure below.\nSpecific comments:\nFigure 1:\nWhy is there a colour wheel in the “Unforced reference simulation” in figure 1a . Is this not forced with\nconstant MIS3 forcing and so should be blank? Or have I missed something?\nThe current Figure 1 could be enhanced with another schematic diagram beside the top panel\n“Unforced reference simulation” (to the left , illustrating how the ice volume builds as a curve spirals\nand then collapses with a Heinrich event. Also it was not initially clear what all the thin black lines\nwere indicating (the ensemble of runs) and how many revolutions they made. Do all the thin lines\nstart at 0 time at some small ice volume? It will take the reader some time to see what is happening\nso a simpler schematic of one simulation with text (possibly) highlighting HEs and DO stadial and\nInterstadial periods.\nThe matrix of plots in figure 1 has \\Delta SMB vs delta \\Delta Temp, shouldn’t this be A_smb vs\nA_temp (the amplitudes given from the equations (3 and 4) in the Methods section. This could be\nchanged and pointed out in the figure caption.\nWhat does the white significance level indicate (< 90th percentile) , maybe state this.\nFigure 2:\nIt would be useful to highlight that the insets to figure 2 are from equation 1 in the Methods section\nand the colours on the curves represent the derivative of the temperature and mass balance changes.\nTo most people I think they would associate the point in Figure 2c as early interstadial. It is so close to\nthe D-O warming event that I don’t think one can possibly argue that this is a stadial\ncondition.\nPg. 3 L 51-54: Is it possible to infer evidence of ice shelves from the DeVernal et al 2000 paper on\nLGM. Is there a strong statement in this paper about no ice shelves existing along the Hudson Strait\nIce Stream before H1?\nPg 4: L67-68: On explaining the observed characteristics of the paleo record:\nThe atmospheric forcing functon(Equation 1 Methods does not have the same quick drop off at the end\nof the interstadial as seen in the ice core records (e.g. Rassmussen et al 2014). In particular if you\nlook at GI-13 the interstadial period remains elevated for some time (800 years) and then quickly\ndrops off. The results here possibly work here with H5, but if you look at GI-9 the interstadial is short\n(< 500 years). The interstadials previous to the other Heinrich events do no fit well into the context of\nthis paper and may have other origins other than the Hudson’s Bay Ice Stream. It would be difficult to\nsay anything using observations to corroborate the findings of this study.\nPg 4. L80: This would be a good place to highlight the amplitude from the methods (A_sbm and\nA_temp) and refer to the methods.\nPg5 L88-89:\nIf you assume that there are Heinrich events at 30 , 39 and 47.5 kaBP you would argue that there is a\nseries of decreasing interstadial lengths previous to the H-Events as observed from the ice core data\n(a sort of damped D-O oscillation).\nI would like to see what kind of effect would result from a damped harmonic oscillator (e.g if you\nadded an exp(-t/tau) in front of the arctan term in Methods equation 1. From the results presented ,\nI assume that the synchronisation would decrease as the amplitude goes down. So I have a hard time\nseeing how this could be corroborated with observations. I didn’t see any strong attempt to argue how\nthe model results compare with observations of D-O variability.\nZieman et al 2019 (CP) include oceanic forcing. Without the influence of ocean effects it is really hard\nto argue that this is “the” mechanism that causes H-Events.\nPg5 L96-97:\nIs this consistent with Zieman et al 2019? How much is this in mSv (50)?\nPg5 102-104: These results suggest that the short D-O cycles would have an H-Event close to the end\nof the interstadial, is this consistent with observations.\nPg 6: L 122-123: Figure 1 suggests that the SMB is the primary driver from figure 1 if I read the\nsignificance correctly.\nPg 14: L219-224: What is the range of the AMOC in the MIS3 simulation?\nIn the extended data figures, what is the change in radial ice volume in meters of RSL (it is only given\nin m^3."}, "author_response": "Response to the reviewers\nWe thank all referees for their thoughtful and thorough reviews of our paper. We appreciate you\ntaking the time to complete these reviews and welcome your helpful comments. In the following\nwe address their concerns point by point. Throughout this response to review document referee\ncomments are provided in regular, non- italic font text, our response comments are provided in red\nfont.\nTo accomodate the main comments of the reviewers, we have made the following main changes to\nthe manuscript:\n• We have substantially expanded the model description section to provide a more complete\npicture of the implementation of key processes.\n• We have added four more simulations in which we vary the DO cycle length stochastically\nfrom one cycle to the other.\n• Reworked the general structure of the manuscript by swapping the order of Figure 1 and\nFigure 2 as well as text adjustments reflecting the modified Figure order.\nReviewer 1\nReviewer Point P1.1 — This study suggests a new hypothesis for Heinrich Events—mass\niceberg evacuations from the Laurentide ice sheet. The cause of Heinrich Events have long been\ndebated because they occur during the coldest portions of the glacial cycle. The original hypothesis\nfor Heinrich Events, proposed by Doug MacAyeal involves a two-stage process called the binge-\npurgecyclewheretheiceisfirstfrozentoitsbed, butthenreachesthepressuremeltingpointatthe\nbase. The transition to the melting point at the base results in a surging behavior of the ice sheet.\nAlthough this hypothesis resulted in a time scale for Heinrich Events roughly comparable to that\nobserved for Heinrich Events, the predicted Heinrich Events were not synchronized with the DO\ncycles. Subsequent hypotheses have largely revolved around the fact that paleo-proxies indicate\nthat ocean warming associated with DO events precedes Heinrich Events and these subsequent\nstudies have assumed that the warming ocean can act as a trigger for Heinrich Events.\nThe authors show here that the original binge-purge hypothesis can be resuscitated if tempera-\nture and snowfall are taken to be a function of DO events. In the authors model, this results in a\nHeinrich event cycle that is phase locked to the DO events. The phase locking proposed here seems\nsimilar to the phase locking found by Mann et al., (2021), although that study assumed an ocean\ntrigger. Moreover, the novelty of this is the syncing of the binge-purge cycle with the DO events.\nBut the classic binge-purge hypothesis proposed seems like it is still at the heart of the mechanism\nproposed.\nOverall, I think the manuscript is generally well written and put together. I include a couple\nof suggestions and questions below.\nReply: We thank the reviewer for the positive assessment of our manuscript.\ncalving in the current model? To be more precise, during the purge phase of the DO cycles the\nice sheet bed reaches the pressure melting point. This causes the ice stream to surge. But this\ncould cause the ice stream to thin and advance, perhaps forming an ice shelf? Presumably the\nsame mechanism would work on land-terminating ice streams, but would not produce an armada\nof icebergs? Fundamentally, what causes the ice to break and the deluge of calving events? The\nsurging part of the model is clear to me, but surging glaciers observed today don’t typically create\narmadas of icebergs even when they terminate in the ocean. Figure 2 looks like a small ice shelf is\nforming, but what controls the extent of the ice shelf?\nReply: We have added a description of how iceberg calving is parameterised in our simulations. The\nreviewer is correct that this mechanism works also for land-terminating ice streams, but they would\nnot produce icebergs but a liquid freshwater flux that would be transported to the ocean via rivers.\nDuring the build-up phase, Mackenzie is land-terminating. For the Hudson ice stream, a small ice\nshelf typically forms during a Heinrich event, but its existence is short-lived and the ice-stream front\nretreats to its pre-surge position. This is also shown in this video from our previous paper: Video link.\nIn comparison, Mackenzie ice stream is land-terminating during the build-up phase and only becomes\nmarine terminating with a small floating ice tongue during the actual Heinrich event.\nReviewer Point P1.3 — A related question is how do the authors treat sliding, till deformation\nand the thermodynamics of the bed? Sliding, heat flow and the dynamics at the ice-till interface\nseems crucial to the model and it would be helpful to see this described in the methods section\nwith a little bit more detail.\nReply: Weagreethatthisstronglyinfluencesthesurgecharacteristics. Wehavesubstantiallyexpanded\nthe Methods section to describe this is more detail.\nReviewer Point P1.4 — I don’t think it is entirely true that other mechanisms proposed aren’t\nable to reproduce behavior of the non-Hudson ice streams. The ocean forcing hypothesis suggests\nwidespread ocean warming which would trigger widespread retreat. The retreat of Greenland and\nother glaciers driven by ocean temperatures is analogous to what is observed currently. Where\nit seems the ocean forced model will result in significant differences is that the model proposed\nby the authors will apply to both terrestrial and marine margins. This is hinted at in the paper,\nbut not elaborated on. I would encourage the authors to take readers through the evidence that\nland-terminating ice streams also surged at the exact same time as marine-terminating ice streams.\nI think this would be powerful evidence that it isn’t an ocean signal that is driving these events\nand would provide strong support for the authors modified binge-purge model.\nReply: We have expanded the paragraph why we think it is an advantage that our mechanism works\nequally well for permanently marine-terminating ice stream, such as Hudson ice stream, as well as for\nice streams that are land-terminating and only become marine-terminating during the actual Heinrich\nevent. We have added evidence of surging from land-terminating ice streams of the Laurentide ice sheet\nas well as evidence from synchronous surges of mountain glaciers, a modern-day analouge for Heinrich\nevents.\nReviewer Point P1.5 — My next comment relates to how the authors have driven their model.\nThe DO events observed in the record are not periodic. There have been some hypothesis that DO\nevents are stochastic poisson processes. What I would like to know is whether the phase locking is\na consequence of the assumption that the DO events are periodic and occur quasi-regularly. What\nhappens if the onset of DO events are set according to the observed records (or stochastically) with\nlong gaps between some of the DO events. Do you still avoid Heinrich Events during the long gaps?\nDoes the phase locking only occur in a small portion of the record? Or does it occur across large\nportions of the record?\nReply: Our analysis here focuses primarily on the MIS3 period during which five of the six recorded\nHeinrich events occurred. During this period, DO cycles are fairly regular and prolonged gaps between\nDO cycles are not present. If gaps would extend to longer than 5,000 - 7,000 years, our model would\nproduce a Heinrich event even in the absence of a DO cycle. To further investigate the robustness of\nour results to stochastic variability in the DO cycle, we performed four additional simulations in which\nthe DO cycle length is randomly varied between 1,000 and 2,000 years (in 100 year intervals). These\nsimulationsshowanexpectedweakenedphaselockingsignalinresponsetotheincreasednoiselevel(see\nExtendedFigure8intherevisedmanuscript). Still, thevastmajorityoftheeventsstilloccursduringthe\ncooling phase of the DO cycle. Interestingly, the onset of many Heinrich events in theses simulations\nhappens much deeper into the stadial phase which is more in accordance with oberservations. It is\nimportant to note that not all of the recorded Heinrich events occurred during the stadial phase of DO\ncycles.\nReviewer Point P1.6 — This is a minor point, but why not use the temperature and mass\nbalance perturbation associated with DO events inferred from the Greenland ice sheet ice cores?\nThis seems like it would be especially useful in a simulation driven with the timing and magnitude\nof observed DO events. I realize that this would not necessarily apply to the Laurentide, but it\nwould provide a magnitude.\nReply: The chosen magnitude ranges are based on Greenland ice core reconstructions as mentioned in\nLine 98 in our initial submission. Because there are substantial variations in signal magnitude from DO\ncycle to DO cycle, we chose to rather encompass this range with our ensemble approach rather than\nchoose a sepcific value. As pointed out by the reviewer and shown in our 2D model forcing (Extended\nFigure 6) in the initial submission, the magnitudes over the Laurentide ice sheet are much lower than\nwhat has been recorded in the Greenland ice cores.\nReviewer 2\nReviewer Point P2.1 — This manuscript on “A mechanism for reconciling the synchroniza-\ntion of Heinrich events and Dansgaard Oeschger cycles” discussed a new modeling-based theory\nfor millennial-scale glacial climate variability and the interactions between two different modes of\nvariability. The central argument is that atmospheric processes at the surface of the Laurentide ice\nsheet allow for synchronization between internally generated ice sheet variability (Heinrich events)\nandoceanvariability(D-Oevents). Thissynchronizationpotentiallyhelpsresolveseveraloutstand-\ning questions in the literature including the occurrence of Heinrich events during D-O stadials and\nthe activation of many Laurentide ice streams during some Heinrich events.\nThe particular idea here (synchronization of Heinrich and DO events through surface temperature\nand SMB) is novel, though it has unacknowledged similarities to previous work on Heinrich-DO\nsynchronization. The manuscript does a good job explaining why these mechanisms explain certain\nenigmas within the paleoclimate records. It does so by building on well established models and\nideas from climate and ice sheet dynamics.\nReply: We thank the reviewer for the positive assessment of our manuscript.\nReviewer Point P2.2 — This manuscript argues that certain aspects of the results improve on\nprior theories in ways which are not clearly supported by the model results.\n(a) In the abstract and later, a distinction is made between a ”trigger mechanism” and ”preocon-\nditioning”. It is not clear to me how these are different. It could be argued that the atmopsheric\nwarming ”triggers” the Heinrich event by causing the basal thermal state to warm. To me, this\ndifference seems to be more rhetorical than substantive. A related point is that in systems of syn-\nchronized oscillators causality isn’t straightforward. The oscillations in the two systems (ice sheet\nand ocean in this case) would happen regardless of the interaction between them, but both systems\ninteract with each other to produce a situation in which the timing of the events is consistent. One\ndoes not ”cause” the other, but their timing is related.\nReply: We have rephrased this in the revised version of the manuscript. We now highlight that the\nunderlying binge-purge mechanism preconditions the ice stream for a Heinrich event. The final trigger\nthat determines the onset of the rapid sliding is however provided by the DO cycle forcing. As a result,\nthis means that Heinrich events occur in our model even without the DO cycle forcing, but the timing\nand periodicity of the events differs between the unforced and forced Heinrich event simulations.\nReviewer Point P2.3 — (b) In the abstract and elsewhere, it is argued that the mechanism put\nforward in this study ”provides a simple explanation for the observational evidence of synchronous\nHeinrich events from different ice sheets”. The MacKenzie and Hudson Strait ice streams are in\nthe same ice sheet, so its no clear that this is in fact demonstrated in this study. Given that the\nmodel configuration in question seemingly includes the Fennoscandian ice sheet, it would seem that\na similar analysis could be done to determine whether the atmospheric mechanism also produces\nsynchronizationwithFennoscandianicestreams. However,ifsuchevidenceislacking,thenitshould\nnot be claimed. Perhaps more accurate would be to say that it ”provides a simple explanation for\nthe observational evidence of synchronous discharge events from different ice streams within the\nLaurentide Ice Sheet”.\nReply: SincetheFennoscandianicesheetisheldfixedinoursimulations,wefollowtherecommendation\nfrom the reviewer and removed references to other ice sheets from the manuscript.\nReviewer Point P2.4 — The central idea in this study is the most similar to a recent study\nby Mann et al. (2021 in Paleoclimatology and Paleoceangoraphy) which is also able to explain the\nsynchronization of Heinrich and DO events, but instead using ice-ocean interactions to generate\nsynchronizationbetweenaflowlineglaciermodelandasimpleoceanmodel. Clearlythismanuscript\nproposes a different synchronization method and uses more realistic models, but that prior study\nprovides many potential points of comparison (including more comprehensive parameter sweeps).\nOnethingthatisnotcleariswhethertheMannetal. ocean-basedmechanismforsynchronization,if\n4\nincludedinamorecomplexmodellikethisone,couldalsoexplainalltheaspectsofthepaleoclimate\nrecord that the authors focus on in this study. The paper would benefit from a more extensive\ncitation and comparison to the previous study by Mann et al.\nReply: We have added a citation to Mann et al as well as a paragraph to the introduction that\nacknowledges that phase locking has also been achieved in simple coupled models of two autnomomous\noscillatorsthroughaweakcouplingifthecouplingisstrongenoughtoovercomethenaturalfrequenceies\nof the two oscillators.\nReviewer Point P2.5 — One surprising result is that surface temperature perturbations from\nDO events have any signficant impact on ice sheet sliding on time scales of a few thousand years.\nThe diffusion time scale of heat through thick ice is tens of thousands of years. This is why in\nthe original papers on the B-P theory for Heinrich events, MacAyeal discounts the role of surface\ntemperaturechangesindrivingHeinricheventsbecausetheywouldtaketoolongtodiffusethrough\nthe ice sheet. It would seem that a more robust explanation is needed here, particularly with some\nplots showing the influence of adding the atmospheric perturbation on the basal heat budget of ice\nstreams.\nReply: The reviewer is correct that the diffusion timescale is much longer than the few hundred years\nof warming during a DO cycle. Rather, the driving process here is that the surface warming results in\nan enhanced ice deformation because the ice deformation law is dependendent on ice temperature. This\nin turn raises the internal frictional heat generation and slightly warms the ice temperature at the ice\nbase. The effect is demonstrated in belows Hovmoeller plots.\nMinor\nReviewer Point P2.6 —Line87: Inwhatsensedoyoumean”statisticallysignificantdifferent”?\nWhat metrics characterizing the Heinrich event are you measuring and how is that difference\nmeasured?\nReply: We have added a short explanation what is being compared and what statistical test we use.\nWe also added a reference to the Methods section where this is also detailed.\nReviewer Point P2.7 — Line 90: Why does it seem that the SMB perturbation has more of an\neffect than the temperature perturbation?\nReply: We showed in an earlier publication [Schannwell et al., 2023] that Heinrich event periodicity\nis primarily affected by surface mass balance because of its climatic setting which is characterised by\ncold temperatures (well below -20°C) and by surface mass balance values between 0.1 - 0.15 m/yr.\nIn such a regime, the system reacts to surface mass balance perturbations much more because the\nrelative magnitude of the perturbation is much larger than for the surface temperature. We have added\na citation to our previous paper.\nReviewer Point P2.8 —Line93: Increasedsurfacetemperaturewillincreasesnowaccumulation\nat high elevation but also increase surface melting at low elevations? Is the latter not important at\nthe edges of the ice sheet?\n600simulationyears. ∆referstodifferencebetweenthemaximumtemperature,minimumSMBperturbation\nsimulation and the unforced reference simulation.\nReply: In reality this is correct, but our applied perturbations here are spatially homogeneous and act\nindependently. WehaveaddedmissinginformationabouttheSMBcomputationtotheMethodssection\n(see point below).\nReviewer Point P2.9 — Line 95: What is the mechanism explaining why ”The onset of the\nHeinrich event is, however, delayed into the early phase of the transition into the stadial of the\nDansgaard-Oeschger cycle”\nReply:Wehaverephrasedthissentencetohighlightthatitisnotsufficientforasinglelocationtoreach\npressure-melting-point and permit rapid basal sliding, but the warm-based region first has to reach a\ncertain size before a Heinrich event is triggered. The lag was referring to the time difference between\nreaching pressure-melting-point and for the region to grow large enough to trigger a Heinrich event.\nReviewer Point P2.10 — Line 99: Similarly, can you explain why it takes 3-5 D-O cycles before\nthe next Heinrich events occurs? Is it simply that H events are paced internally by ice streams?\nDo the Hudson Strait and MacKenzie ice streams have the same natural periodicity for discharge\nevents?\nReply: This is because the natural frequency dictated through the binge-purge is roughly of that or-\nder. Of course the natural frequency is modified by different background climate conditions and other\nperturbations. As shown in Schannwell et al. [2023], Hudson (7,200 years) and Mackenzie ice stream\n6\n(4,200 years) indeed do have very different natural surge frequencies. We have added a short mention\nof this to the main text.\nReviewer Point P2.11 — Line 140: Can you say with significance that an increase from 8% to\n17% or to 22% would not just happen by chance? How many more Heinrich events is this? Is it\nmore than one?\nReply: We have replaced “significant” with “substantial” to remove the statistical connotation of the\nsentence. We also added the information of how many more synchronous Heinrich events occurred\nduring the simulations. Provided that the signal is monotonously increasing with higher temperature\nperturbations and consistent across search windows from 100 years to 800 years, we are very confident\nthat the findings are robust.\nReviewer Point P2.12 — Line 158: What makes this an advantage? Why do you need this\nmechanism to apply to land-terminating glaciers? It wouldn’t help explain anything about the\nmarine sediment record which is only informed by marine-terminating glacier discharge.\nReply: We have clarified that we think that our mechanism provides an advantage in that it works\nequally well for permanently marine-terminating ice streams like the Hudson ice stream, but also for\nice streams that are land-terminating and only become marine-terminating during the actual Heinrich\nevent.\nReviewer Point P2.13 — L163: related to point #3 above, it would be helpful to explain\nhere how internal temperature and SMB are modeled in mPISM-VILMA since those processes are\ncentral to the mechanism in this study.\nReply: We have added the missing information to the Methods section.\nReviewer 3\nReviewer Point P3.1 — This paper presents a novel mechanism that is introduced to explain\nthe synchronous occurrence of Heinrich events and the stadial phase of Dansgaard-Oeschger cycles\nduring the last glacial period. Heinrich events, quasi-episodic iceberg discharge events from the\nLaurentide ice sheet into the North Atlantic, have been a significant feature of the northern hemi-\nspheric climate during the last ice age cycle and is a difficult process to study because of the large\nrange of temporal and and spatial scales that are involved and the small number of Heinrich events\nobserved. Using simulations with a coupled ice sheet-solid earth model, the researchers present a\nnovel mechanism that eliminates the need for a trigger mechanism during the stadial phase. In-\nstead, it is proposed that the interstadial phase of the Dansgaard-Oeschger cycle, characterised by\natmospheric warming, preconditions the ice sheet for a Heinrich event, causing it to occur during\nthe subsequent stadial phase.\nThis paper re-visits the binge-purge hypothesis and presents a fresh re-interpretation on the\nmechanisms behind Heinrich events. While I personally resonate with the subsurface ocean warming\n7\nmarineicesheetinstabilityhypothesiswhichoccursduringstadials, itisimportanttoprovideother\npossible theories that may help explain a dominant feature of ice age climate. I feel the paper is\nwell written and the study rests on some solid previous work by the group. I would like to see some\nof my concerns addressed which I list below.\nReply: We thank the reviewer for the positive assessment of our manuscript.\nReviewer Point P3.2 — General Structure of the paper:\nI like Figure 2 and I think this nice graphic should be the highlight of the paper. I would suggest\nmoving this to Figure 1. It should not be too hard to rework the text on pages 4 and 5. Start with\nPg 5 L92. I have more specific comments on the figure below.\nReply: Yes, we agree with that and have swapped the order of Fig. 1 and Fig. 2 as well as adjusted\nthe text to reflect the adjusted Figure order.\nMinor comments\nReviewer Point P3.3 — Figure 1: Why is there a colour wheel in the “Unforced reference\nsimulation” in figure 1a. Is this not forced with constant MIS3 forcing and so should be blank?\nOr have I missed something?\nReply: The reviewer is correct. We removed the colour wheel.\nReviewer Point P3.4 — The current Figure 1 could be enhanced with another schematic\ndiagram beside the top panel “Unforced reference simulation” (to the left, illustrating how the ice\nvolume builds as a curve spirals and then collapses with a Heinrich event. Also it was not initially\nclear what all the thin black lines were indicating (the ensemble of runs) and how many revolutions\nthey made. Do all the thin lines start at 0 time at some small ice volume? It will take the reader\nsome time to see what is happening so a simpler schematic of one simulation with text (possibly)\nhighlighting HEs and DO stadial and Interstadial periods.\nReply: We have added an annotated explainer subplot to the revised Figure.\nReviewer Point P3.5 — The matrix of plots in figure 1 has ∆SMB vs ∆Temp, shouldn’t this\nbe A vs A (the amplitudes given from the equations (3 and 4) in the Methods section. This\nsmb temp\ncould be changed and pointed out in the figure caption.\nReply: Since we think that ∆SMB vs ∆Temp are better understandable for the reader, we decided to\nswap symbols in equations 3 and 4 instead.\nReviewer Point P3.6 — What does the white significance level indicate (< 90th percentile) ,\nmaybe state this.\nReply: We have added this to the Figure caption.\nReviewer Point P3.7 — Figure 2: It would be useful to highlight that the insets to figure 2 are\nfrom equation 1 in the Methods section and the colours on the curves represent the derivative of\nthe temperature and mass balance changes.\nReviewer Point P3.8 — To most people I think they would associate the point in Figure 2c\nas early interstadial. It is so close to the D-O warming event that I don’t think one can possibly\nargue that this is a stadial condition.\nReply: We have changed the term to “late interstadial” in the revised version.\nReviewer Point P3.9 — Line 51-54: Is it possible to infer evidence of ice shelves from the\nDeVernal et al 2000 paper on LGM. Is there a strong statement in this paper about no ice shelves\nexisting along the Hudson Strait Ice Stream before H1?\nReply: We rephrased that sentence to state that “there is no firm evidence of the existance of such a\npersistent ice shelf”\nReviewer Point P3.10 — L67-68: On explaining the observed characteristics of the paleo\nrecord: The atmospheric forcing functon (Equation 1 Methods does not have the same quick drop\noff at the end of the interstadial as seen in the ice core records (e.g. Rassmussen et al 2014). In\nparticular if you look at GI-13 the interstadial period remains elevated for some time (800 years)\nand then quickly drops off. The results here possibly work here with H5, but if you look at GI-9\nthe interstadial is short (< 500 years). The interstadials previous to the other Heinrich events do\nnot fit well into the context of this paper and may have other origins other than the Hudson’s Bay\nIce Stream. It would be difficult to say anything using observations to corroborate the findings of\nthis study.\nReply: We are aware that our forcing function is idealised. To capture the effect of temporally varying\nDOcycles,weperformedadditionalsimulationsinwhichtheDOcyclelengthisrandomlyvariedbetween\n1,000 and 2,000 years (in 100 year intervals).\nReviewer Point P3.11 — L80: This would be a good place to highlight the amplitude from the\nmethods (A and A ) and refer to the methods.\nsmb temp\nReply: Yes, we added this.\nReviewer Point P3.12 — If you assume that there are Heinrich events at 30 , 39 and 47.5 kaBP\nyou would argue that there is a series of decreasing interstadial lengths previous to the H-Events\nas observed from the ice core data (a sort of damped D-O oscillation). I would like to see what\nkind of effect would result from a damped harmonic oscillator (e.g if you added an exp(-t/tau) in\nfront of the arctan term in Methods equation 1. From the results presented , I assume that the\nsynchronisation would decrease as the amplitude goes down. So I have a hard time seeing how this\ncould be corroborated with observations. I didn’t see any strong attempt to argue how the model\nresults compare with observations of D-O variability.\nReply: We agree that the synchronisation will decrease with a decreasing perturbation magnitude. As\nmentioned above, we did add simulations where we investigated the effect of temporally varying DO\ncycle lengths. However, from the modelling perspective there is no strong evidence that DO cyles are\n9\nthe result of a damped oscillator (see Klockmann et al. [2018], Vettoretti and Peltier [2018],\nVettoretti et al. [2022]).\nReviewer Point P3.13 —Ziemenetal2019(CP)includeoceanicforcing. Withouttheinfluence\nofluence\nof ocean effects it is really hard to argue that this is “the” mechanism that causes H-Events.\nReply: We have added to the Methods section that we also force the model with constant ocean\nconditions. Please note that we tested the sensitivity to the ocean forcing in our model setup in\nSchannwell et al. [2023] and found that Heinrich events from the Hudson ice stream were unaffected\nby changes in ocean temperature or sea level forcing.\nReviewer Point P3.14 — L96-97: Is this consistent with Ziemen et al 2019? How much is this\nin mSv (50)?\nReply: Yes this is consistent wit Ziemen et al. 2019. We added the peak discharge rate to the sentence.\nReviewer Point P3.15 — L102-104: These results suggest that the short D-O cycles would\nhave an H-Event close to the end of the interstadial, is this consistent with observations.\nReply: For some of the recorded Heinrich events such as HE6 and HE4, it is consistent, but for others\nit is not. The aim of our study is however more a proof of concept of the mechanism rather than aiming\nto reproduce the observational record exactly. Having said this, we believe that with more model tuning,\nit is certainly possible to match the observational record better than we do with our idealised forcings.\nReviewer Point P3.16 — L122-123: Figure 1 suggests that the SMB is the primary driver from\nfigure 1 if I read the significance correctly.\nReply: Yes, this is correct. See reply to Reviewer 2.\nReviewer Point P3.17 — L219-224: What is the range of the AMOC in the MIS3 simulation?\nIn the extended data figures, what is the change in radial ice volume in meters of RSL (it is only\ngiven in m3.\nReply: We have added these numbers to the revised manuscript.\nReferences\nM. Klockmann, U. Mikolajewicz, and J. Marotzke. Two amoc states in response to decreasing\ngreenhouse gas concentrations in the coupled climate model mpi-esm. Journal of Climate, 31\n(19):7969–7984, Oct. 2018. ISSN 1520-0442. doi: 10.1175/jcli-d-17-0859.1. URL http://dx.\ndoi.org/10.1175/JCLI-D-17-0859.1.\nC. Schannwell, U. Mikolajewicz, F. Ziemen, and M.-L. Kapsch. Sensitivity of heinrich-type ice-\nsheetsurgecharacteristicstoboundaryforcingperturbations. Climate of the Past,19(1):179–198,\nJan. 2023. ISSN 1814-9332. doi: 10.5194/cp-19-179-2023. URL http://dx.doi.org/10.5194/\ncp-19-179-2023.\nrelaxation oscillation. Journal of Climate, 31(9):3423–3449, Mar. 2018. ISSN 1520-0442. doi:\n10.1175/jcli-d-17-0559.1. URL http://dx.doi.org/10.1175/jcli-d-17-0559.1.\nG. Vettoretti, P. Ditlevsen, M. Jochum, and S. O. Rasmussen. Atmospheric co2 control of spon-\ntaneous millennial-scale ice age climate oscillations. Nature Geoscience, 15(4):300–306, Apr.\n2022. ISSN 1752-0908. doi: 10.1038/s41561-022-00920-7. URL http://dx.doi.org/10.1038/\ns41561-022-00920-7.\nThe manuscript on \"A mechanism for reconciling the synchronisation of Heinrich events and\nDansgaard-Oeschger cycles\" has been substantially revised in response to prior reviews. In general, I\nthink the manuscript has substantially improved and that the authors have been responsive to my\nprior suggestions and those of the other reviewers, particularly in clarity of the explanation, the claims\nbeing made and the relationship of this study to prior literature.\nIn this review, I can hone my critique to one critical aspect of the discussion that I do not feel has\nbeen resolved by the revision. In particular, this related to my previous review point (P2.5 in the\nresponse to reviewers document) on explaining the role of surface temperature perturbations in the\nphase locking mechanism described in the paper. Since this is the central aspect of this study which is\nnovel (the atmospheric mechanism that is), I think it is quite important that it be clear that surface\ntemperature perturbations do indeed play a role in the phase locking. I have two outstanding\nquestions:\n1. It is not clear to me that the surface temperature perturbations are strictly necessary (separate\nfrom the SMB perturbations) to achieve phase locking. The revised text says the following (my\nasterisk emphasis added): \"The warmer surface temperatures lead to increased ice deformation and\nenhanced frictional heat generation, while the higher surface mass balance flux increases ice stream\nthickening. *Together* these processes ensure that a critical ice thickness and basal temperature\nthreshold, required for the onset of the internal instability, is reached towards the end of the\ninterstadial phase\". Reading this and looking at figure 2, what I cannot tell is whether the surface\ntemperature perturbation is strictly necessary to produce phase locking. What would be helpful to\nshow is whether phase locking occurs in experiments with zero surface temperature perturbations and\nnon-zero SMB perturbation, or conversely in experiments with non-zero surface temperature\nperturbations and zero SMB perturbation.\n2. Part of the reason for my raising the issue above, is that I am still having a hard time\nunderstanding the mechanism proposed for surface temperature perturbation influence on basal\nsliding on time scales of less than 1000 years. The response says \"the driving process here is that the\nsurface warming results in an enhanced ice deformation because the ice deformation law is dependent\non ice temperature. This in turn raises the internal frictional heat generation and slightly warms the\nice temperature at the ice base.\" Where I have the disconnect is that the temperature effect and effect\nof ice deformation and internal heat generation should be confined to the upper 100 or so meters on\ntime scales of less than 1000 years (see MacAyeal 1993 in Paleoceanography for this argument\nexplicitly). On that same time scale, the change in temperature at the bed is likely to be minuscule\n(I'd estimate order 10^-11 C). Thus, even with feedbacks associated with changes in ice flow and\ninternal heat generation, those effects would still be confined to the upper part of the ice sheet. Given\nhow cold the ice is near the surface, the surface temperature perturbations aren't likely to produce\ninternal melt (which could conceivably drain to the base), and so I remain at a loss for how the\ntemperature signal (absent SMB forcing) is meant to propagate to the base of the ice sheet in less\nthan 1000 years. It is possible that in performing the vertical integration to solve the momentum\nbalance with the hybrid scheme described here, the change in temperature at the surface effectively\ngets included in the basal sliding momentum balance. However, if that is the case, then this effect is\nan artifact of the model formulation, not a physical mechanism for propagating surface temperature\nperturbations to the base.\nOther than these issues, I am happy with the manuscript and where it is, but given that the\noutstanding issues related to the central mechanism novel to this paper, it bears ensuring that it is (a)\nwell described, (b) shown convincingly to be occuring using careful experiments, and (c) due to real\nphysical processes and not model formulation.\nReviewer #3 (Remarks to the Author):\nThe Authors have addressed my comments and have done an excellent job in revising the manuscript.\nI have no further comments.\nResponse to the reviewers\nWe thank all referees for their thoughtful and thorough reviews of our paper. We appreciate you\ntaking the time to complete these reviews and welcome your helpful comments. In the following\nwe address their concerns point by point. Throughout this response to review document referee\ncomments are provided in regular, non- italic font text, our response comments are provided in red\nfont.\nTo accomodate the main comment of the reviewer, we have made the following main changes to\nthe manuscript:\n• We included results from six additional sensitivity simulations which show, if the forcing is\nlarge enough, one of the two considered atmospheric perturbations is sufficient to achieve\nphase locking.\nReviewer 1\nReviewer Point P1.1 — The manuscript on ”A mechanism for reconciling the synchronisation\nof Heinrich events and Dansgaard-Oeschger cycles” has been substantially revised in response to\nprior reviews. In general, I think the manuscript has substantially improved and that the authors\nhavebeenresponsivetomypriorsuggestionsandthoseoftheotherreviewers,particularlyinclarity\nof the explanation, the claims being made and the relationship of this study to prior literature.\nIn this review, I can hone my critique to one critical aspect of the discussion that I do not feel\nhas been resolved by the revision. In particular, this related to my previous review point (P2.5 in\nthe response to reviewers document) on explaining the role of surface temperature perturbations\nin the phase locking mechanism described in the paper. Since this is the central aspect of this\nstudy which is novel (the atmospheric mechanism that is), I think it is quite important that it be\nclear that surface temperature perturbations do indeed play a role in the phase locking. I have two\noutstanding questions:\nReply: We thank the reviewer for the positive assessment of our revised manuscript. We address the\nremaining concerns in detail below.\nReviewer Point P1.2 — It is not clear to me that the surface temperature perturbations are\nstrictly necessary (separate from the SMB perturbations) to achieve phase locking. The revised\ntext says the following (my asterisk emphasis added): ”The warmer surface temperatures lead to\nincreased ice deformation and enhanced frictional heat generation, while the higher surface mass\nbalance flux increases ice stream thickening. *Together* these processes ensure that a critical ice\nthickness and basal temperature threshold, required for the onset of the internal instability, is\nreached towards the end of the interstadial phase”. Reading this and looking at figure 2, what I\ncannot tell is whether the surface temperature perturbation is strictly necessary to produce phase\nlocking. What would be helpful to show is whether phase locking occurs in experiments with zero\nsurface temperature perturbations and non-zero SMB perturbation, or conversely in experiments\nwith non-zero surface temperature perturbations and zero SMB perturbation.\nperformed in total six new simulations. In three of these simulations, only the surface mass balance was\nperturbed and in the other three simulations only surface temperature was perturbed. The results show\nthat if the forcing is strong enough, one of the two forcings is sufficient to achieve phase locking. We\nhave added the corresponding figures to the extended figures section and added the following text to\nthe manuscript: “Additional sensitivity simulations reveal that for the stronger perturbation scenarios\n(∆Temp ≥ 10°C and ∆SMB ≥ 0.075 m/yr) only one of the two atmospheric forcings is required to\nachieve phase locking. Hereby, perturbations to only the surface mass balance result in Heinrich events\noccuring deeper into the stadial phase of the Dansgaard-Oeschger cycle.”\nReviewer Point P1.3 — Part of the reason for my raising the issue above, is that I am still\nhaving a hard time understanding the mechanism proposed for surface temperature perturbation\ninfluence on basal sliding on time scales of less than 1000 years. The response says ”the driving\nprocess here is that the surface warming results in an enhanced ice deformation because the ice\ndeformation law is dependent on ice temperature. This in turn raises the internal frictional heat\ngeneration and slightly warms the ice temperature at the ice base.” Where I have the disconnect\nis that the temperature effect and effect of ice deformation and internal heat generation should\nbe confined to the upper 100 or so meters on time scales of less than 1000 years (see MacAyeal\n1993 in Paleoceanography for this argument explicitly). On that same time scale, the change in\ntemperature at the bed is likely to be minuscule (I’d estimate order 10−11°C). Thus, even with\nfeedbacks associated with changes in ice flow and internal heat generation, those effects would still\nbeconfinedtotheupperpartoftheicesheet. Givenhowcoldtheiceisnearthesurface, thesurface\ntemperature perturbations aren’t likely to produce internal melt (which could conceivably drain to\nthe base), and so I remain at a loss for how the temperature signal (absent SMB forcing) is meant\nto propagate to the base of the ice sheet in less than 1000 years. It is possible that in performing\nthe vertical integration to solve the momentum balance with the hybrid scheme described here,\nthe change in temperature at the surface effectively gets included in the basal sliding momentum\nbalance. However, if that is the case, then this effect is an artifact of the model formulation, not a\nphysical mechanism for propagating surface temperature perturbations to the base\nReply: We agree with the reviewer that the effect of the surface temperature on the basal temperature\nisnotmassiveonthe1,500yearcycleconsideredhere. However,thereisanimportantdifferencebetween\nMacAyeal’sworkandtheworkwepresenthere. MacAyealassumesaspatiallyhomogenousicethickness\nof 3,000 meters for Hudson ice stream. While our ice stream shows a similar ice thickness near the\ndrainage divide, the Heinrich events are initiated close to the ice sheet margin and propagate upstream.\nIn this region, ice thickness is around 700 m prior to a Heinrich event in our model simulations. While\nthe reviewer is correct that the e-folding depth, assuming no vertical advection, is 75 m and 145 m\nfor 400 year and 1,500 year cycle respectively, the temperature warming at the base of the region\nstill amounts to 0.002°C and 0.1°C at 700 m depth for the largest of our temperature perturbation\n(14°C), respectively. Naturally, these numbers would only increase if vertical advection and frictional\nheat generation would also be considered. Overall, this certainly presents not a massive warming, but\nbecausethebaseoftheicesheetisalreadyveryclosetothepressure-melting-point,itisenoughtonudge\na few cells across the pressure-melting-point threshold, resulting in a different timing of the Heinrich\nevent. However, as shown in our additional sensitivity simulations, this effect only comes into play for\nsufficiently large enough surface temperature perturbations.\nit is, but given that the outstanding issues related to the central mechanism novel to this paper,\nit bears ensuring that it is (a) well described, (b) shown convincingly to be occuring using careful\nexperiments, and (c) due to real physical processes and not model formulation.\nReply: We hope that our above explanations and additions to the manuscript address the concerns of\nthe reviewer.\nReviewer 2\nReviewer Point P2.1 — TheAuthors haveaddressedmycommentsandhavedoneanexcellent\njob in revising the manuscript. I have no further comments.\nReply: We are happy that the reviewer is satisfied with our revisions to the initial submission.\nI am satisfied by the new simulations and comments made by the reviewers addressing my concerns.\nI have no further suggestions.\nReviewer 1\nReviewer Point P1.1 — I am satisfied by the new simulations and comments made by the\nreviewers addressing my concerns. I have no further suggestions.\nReply: We are happy that the reviewer is satisfied with our final revisions to the previous submission.\n1"}]

Decision letter