Sensitivity of Heinrich-type ice-sheet surge characteristics to boundary forcing perturbations

Clemens Schannwell, Uwe Mikolajewicz, Florian Ziemen, and Marie-Luise Kapsch

 General comments

This is an interesting paper that seeks to advance understanding of the North Atlantic’s Heinrich Events (HEs), which are well documented sediment layers strongly suggesting large-scale periodic surges from the former Laurentide ice sheet. The paper is a significant advance in this field and will be important to communities working on paleoclimate as well as the ice-sheet instability. It investigates the impact of ice sheet external forcing and boundary conditions on surge characteristics by way of a series of sensitivity experiments using a standalone ice-sheet model. It contains a description and analysis is the basic surge mechanism from a control experiment, as well as analysis of paired (high-low) experiments in which the magnitude and frequency of external forcing is varied. The text is, on the whole, clearly written and structured, and the figures well presented and informative.

My comments are mostly around the interpretation and discussion of results. I have six major concerns that require either additional/revised experiments or appropriate rewriting.

1. The authors should consider refocussing the Introduction to better reflect the remainder of their paper, which focusses on internally generated surges as opposed to mechanisms that rely on interactions with other components of the Earth system. This refocussing might be informed by the following suggestion.
2. One of the more interesting aspects of the work is this discovery of two types of ice-sheet surging. One (Mackenzie) is closely allied to the original binge-purge concept of MacAyeal, while the other (Hudson) has a very different mechanism linked to the upstream propagation of sliding. Much more could be made of the differences between these two surge mechanisms, such as why they occur in their specific geographical locations (perhaps linked to channels in the subglacial topography) and the contrasting impact of changes in forcing on their dynamics.
3. There are issues with the design of the anomaly experiments. I believe the anomalies are applied instantaneously at the start of each experiment so that there will be a long-term trend through the experiment as the modelled ice sheet gradually responds to the often large changes in forcing. These experiments may not therefore measure the full impact of the anomaly; this would require a separate spin-up for each anomaly experiment. This needs to be made very much clearer in the text and the implications for the results discussed. Alternatively, the anomaly experiments could be rerun with appropriate spin up.
4. The surging behaviour observed is not regular and exhibits variations in duration of active/quiescent period, as well as magnitude, in a single experiment. This means the results need to be treated very carefully. The situation is compounded by the small sample size of surges within each experiment (typically ~10), which makes statistical analysis difficult. It would be interesting to investigate whether a series of experiments with a particular type of anomaly (for instance varying from Smb- through Ctrl to Smb+) would support the conclusions drawn here or whether surge irregularity and sampling issues would create a more complicated picture. Finally, the ranges chosen for each pair of anomaly experiments seem fairly arbitrary so that it is problematic to compare their impacts on surging (apples and oranges issue). Linking these ranges to the variability exhibited in the original MIS3 experiment may be one way of tackling this issue.

5. The discussion concerning phase locking and synchronicity (both with other elements of the Earth system and between individual ice streams) is a little naïve. To my mind, there are genuine questions concerning whether phase locking is possible given the small number of surges/oscillations possible within the lifetime of the Laurentide during a single glacial period. None of the experiments presented unequivocally demonstrates phase locking (which would require for instance the phasing of the two ice streams’ surges to evolve toward synchroneity through time) and much of the discussion is therefore speculative.
6. The discussion on the interaction of the two surging ice streams by migration of the ice divide is interesting but lacks detail on the mechanisms involved. In particular, it is unclear how this type of interaction could occur given that the time taken for ice to flow from the divide to an ice stream is well in excess of the period of the modelled surges.

Specific comments (with line number)

