From Heinrich Events to cyclic ice streaming : the grow-and-surge instability in the Parallel Ice Sheet Model

Here we report on a cyclic, physical ice-discharge instability in the Parallel Ice Sheet Model, simulating the flow of a three-dimensional, inherently buttressed ice-sheet-shelf system which periodically surges on a millennial timescale. The thermo-mechanically coupled model on 1 km horizontal resolution includes an enthalpy-based formulation of the thermodynamics, a non-linear stress-balanced based sliding law and a very simple sub-glacial hydrology. The simulated unforced surging is characterized by rapid ice streaming through a bed trough, resulting in abrupt discharge of ice across the grounding 5 line which is eventually calved into the ocean. We identify and visualize the central feedbacks that dominate the sub-sequent phases of ice build-up, surge and stabilization which emerge from the interaction between ice dynamics, thermodynamics and the sub-glacial till layer. A reduction in the surface mass balance or basal roughness yields a damping of the feedback loop which suggests that thinner ice sheets may be less susceptible to surging. The presented mechanisms underlying our simulations of self-maintained, periodic ice growth and destabilization may play a role in large-scale ice-sheet surging, such as the 10 surging of the Laurentide Ice Sheet, which is associated with Heinrich Events, and ice-stream shut-down and reactivation, such as observed in the Siple Coast region of West Antarctica.

Editor Decision: Reconsider after major revisions (08 Jun 2017) by Andreas Vieli Comments to the Author: Dear Authors, After 2 in general positive reviews of the first version of the manuscript and with some rather substantial suggestions for changes/revisions of one of the reviewers, the authors undertook major revisions and added additional experiments and modelling investigations and addressed the more technical points as well.As the revisions were rather substantial the paper has been sent to 2 further re-reviews, of which one of the reviewers had already reviewed the first version of the manuscript and was happy with the revised version of the manuscript.The second new reviewer had however some additional in my view valid comments that should be addressed before publication and which I outline below (the detailed review is given further below).The main points to be addressed are in brief: 1. to clarify the difference in the model and results of this study in relation to an earlier study by VanPelt and Oerlemans (2012) who used almost the same model.This should include -a clarification in the model difference -some better discussion/integration of the own and maybe new findings against/with Pelt and Oerlemans ( 2012) -perhaps a clearer discussion of what the implications of 'ocean-terminating' really changes /means (compared to the land-based VanPelt Oerlemans (2012)

study).
We agree with the Editor that a better highlighting of the differences between the models and the applied experimental setups is of great help in order to clarify the difference between the two studies.Accordingly, we added a paragraph in the introduction referring to vanPelt and Oerlemans (2012) who also used PISM, briefly mentioning the main differences (p. 2, ll. 31-35).We detail these differences in the model section (p. 4, ll. 2-5, ll. 8-9) and the experimental setup section (p. 5, ll. 1-9).The main differences that we carve out lie in the friction law (linear vs. non-linear friction coefficient) and in the model setup (land-terminating glacier vs. marine ice-sheet-shelf system), besides other differences (spatially uniform vs. non-uniform surface mass balance and surface temperature, horizontal diffusion of basal melt water).In this respect, we also revised the discussion/conclusion section, adding a more detailed comparison between our results and those from vanPelt and Oerlemans (2012), discussing the influence of the differences in the model/setup mentioned above on the results in three new paragraphs (p. 10, l. 16 -p.11, l. 6).This also includes the implications from the qualitatively different experimental setups, i.e., land-based glacier with unconfined tongue-shaped terminus vs. ocean-terminating ice-sheet-shelf system with the implication of non-buttressed vs. buttressed ice flow (p. 10, ll. 23-28).
2. Related, this study should be embedded slightly better into the surge-modelling literature in the introduction and discussion and therefore VanPelt and Oerlemans should likely be mentioned in the introduction, as their study is very close (with the difference of being land-based).
We agree that this study was mentioned quite late in our previous version of the manuscript and that it makes a lot of sense to refer to it much earlier in the text.Following the suggestion by the Editor, now we also mention it in the introduction in the list of similar surge-related studies (p. 2, l. 11), as well as in a separate paragraph (p. 2, ll.31-34, see our comment above).Furthermore, as suggested by the Reviewer, we included a couple of references to surging mountain glaciers, extending our introduction which only considered cyclic ice streaming and large-scale ice-sheet surging before (p. 1, ll. 15-18).
3. link to HE-events (overselling): the used geometry seems pretty generic and perhaps together with the results not too close to the real HE-case (the modelled cycle-periods also do not match and in reality there was likely a slight retrograde slope into Hudson bay which would change things quite a bit regarding groundingline motion and feedbacks).Thus although a motivation by the HEevents seems ok, maybe the title could be rethought and be less strongly linked to HE-events (HEevent 'like', 'towards' HE-event cycles….???).At the end it is the choice of the authors of how to label the manuscript but i think one should be careful to also deliver to the raised expectations.Note, in the context of HE-event mechanisms and for completeness, in the introduction perhaps the recent literature of Alvarez-Solas (2013) and Bassis et al (2017) (see detailed references further below) are may be relevant (although not directly related to the mechanism investigated here, so not crucial).
We see the Editor's point and at the same time are glad that he understands our idea of using Heinrich Events as a motivation.Following his suggestion we changed the title of the manuscript to "From cyclic ice streaming towards Heinrich-like events: the grow-and-surge instability in the Parallel Ice Sheet Model" in order to have a weaker link to Heinrich Events while at the same time keeping the motivation via Heinrich Events.We are glad to include the Editor's references to further studies that investigate externally-forced periodic surging and also added another study that proposed that Heinrich Events might not be triggered by purely ice intrinsic mechanisms (p. 2, ll.3-5).