18. The Introduction suggests that there are two main schools of thought on the cause of HEs: internal-generated ice sheet oscillations and oscillations that rely on the interaction of the ice sheet with other components of the Earth system. The paper focuses of the former and is based around a model that exhibits internal oscillations of this type. The Introduction, however, mainly focusses on interactions. This seems odd and the Introduction would benefit from refocussing toward a fuller discussion of the literature around internal oscillations.
19. “implicated” need to make the link between the sliding perturbation of preceding sentence and ocean warming. Why do changes in sliding necessarily suggest that ocean forcing is important?
20. “only” see comment about strengthening review of literature on internal oscillation mechanism. For instance, Roberts et al. 2016 (Climate of the Past) perform a similar sensitivity analysis. There are also several papers by Marshal and Clarke that are also relevant.
21. My understanding of the Holland Jenkins (1999) three-equation model is that they use the temperature of the ocean sublayer immediately under the ice shelf as opposed to the ambient ocean temperature. How is this sublayer temperature found?
22. The spatial distribution of basal heating over a wider area is crucial to the propagation of the surges so that this detail is likely to be very important to the results reported elsewhere in the paper. The authors should explain and justify this choice and, if possible, describe its impact on their results.
23. “better” this Is very subjective. In what way is the simulation better?
24. Need to say how long the experiments were? I guess 60 kyr based on the figures but would be good to say so.
25. Hudson ice stream is marine terminating. The movement of its grounding lines may therefore be important. Please include information on how the grounding line was modelled.
26. The initial state was taken from a transient simulation (MIS3) so presumably is not in equilibrium with the climate forcing at that time (36 ka in MIS3). Is this an issue and is there any evidence of model drift in additional to the internal oscillations?
27. It would be good to summarise the forcing strategy here as well. My understanding is that forcing for the period 36 ka to 23 ka (36 ka minus 13 kyr) is taken from the MIS3 simulation as a starting point and then adapted by taking the temporal mean or applying anomalies etc. This 13 kyr (as two 6.5 kyr cycles) is then repeated for the duration of the experiment (60 kyr).
28. “favourable” this is again subjective. Please define what is meant by this and explain why it is important.

121-129. I found this description confusing. Figure 2 helps but I think there is room for a clearer explanation. Perhaps it would be helpful to explain why this seemingly complicated process is necessary.

139. 3D fields. Not clear to me what is meant here. Surface temperature and SMB are inherently 2D fields. Is the 3^rd dimension time or vertical dependence (perhaps used in an interpolation downscaling scheme)?
140. Ocean forcing. This relates back to the use of the three-equation model.
141. Is the same initial condition used in the anomaly experiments as the control? This implies that the forcing (after anomaly applied) will be out of equilibrium with the initial state (even more so than the issue on line 119). Is this a problem? Presumably there will be an overall trend in the anomaly experiments as the ice sheet responds. How does this relate to the internal oscillations? I don’t see much evidence that their character changes through each individual experiment.

This also raises the issue of what actually is being done in these anomaly experiments. If the control initial condition is being used in the anomaly experiments, the anomalies may not have enough time during the experiment to affect much change (to say ice thickness or internal temperature field). The analysed impact of a particular anomaly might therefore be its instantaneous impact as opposed to its larger, longer-term impact.

Table 1. Unclear what “experiment type” refers to.

147. I think it would be worth being even more explicit about the forcing here: there is no temporal variation in forcing so that any surging behaviour observed must be due entirely to internal mechanisms.

Figure 3. Why is the unit m²/yr? Is this the average flux through the gate per m width? I would have expected m³/yr.

163. This is one of the most interesting parts of the paper – the difference types of surge in the Mackenzie and Hudson ice streams. The authors should consider comparing their surge mechanism with the paper that first described this style of surging (Payne et al 1995 JGR) which contains a more process-based analysis of the mechanism than Calov et al. (2002) including the roles of ice temperature and stress concentrations in the initiation and termination of a surge.

This style of surge is related to but distinct from the original binge-purge theory of MacAyeal (1993). In the latter, the bed of the ice sheet interior warms leading to a pocket of fast flow that propagates towards the margin. In the former, stress concentrations at interface between fast and slow flow leads to local warning and the upstream migration of this interface. It is fascinating that both types of surge are found in different parts of the ice sheet, and the text could usefully be refocussed around discussing this distinction.

169. Grows is imprecise – is it ice thickening or later extension that is important?
170. It is not clear why delta stress (driving stress minus basal shear stress) is a meaningful metric. Heat generation by friction within the ice sheet and at its bed are both functions of stresses and strain rates. I don’t see how delta stress helps in this analysis.

The links between this description (163-177) and Figure 4 could be tighter (only three of the panels are referenced).

Figure 4. X axis title and labelling missing.

178. Not clear what “ice-stream front” refers to. Marine margin? Use of the term front also crops up later as well and should be tightened (they are a number of ‘fronts’ of different types in the model).
179. It don’t think it is clear from Figure 6 that propagation is downstream, in fact panel e suggests warming at the margin before warming in the interior. This might be a case of selecting a better time for the pre-surge time slice. Also giving times for each phase would help, i.e. what is the time difference between panels e and h?

Here again it would be worth emphasizing that this mechanism is closer to the original binge-purge concept. Also worth mentioning that Mackenzie is land terminating while Hudson is marine terminating.

184. I don’t think there is enough evidence to state that the warmer air temperatures are the “likely” cause. Can you exclude other components of the heat budget such as variations in geothermal heat flux and/or reduced vertical advection?

Figure 5. Why does surface gradient have units of m^-1. Should this not be nondimensional? Need to clarify relation between distance on x axis and Figure 1 (i.e., is it measured from upstream or downstream end). Similarly, clarify what “ice stream front” means. Also need to explain and justify the choice of surface gradient as a metric. Is this because of its role in controlling gravitational driving stress/frictional heat production? If it is simply as a map of the extent of the surge then why not use ratio of surface to basal velocity (of just the surface velocity itself).

196. The ranking of different boundary forcing anomalies seems arbitrary. The individual +/- experiments are interesting but it is hard to compare them to make statements about which type of forcing is more important. Some attempt to justify the anomaly range is made for each forcing variable, however I do not think this is precise enough to allow comparison between them. For some of the perturbations it would make sense to make the +/- range a fraction of the variability seen in the MIS3 simulation (e.g., SMB, surface temperature, sea level, ocean temperature), although this would not work for geothermal heat flux.

196 Another issue with this comparison is that the surges are at least in part chaotic, i.e. the magnitude and active and passive duration of individual surges is not the same even in one experiment. The length of each experiment is also short compared to a surge so that each contains only ~10 surges. It is therefore difficult to compare details. Performing a number of experiments in each +/- forcing range would shed light on this issue.

Figure 4 and 6. Need an indication of the time for each phase so that they can be related to Figure 3.

214. Is the geothermal heat flux a constant? Need to clarify.
215. You need to be very careful with statements such as “most sensitive to changes in SMB” this really assumes that all of the perturbations have the same ‘strength’, which I do not think you are able to do. The only way of doing this would be to use a common metric, i.e., a fraction base on their variability (see 196).

Figure 8. I estimate that each experiment contains around 10 surges so that illustrating their variability using standard deviation is overkill. Why not show all individual surges in an experiment as a cloud of separate points (with same colour)? Also need to indicate how period is calculated/defined. Difference in time between maximum and minimum ice volume?

265. Worth stressing here that the ice streams have two very different surge mechanisms hence forcing anomalies likely to affect them in different ways.
266. Comparison between the present day and LGM is not a good a measure of the validity of these anomalies (the Laurentide ice sheet does not exist at the present day). Better to compare against observed or modelled changes for times when the ice sheet existed at roughly the size that you are simulating (i.e. variability during the glacial period).
267. This may also link to the saddle collapse mechanism of Gregoire et al. (2012).

Figure 9. These plots strongly suggest that there is a trend in the model’s response in the anomaly experiments (i.e., that the anomaly is imposed instantaneous at the start of the experiment). As indicated elsewhere, the ice sheet will require time (thousands of years) to achieve its full response to all of the anomalies by either evolving its internal temperature field or its geometry. The use of instantaneous anomalies does not therefore measure the full response of the ice sheet to the anomaly. For instance, increased surface air temperature by itself will have little immediate effect but might have a strong effect after this warming has had time to propagate down through the ice to impact basal thermal regime. The same is true of all of the anomalies with the possible exception of geothermal heat flux (because it is applied at the bed), although even here there is likely to be a difference between the instantaneous and long-term response.