the choice and implications of a rather low reference velocity (1m/y) should be addressed/ or clarified in relation to the variable flow law exponent 'q' (and related issues, see reviewer comment)
Following the advice of the Reviewer, we included a paragraph discussing the choice of the parameter value, its meaning for the q-phi parameter study and also its implications regarding the differences in results between the two PISM studies (p. 10, l. 34 -p. 11, l. 5).

I agree with the reviewer that the 'averaged grounded' visualization is not that ideal, and makes
it also difficult to compare to other studies (e.g. VanPelt and Oerlemans 2012).At least for some of the experiments maybe a average or max.centreline value could also be shown (as additional dashed lines or in appendix).
We also agree that it makes a lot of sense to give an alternative to the averaged values.For this purpose we generated two figures that show the centerline average as well as two point values on this centerline (as also suggested by the Reviewer), analogous to Figs. 3 and A1 (now A2), respectively.As suggested by the Editor, we added these Figures to the Appendix (Figs.A1 and A3)

For details see reviewer comment below Editor Andreas Vieli
Reviewer comment to be addressed: In their manuscript From Heinrich events to cyclic ice streaming, the grow-and-surge instability in the Parallel Ice Sheet model, Johannes Feldmann and Anders Levermann nicely describe their findings from surge-type experiments with the Parallel Ice Sheet Model.In this respect, this manuscript is very similar to the Van Pelt and Oerlemans 2012 publication.There is nothing fundamentally wrong with this.Testing previously obtained results is an important process in science, and it is interesting to see that and under which parameters the current version of the Parallel Ice Sheet Model shows surge-type behavior.
We would like to thank the Reviewer for taking the responsibility of reviewing our manuscript, the helpful advice and constructive criticism.We think that the Reviewer's valuable comments and suggestions, which we address below, were most useful in improving our manuscript.
However, it then needs to be made very clear, where previous results were (not) confrmed and where new results are obtained.In contrast to the Author's claim, Van Pelt and Oerlemans also used a model version that included modeled pore water as enthalpy (page 348, second paragraph).The friction laws in the two publications look very similar to me.Please point out where the 'simpler friction law` was improved.The visible main differences between the publications are the investigation of a water-terminating glacier instead of a land-terminating one (this has no obvious effect on the surge cycle as drawn in Figure 2), the parameters that were varied, and slight changes in the treatment of basal water.
We agree with the Reviewer that the differences between the models and also the obtained results were not made very clear in the previous version of our manuscript.Thus in the revised version we added a paragraph to the introduction, briefly mentioning the main differences between the models and experiments (p. 2, ll. 31-35).Extending the methods section, now we detail the differences in the friction law (linearly vs. non-linearly dependent on basal till water [p. 4, ll. 2-5] and the treatment of the basal till water [p. 4, ll. 8-9]).We also describe the differences in the experimental setup (boundary conditions and in particular the bed topography) in more detail and its consequences for the resulting qualitatively different types of modeled ice bodies (p. 5, ll.1-9).Revising the discussion section, we added three paragraphs that discuss the main differences between 1) the models, 2) the experimental setups and 3) the velocity scaling parameter value and their effect on the results (p. 10, l. 16 -p.11, l. 6).We removed our statement on the enthalpy scheme.
The manuscript left me a bit concerned regarding the appreciation of the authors for previous work, and a tendency towards over-selling their own results.This begins in the title where Heinrich events are mentioned right at the start, while the manuscript deals with cyclic surging in a pretty generic ocean-terminating glacier, that is not obviously set up to resemble the situation in duthe Hudson Bay / Hudson Strait area, the source region of Heinrich events.I would suggest cutting the title to The grow-and-surge instability in the Parallel Ice Sheet model or something similar.
We understand the the Reviewer's concern that the title might promise too much since our study is indeed not designed to reproduce Heinrich Events.Our idea of including Heinrich Events in the title was to use it as a motivation and also to reach a broader community.In order to weaken the link to Heinrich Events we followed the suggestion by the Editor and changed the title to: "From cyclic ice streaming towards Heinrich-like events: the grow-and-surge instability in the Parallel Ice Sheet Model".
While continuing through the manuscript, I would have enjoyed seeing more references to previous publications on glacier surging, and the similarities and differences in the underlying theory and results.This has already substantially improved since the first version of the manuscript, but I think it can be improved further.One example is the introduction of Figure 2 in Section 3.1, where a comparison of the model presented here to the binge-purge cycles of MacAyeal 1993(and subsequent publications, e.g, Roberts et al, 2016) could make it easier for the reader to appreciate the additions presented in this manuscript.There also is a substantial body of research on growthand-surge cycles in the context of mountain glaciers that could be referenced.
Following the Reviewer's suggestion we extended the first paragraph in Section 3.1, setting our study into context of previous work (p. 5, ll. 18-28).In this respect, we would also like to point to the introduction where we give a brief overview on previous surge-related work and also highlight the additions made by our study (p. 2, ll. 3-30).Regarding the literature on surging mountain glaciers we were glad to add several references to the introduction (p. 1, ll. 15-18) to also mention this type of surging (besides large-scale ice-sheet surging and cyclic ice streaming).
The newly introduced comparison of different values for the basal sliding exponent q suffers from an ill-chosen reference velocity uc = 1m/s.This is the velocity for which the basal friction _b is independent of the exponent q.For comparing different flow law exponents, this value has to be in the range of typical sliding velocities.Using numbers obtained from the surge phase in Figure 3, the basal shears stress for ice sliding at ub = 1000 m/a, and a _c value of 300 kPa varies over more than six orders of magnitude for the different values of q.It is highest for the plastic flow law with _b = 300 kPa (q=0) and decreases via _b = 9500Pa (q=1/3 ) to _b = 9:5Pa (q=1).For lower velocities relevant for surge initiation, this contrast is even more extreme.Van Pelt and Oerlemans chose uc = 100m/a, a much more realistic value.This might explain some of the differences in the results of the two studies.
We would like to thank the Reviewer for pointing to this important issue.Consequently, we included a paragraph in the discussion/conclusion section, discussing the implications of our choice of the velocity scaling parameter on the q-phi parameter study and on the comparability of the results mentioned by the Reviewer (p. 10,l. 5).Figs. 3,[6][7][8][9]A1, considering that at least two thirds of the grounded domain are not affected by the surge at all.I would probably prefer one or two point measurements in the center line, or a line-average of this line.

I am not convinced by averaging the values for
We agree with the Reviewer that values along the centerline give deeper insight into the surge dynamics.Following the Reviewer's suggestion now we also show curves for the centerlineaverage as well as two point measurements on the centerline (Figs.A1 and A3).

Model
We use the open-source Parallel Ice Sheet Model (PISM; Bueler and Brown, 2009;Winkelmann et al., 2011;PISM authors, 2017), version stable07 (https://github.com/pism/pism/).The thermo-mechanically coupled model applies a superposition of the shallow-ice approximation (SIA; Morland, 1987) and the shallow-shelf approximation (SSA; Hutter, 1983) of the Stokes stress balance (Greve and Blatter, 2009).In particular, the SSA allows for stress transmission across the grounding line and thus accounts for the buttressing effect of laterally confined ice shelves on the upstream grounded regions (Gudmundsson et al., 2012;Fürst et al., 2016).The ice rheology is determined by Glen's flow law (Cuffey and Paterson, 2010).An energyconserving enthalpy formulation of the thermodynamics in particular allows for an advanced calculation of the basal melt rate for polythermal ice (Aschwanden et al., 2012).The model applies a linear interpolation of the freely evolving grounding line and accordingly interpolated basal friction, and uses one-sided differences in the driving stress close to the grounding line (Feldmann et al., 2014).

Experimental setup
The three-dimensional setup is designed to model a marine ice sheet, which drains through a bed trough, feeding a bayshaped ice shelf which calves into the ocean.The idealized bed topography (Fig. 1) is a superposition of two components: , is an inclined plane, sloping down towards the ocean (Fig. 1b).The component in y direction, b y (y), has channel-shaped form (Fig. 1c) and is a widened version of the one used in the MISMIP+ experiments (Asay-Davis et al., 2016, here with adjusted parameters for domain width and channel side-wall width, see Table 1).The superposition of both components yields a bed trough which is symmetric in both x and y directions (symmetry axes x = 0 and y = 0).While the main ice flow is in x direction (from the interior through the bed trough towards the ocean) there is also a flow component in y direction, i.e., from the channel's lateral ridges down into the trough.Resulting convergent flow and associated horizontal shearing enable the emergence of ice-shelf buttressing, which leads to a groundingline position further downstream than in the absence of an ice shelf.Ice is cutoff from the ice shelf and thus calved into the ocean beyond a fixed position (Fig. 1a).There exist more sophisticated methods to represent calving in a numerical model (e.g., Nick et al., 2010;Levermann et al., 2012;Pollard et al., 2015) leading to more realistic calving behavior and calving-front geometry.Our simple approach of prescribing a calving front (fixed in space and time) far enough from the region of the confined ice shelf makes sure that the calving front does not interact with grounding-line migration in the course of the surge cycle.:::: Due :: to ::: the :::::::: symmetry :: of ::: the ::::: setup ::: we :::: only ::::::: consider ::: the ::::::::: right-hand ::: half ::: of ::: the :::::: domain ::::::::: throughout ::: our :::::::: analysis.(in contrast to Van Pelt and Oerlemans, 2012, where surface mass balance and temperature are parameterized by the bedrock elevation, pre There is no melting beneath the ice shelf.Glacial isostatic adjustment is not accounted for in the experiments.The simulations are initiated with a block of ice from which an ice-sheet-shelf system evolves while ice flow, basal mechanics and till-stored water content adjust.This spinup lasts a few 1000 yr and thus we focus on the time after this phase.Due to the symmetry of the setup we only consider the right-hand half of the domain throughout our analysis.
The model is run using finite differences and a regular grid of 1 km horizontal resolution.An initial examination of the flow field reveals that the SIA velocities are small compared to the SSA velocities in our simulation.Despite this fact, considering the SIA in the simulations in particular allows for the representation of a three-dimensional temperature field.
On the slow time scale, the ice sheet tends to grow toward an equilibrium thickness, which is determined by the balance between snowfall and ice flux (velocity).If this equilibrium ice thickness is too large to be sustained by the basal conditions, this build-up (negative gray feedback loop in Fig. 2) is interrupted by an abrupt surge event with a rapid, self-enforcing speedup of the ice flow (positive feedback loop in red).The associated large-scale ice discharge into the ocean eventually leads to a stabilization of the shrunken ice-sheet-shelf system (negative feedback loop in blue), which again tends to restore a balance thickness before a new surge event kicks in.
At the beginning of the modeled surge cycle the negative feedback loop of slowing-down ice growth is dominant: the basal till water content drops close to zero and basal friction is high, allowing gradual thickening of the ice sheet (Figs. 3,4 and A1).The thickening causes an increase in basal melt water production due to the lowering of the pressure melting point at the ice base.The increasing water content in the basal till attenuates further increase in basal friction (which still increases due to the effect of ice thickening, i.e., growing overburden pressure P o , see Eq. 2), leading to an increase in ice discharge and thus reducing further thickening.In the absence of any other mechanisms, the ice sheet would hence reach a steady state as ice thickening would approach zero, eventually.
However, the continuous accumulation of water in the sub-glacial till during the slow build-up initiates a surge event before the equilibrium thickness is reached.The self-enforcing feedback of rapid ice speed-up becomes dominant: lowered friction at the well lubricated ice-sheet base leads to an acceleration of ice flow through the bed trough (Fig. 3 ::: and :: ??).In turn, this causes an increase in strain and frictional heating due to enhanced shearing inside the ice sheet and sliding of the ice over the bed, respectively (Fig. A1 :: and ::: ??).The resulting additional melt water production further lubricates the ice base, leading to even more speed-up (termed "hydraulic runaway" by Fowler and Johnson, 1995).Inside the bed trough, the previously relatively stagnant ice flow has entered a state of rapid ice streaming (velocities at several km yr −1 , Figs. 1a, b and 4d).The ice streaming is additionally fostered by the effect of strain heating at the side margins of the trough (Fig. 4): faster flow causes stronger shearing of the ice, resulting in more heat production which in turn softens the ice, allowing for more shearing and thus flow acceleration (so called "creep instability", Clarke et al., 1977, see positive feedback loop in our Fig.5).
The ice streaming inside the bed trough leads to enhanced downstream advection from the ice sheet's thick interior into the ice shelf, manifesting a pronounced peak in ice discharge and iceberg calving, respectively (Figs. 3d and e).The associated damping feedback between ice velocity (ice flux) and thickness (blue stabilization loop in Fig. 2) counteracts the self-enforcing feedback between ice velocity and till water (red surge loop).On the long term, this discharge-related thinning of the ice sheet leads to the end of the surge as melt water production decreases, basal friction increases and ice flow decelerates.When the ice sheet has become too thin to maintain insulation of its base from the cold atmosphere then the basal melt rate drops.The associated decrease in till water now is amplified through the same mechanism that was responsible for the till water increase during the surge phase (red loop) and the ice stream shuts down.At some point basal refreezing sets in, consuming further water from the till layer (Fig. A1).As the water content in the drained till drops close to zero and thus bed friction quickly increases, the ice sheet can build up again.The period duration of a whole surge cycle is of about 1800 yr, from which the slow build-up phase takes more than ::::: about 80 %.
We would like to note that the domain-averaged yield stresses resulting from our simulations (order of ∼ 100 kPa, see Fig. 3c) are relatively large compared to values from in situ and laboratory experiments (order of ∼ 1 kPa to ∼ 10 kPa, see Table 7.5 in Cuffey and Paterson, 2010).In our experiments the highest values occur during the build-up phase which is when the water content in the till is very close to zero (s ≈ 0) for which Eq. ( 2) yields τ c (s = 0) ∼ 10 5 kPa.Though in the model τ c is limited by the overburden pressure, the maximum possible value in our experiments is still on the order of ∼ 10 3 kPa (assuming an ice thickness of 1000 m; for a visualization see Fig. 1 in Bueler and van Pelt, 2015).Such large values occur predominantly in the regions outside of the bed channel and in the thick interior of the ice sheet where the basal till layer is continuously dry and ice flow is stagnant, biasing the domain-average towards high values.In contrast, inside the lubricated bed channel the simulated yield stresses are much lower, especially during the phases of ice streaming (on the order of ∼ 10 kPa), lying within the observatory range.This is in accordance with the fact that the observational values were inferred from till samples stemming from regions of relatively fast ice(-stream) flow (e.g., Truffer et al., 2000;Tulaczyk et al., 2000a;Kamb, 2001).

Surge damping
Varying the bed strength in our simulations, we find that surging is maintained in a cyclic manner (oscillatory equilibrium) only if the bedrock roughness allows the evolution of an ice sheet of medium thickness.For rather slippery basal conditions, realized by low values of the till friction angle, φ ≤ 8 • , and thus rather thin ice sheets, surging occurs initially but then is damped such that on the long term the ice sheet reaches a non-oscillating stable equilibrium state (Fig. 6).Decreasing the value of φ within this regime leads to faster damping and a shorter cycle duration.For sufficiently lubricated (thin) ice sheets no surging takes place at all.In contrast to the case of maintained cyclic surging, the ice flow enters a state of continuous streaming at velocities of several 100 m yr −1 with stable till water thickness (Figs.6b and c).
Vice versa, increasing the friction angle towards large values yields rougher beds, promoting the evolution of thicker ice sheets which surge at larger magnitude and lower frequency (Fig. 7).For sufficiently strong beds with φ ≥ 60 • (and thus comparatively thick ice sheets) initial surging is damped, similarly to the case of relatively low values of φ discussed above.
Consequently, surging is maintained only in a regime of medium bed strength (medium values of φ) that promote ice sheets of medium thickness.Damped surging occurs on both ends of this regime (above an upper and below a lower critical threshold of φ), i.e., for relatively strong and weak beds.
We investigate changes in the ice-flow characteristics close to the lower regime boundary in response to a small modification of the basal roughness.For this purpose we perturb the above equilibrium ice sheets of oscillatory (φ = 10 • ) and non-oscillatory (φ = 8 • ) type by decreasing/increasing the value of φ.Our results show that when the friction is lowered from φ = 10 • to values of φ ≤ 8 • then the originally surging ice sheet undergoes damping and enters a stable steady state, eventually (Fig. 8a).Hence, the flow characteristics of the perturbed ice sheet are more or less the same as in the spin up experiments when using the same values of the friction angle (compare Figs. 6a and 8).In contrast, perturbing the system in the other direction, i.e. increasing the friction angle from φ = 8 • (Fig. 8b) then maintained surging does only occur for values of φ ≥ 20 • (compared to φ = 10 • for the case of ice-sheet spinup).For lower values of φ the ice flow starts to oscillate initially but then goes back into a state of stable flow at the same velocity as before, whereas now ice-sheet thickness and till water content are larger (both increasing for increasing φ).Thus, the ice sheet in stable equilibrium requires a comparatively large perturbation of the basal conditions in order to turn into a state of maintained surging.In contrast, a small perturbation is sufficient to bring the continuously oscillating ice sheet into a stable steady state.
The finding from above, that thin ice bodies are less likely to surge than ice sheets of medium thickness, is supported by additional experiments with reduced surface accumulation a.According to these simulations, lower accumulation results in thinner ice sheets, a weaker surge amplitude and a longer surge-cycle duration (Fig. 9).The longer cycle duration can be explained by the fact that less snowfall causes the ice sheet to take longer to grow thick enough to trigger a surge event.At the same time the formation of basal till water during build-up takes longer and the kick-off of the surge event requires a smaller amount of till water (and thus a thinner ice body).Below a threshold of a fifth of the default value (a = 0.075 myr −1 ) a rather thin steady-state ice sheet forms and surging is not existent anymore.

Role of basal sliding law
The above results show cyclic or damped surging for a confined set of parameter values (default surface accumulation a and till friction angle φ given in Table 1 are only slightly varied).These simulations use a particular non-linear sliding law, determined by a basal sliding exponent of q = 1 3 (Eq.1).In general, this exponent can range from q = 0 (purely plastic sliding law) to q = 1 (linear sliding law).To explore the influence of the basal sliding law on the ice flow behavior, we conduct further simulations, sampling q between 0 and 1 at an interval of 1 12 .For each applied parameter value of q the till friction angle φ (Eq.2) is varied between 5 • and 85 • , spanning a wide range from relatively slippery to very rough bed conditions, respectively.The resulting q − φ parameter space is explored in terms of surge-cycle duration and ice-sheet volume (Fig. 10).Due to the large number of simulations the experiments are carried out on a grid of 5 km horizontal resolution (in contrast to 1 km used in the default simulations).
The results show that (damped) surging occurs in a range from q = 1 12 to q = 3 4 (circles and triangles in Fig. 10).Within this regime larger values of q correspond to higher friction angles φ, i.e., going towards a more linear friction law requires a rougher bed in order to observe (damped) surging.Maintained surging occurs in a smaller range, i.e., from q = 1 6 to q = 5 12 .This regime is embedded such that the transition from the oscillatory state into the stable regime in most cases leads through the damped regime.Generally, decreasing q or increasing φ yield a longer period duration of the surge cycle while the mean grounded ice mass increases.This can be explained by considering the relevant acting stresses: a larger value of φ leads to a larger magnitude of the basal yield stress (Eq.2) and thus a stronger basal shear stress (Eq.1).In the shallow-shelf approximation the driving stress (due to surface slope of the ice sheet) is balanced by a combination of the membrane stresses (responsible for ice-flow acceleration) and the basal shear stresses (see Eq. 17 and the following paragraph in Bueler and Brown, 2009).An increase in basal shear thus slows down ice speed-up, promoting a longer period of ice-sheet build-up and larger ice-sheet thickness.Decreasing the exponent q leads to the same results because also here the magnitude of the basal shear stress increases.This becomes evident from Eq. ( 1) where the fraction (|u b |/u 0 ) q ::::::::: (|u b |/u 0 ) q increases with decreasing q since |u b |/u 0 < 1 ::::::::: The duration of the surge cycle ranges from about 1800 yr to 2700 yr for maintained surging (mean ≈ 2200 yr) and from 800 yr to 5700 yr for damped surging (mean ≈ 2600 yr).
For the particular cases of purely plastic (q = 0) or linear sliding (q = 1) no surging occurs in our simulations, which is independent of φ.In the vicinity of q = 0 most of the experiments produce a stable and rather thick ice sheet (squares in Fig. 10), whereas around q = 1 the ice sheets become very thin (a few 10 m of thickness).In some cases these very thin ice sheets do not have any grounded ice inside the bed channel and thus lack comparability (marked by an "x" in Fig. 10).
In general, very small values of φ cause continuous streaming of a rather thin ice sheet on a slippery bed, whereas large φ values lead to rough basal conditions allowing the evolution of a comparatively thick steady-state ice sheet (Fig. 7).Thus, only those ice sheets which are not too large or small show surging behavior.This confirms and generalizes our specific results from Sec. 3.2 that there is a thickness regime in which surging occurs whereas too thin or too thick ice sheets reach a stable equilibrium.

Discussion and conclusions
We model the cyclic surging of a three-dimensional, inherently buttressed, marine ice-sheet-shelf system (Fig. 1).Periodically alternating ice growth and surge are unforced and emerge from interactions between the dynamics of ice flow (evolution of velocity, internal and basal stresses, ice thickness), its thermodynamics (heat conduction, strain and basal frictional heating, melt-water production) and the subglacier (melt-water storage and drainage).
We identify three consecutive phases throughout the surge cycle (ice build-up, surge and stabilization), each characterized by a dominating feedback mechanism which we visualize in a feedback-loop scheme (Fig. 2).These feedbacks of slowing-down ice thickening, rapid ice speed-up and discharge, and decelerating ice thinning (Figs. 3 and A1) can explain central processes that likely prevailed during repeated large-scale surging of the Laurentide Ice Sheet and the associated Heinrich Events of global-scale impact.During the surge phase mainly the process of hydraulic runaway (positive feedback between basal melt water production and flow acceleration; Fowler and Johnson, 1995) is in effect.It is complemented by creep instability (positive feedback between strain heating and ice deformation; Clarke et al., 1977), which additionally promotes rapid ice streaming (Figs. 4 and 5).The modeled cyclic alternation of ice streaming and stagnation provides a simple example of ice-stream shutdown and re-activation, a phenomenon which is characteristic for the dynamics of some of the Siple Coast outlets in West Antarctica.
Our results suggest that medium-sized ice sheets are more susceptible to cyclic surging than rather thin or thick ones.We find a transition from surge to non-surge behavior (surge damping) of the ice flow when decreasing/increasing the thickness of the surging ice body in our simulations, realized by applying lower/larger basal roughness or surface mass balance (Figs. 6   and 9) or by a variation of the friction exponent in the sliding law (Fig. 10).This is consistent with the existence of a critical minimum ice thickness found by Schubert and Yuen (1982).According to their results, exceeding this thickness threshold enables the occurrence of creep instability, potentially leading to rapid surging.Furthermore, our results reveal that an ice sheet in stable equilibrium requires a comparatively large perturbation of the basal conditions in order to turn into a state of maintained surging, whereas a small perturbation is sufficient to bring the continuously oscillating ice sheet into a stable steady state (Fig. 8).
Compared to the observed interval of about 7, 000 yr at which Heinrich Events re-occured during the last glacial period (Hemming, 2004), our modeled surge-cycle period of ∼ 2, 000 yr is much shorter.This is not surprising, given that our idealized model setup on a synthetic bed geometry is not designed and the parameters are not tuned to represent conditions that prevailed for the prehistoric Laurentide Ice Sheet.Thus, we refer to studies designed to model this ice sheet when it comes to the proper representation of the characteristic surge frequency of Heinrich Events (e.g., Marshall and Clarke, 1997;Calov et al., 2002;Papa et al., 2006;Roberts et al., 2016).Our model results are closer to results from conceptual studies which also use an idealized geometry (e.g., Bougamont et al., 2011;Van Pelt and Oerlemans, 2012;Robel et al., 2016).These studies all yield a surge-cycle duration of ∼ 1, 000 − 2, 000 yr, despite considerable differences in degree of physical approximations, parameterizations and complexity in setup geometry.However, all of them use a Weertman-type, stress-balanced based sliding law (Eq.1) and are based on the same (though individually modified) sub-glacial model (Tulaczyk et al., 2000a, b), suggesting that both have a strong imprint on the surge-cycle duration.
Several other parameters in our model likely have an effect on the occurrence of surging and its dynamics (e.g., the overburden-pressure fraction δ in Eq. 2, the till drainage rate C d in Eq. 3, as well as surface temperature, geothermal heat flux and bed slope).However, further investigation of the parameter-dependency of the surging behavior (e.g., as done for surface temperature and geothermal heat flux in Robel et al., 2014) is beyond the scope of this study.In fact, it aims at reporting on the realization of cyclic surging/ice-streaming of an ice-sheet-shelf system in the Parallel Ice Sheet Model based on suitable model components and justified set of parameters.

Figure 2 .Figure 3 .Figure 5 .Figure 6 .
Figure2.Schematic visualizing the three main feedback mechanisms, each of them dominating one of the three sub-sequent phases of slow ice build up (gray), abrupt surging (red) and stabilization (blue), forming a full surge cycle.The sign next to an arrow pointing from variable A to B indicates whether a small increase in variable A leads to an increase (+) or decrease (-) in variable B. According to this convention one can deduce from counting the negative links of a full loop whether this loop describes an amplifying (positive) or stabilizing (negative)feedback.An even number of negative links indicates a positive feedback loop (large +) whereas an odd number of negative links indicates a negative feedback loop (large -).

Figure 10 FigureFigure A2 .Figure
Figure 10.(a) Surge-cycle duration and (b) mean grounded ice mass for the q − φ parameter space.Each colored rectangle represents a simulation characterized by either oscillatory surging (white circles), damped surging (triangles) or stable equilibrium (squares).White rectangles with an "x" denote parameter combinations for which no grounded ice forms inside the bed trough and are thus not considered in the analysis.Since the simulations of stable ice flow do not exhibit periodicity by definition the associated rectangles in panel (a) are colored in gray.The default simulation with parameters of q = 1 3 and φ = 10 • is highlighted by a purple circle (see Figs. 3 and A1, gray curve in Figs. 4, 6, 7, 8a and 9).

Table 1 .
Physical constants and model parameters