Figure 9. Would be good to include experiments in figure caption (also for Figures 11 and 12).

292. Give logic for why this was only done for SMB and GIA.
293. My understanding of phase locking is that many oscillations are required for it to develop. I do not see how it can happen given the small number of ice sheet surges that are possible during a glacial period. The oscillations reported here have period ~10 kyr so at most 10-20 oscillations are possible during the lifetime of the Weichselian Laurentide ice sheet. Phase locking would require a mechanism that can continue through interglacials in the absence of the ice sheet.
294. Comments about the number of oscillations required for phase locking are also relevant here (even more so since synchronisation of surges can definitely only happen while the ice sheet exists).
295. Potential interaction between the Hudson and Mackenzie by ice divide migration is interesting, however the present discussion is a bit simplistic and ignores the time required for a change in upstream catchment area to impact the ice flow within the ice stream. The ice streams are typically 1000 km apart (Figure 1) suggesting a travel distance of 500 km from divide to ice stream. Given ice velocities of ~10 m/yr upstream of the ice streams (Figures 4 and 6), a travel time of ~50 kyr is required. This is similar to the period of the surge limiting the potential for interaction.
296. This paragraph is speculative. Finding two experiments (of 17) that have the same phase is not evidence for the possibility of phase locking. Stronger evidence could come from a series of experiments based on Geo- and St+ that show the development of locking (i.e., oscillations that are initially out of phase but become locked). This also relates to conclusion on line 382.
297. More generally, plotting the ice speed (or flux) for Mackenzie and Hudson flux gates against one another would be a good way of displaying potential synchronicity. This type of phase diagram would show synchronicity as a closed loop.
298. “positive climate perturbation” is hard to understand – needs to be more specific.

Technical corrections

28. glacier/s/
29. because of /the/
30. i/n/crease
31. missing full stop.
32. sensitive/ly/
33. sensitive/ly/

This authors present experiments performed with an ice-sheet model coupled with a GIA model over the North American domain pertaining to representative MIS-3 boundary conditions. As built, the ice sheet exhibits very large magnitude, quasi-periodic oscillations primarily in the Hudson Strait ice stream and the Mackenzie ice stream that are supposed to correspond to Heinrich Events. Along with an in-depth analysis of the dynamics of these ice surges, different purturbations in the boundary conditions are tested to assess their influence on the characteristics of the surges. The topic is interesting and open in the literature, and the contribution here with a state of the art model setup is valuable. The paper is also well written and flows nicely. Nonetheless, I believe framing of the context of these experiments needs to be refined and the credibility of the model setup producing the surges themselves needs more support. 

Framing. The Abstract and Introduction is very nicely written. However, the description of Heinrich events here based on strong assumptions resulting from the authors' own modeling work. Not all characteristics are known to be true, and thus should not be presented as definitively representative of reality. This is particular of the sentence on L3 in the Abstract and the paragraph on L43-56 in the Introduction. It is well known that the AMOC collapse related to Dansgaard-Oeschger Events occurs prior to Heinrich Events in sediment records (e.g., Barker et al., 2015, http://dx.doi.org/10.1038/nature14330). In other words, glacial stadials and the dominant climatic impacts recorded globally, are mainly related to ocean circulation changes. It may be that freshwater discharge during Heinrich Events can act to extend the duration of the stadial, but this is far from clear in the paleo records. Thus, discussing all of the "impacts" of Heinrich Events is not really relevant, at best, and rather misleading at worst. I would simply recommend removing this sentence and this paragraph from the manuscript, as it is not necessary for motivating the study of Heinrich Events, which are large scale and enigmatic in any case. 

Another example is the phrase "during which large amounts of ice are discharged" on L2 in the Abstract. We have no constraints from reconstructions on the total ice mass that was discharged during HEs. Some modeling studies do not show such large fluctuations in volume (e.g. Alvarez-Solas et al., 2013), while the model the authors use does. So, again, I would suggest more precision when discussing what we know versus what comes out of this work. This will facilitate putting the results in context of the open questions on this topic, and will increase the value of the work. 

Experimental setup. The model components used appear to be state of the art and include many/ most processes one would be interested in for studying this problem. However, the results exhibit ice velocity values that sound simply incredible. In the surge phase, it appears that maximum values of up to 40,000 m/yr are reached, with inland values over a distance of 2000 km of > 1,000 m/myr (see Fig. 4, bottom-center panel and Fig. 5b). This seems to me extremely implausible, and at the very least should be acknowledged as such in the text, or explicitly justified and explained.

Furthermore, I think such a result is predominantly dicated by the choices regarding basal friction, which are only described in a very cursory way here. In particular, I think the reader should see more information about how τ_c is calculated in the basal friction law, as well as the relation used to calculate effective pressure, which I assume is an input to τ_c. Rather than "basal sliding is discouraged", I think actual values should be given and any spatial maps used should be included here. The point is, not only do I think this is valuable information for understanding the nature of the surges presented here, but all information should be provided to allow others to reproduce the results. 

I am also surprised by the choice of a constant value of GHF of 42 mW/m2 everywhere. In the motivation for the range of the Geo+/- experiments, the authors cite a spatial GHF map (Lucazeau, 2019) of which several exist. Why not apply such a map as the default distribution and scale from that? 

Specifically with regards to Geo-, GHF values of < ~40 mW/m2, and particularly of 0 mW/m2 do not seem likely to exist in reality. The motivation of "Geo-" is therefore difficult to understand. Since a very low value of GHF is imposed by default, I would expect only a "Geo+" purturbation to be needed.

With regards to setting the imposed value of GHF to -1000 mW/m2, I am very surprised that this did not cause numerical problems for the model. It seems that the authors were able to achieve their goal of deactivating streaming, but out of curiosity, what does the vertical temperature profile of ice look like in these regions? To keep conditions within more realistic bounds, the authors may consider simply modifying the spatially-defined friction coefficient in these regions for these tests. 

Finally, I think overall some additional figures would help. The readers should see some maps of the whole ice sheet being simulated in Ctrl (surface elevation, velocity, basal temperatures, τ_b) and possibly in the different phases (Quiescent, Pre-surge, Surge). For example, a supplementary movie, showing the surge behavior in the control simulation would be very valuable as well. 

In summary, with a change in framing combined with additional information, this paper will be a valuable contribution to the literature on the still open question regarding the nature of Heinrich Events. Note that some additional specific comments are listed below. 

Specific comments

L1 and L18: "among the most dominant" <= What makes them most dominant? Consider rephrasing simply to "prominent"

L35: More recent theories propose => A more recent theory proposes

L38: tidewater glacier => tidewater glaciers

Paragraph L43-56: As mentioned above, this paragraph gives a misleading characterization of HEs. Climate impacts are largely due to AMOC shutdown. In addition, MWP1a is a very specific event during the deglaciation (L53-56). It cannot be considered as representative of Heinrich Events in general. Thus citing it here is not really appropriate. 

L74: cylce => cycle 

L106: 36 ka => 36 ka ago [or before present?]

L109: bathymetrie => bathymetry

L164: Consider also citing Feldmann and Levermann (2017), as they describe this mechanism in detail with experiments using PISM. 

L184-185: What about the impact greater thickness, insulating and warming the base upstream? This seems like it would also contribute significantly.

L214-215: Map of Lucazeau (2019) only goes down to values of about 40 mW/m2. 

L252-253: It is hard to compare magnitude sensitivity, given that the units are totally different. Consider removing this sentence, as it is not necessary, or at least reformulating. Perhaps you can quantify the sensitivity in terms of Δt surge time as a function of percent change in boundary variable. Maybe this would allow a more appropriate comparison.

L258: sensitive => sensitively

L293: What kind of frequencies are imposed on GHF simulations? Are these not constant? Or rather this is GHF at the based of the ice sheet, which can evolve transiently, as opposed to GHF deep in the bedrock? Please clarify here more explicitly.

L298: sensitive => sensitively

L367-370: Again, this statement may be true, but it cannot be justified by the units used. Is a perturbation of +100 mm/yr of SMB equivalent to a surface temperature increase of +5 K? Find a way to rephrase, or rescale.