Hysteretic evolution of ice rises and ice rumples in response to variations in sea level

. Ice rises and ice rumples are locally grounded features found in coastal Antarctica and are surrounded by otherwise freely ﬂoating ice shelves. An ice rise has an independent ﬂow regime, whereas the ﬂow regime of an ice rumple conforms to that of the ice shelf and merely slows the ﬂow of ice. In both cases, local highs in the bathymetry are in contact with the ice shelf from below, thereby regulating the large-scale ice ﬂow, with implications for the upstream continental grounding line position. This buttressing effect, paired with the suitability of ice rises as a climate archive, necessitates a better understanding of the 5 transition between ice rise and ice rumple, their evolution in response to a change in sea level, and their dynamic interaction with the surrounding ice shelf. We investigate this behaviour using a three-dimensional full Stokes ice ﬂow model. The simulations span end-member basal friction scenarios of almost stagnant and fully sliding ice at the ice-bed interface. We analyse the coupling with the surrounding ice shelf by comparing the deviations between the non-local full Stokes surface velocities and the local shallow ice approximation (SIA). Deviations are generally high at the ice divides and small on the lee sides. On the 10 stoss side, where ice rise and ice shelf have opposing ﬂow directions, deviations can be signiﬁcant. Differences are negligible in the absence of basal sliding where the corresponding steady state ice rise is larger and develops a fully independent ﬂow regime that is well described by SIA. When sea level is increased and a transition from ice rise to ice rumple is approached, the divide migration is more abrupt the higher the basal friction. In each scenario, the transition occurs after the stoss side grounding line has moved over the bed high and is positioned on a retrograde slope. We identify a hysteretic response of ice 15 rises and ice rumples to changes in sea level, with grounded area being larger in a sea level increase scenario than in a sea level decrease scenario. This hysteresis not only shows irreversibility following an equal increase and subsequent decrease in sea level, but also has important implications for ice ﬂow model initialisation. The initial grounded area needs to be carefully considered, as this will determine the formation of either an ice rise or an ice rumple, thereby causing different buttressing effects.

Ice rises and ice rumples are locally grounded features surrounded by floating ice shelves and play a dual role in this context. Firstly, ice rises and ice rumples regulate the flow of ice towards the ocean through their buttressing effect (Favier and Pattyn, 2015; Barletta et al., 2018;Reese et al., 2018;Still et al., 2019;Still and Hulbe, 2021;Schannwell et al., 2020) and control the 30 migration of the continental grounding line . Secondly, past adjustments in local ice shelf flow dynamics can be inferred from ice rises by investigating, for example, isochronal structure and the development of features such as Raymond arches within ice rises (Raymond, 1983;Martín et al., 2006;Gillet-Chaulet and Hindmarsh, 2011;Hindmarsh et al., 2011;Drews et al., 2013Drews et al., , 2015Schannwell et al., 2019;Goel et al., 2020). The importance of ice rise formation and decay for continental ice sheet evolution (e.g., due to glacial isostatic uplift or changes in sea level) have been recognised in a number of 35 scenarios and show the key role that ice rises play in large-scale grounding line migration patterns over glacial cycle timescales (Bindschadler et al., 1990(Bindschadler et al., , 2005Barletta et al., 2018;Kingslake et al., 2018;Wearing and Kingslake, 2019).
In adopting terminology from Matsuoka et al. (2015), we identify ice rises as prominent grounded features with a distinct local radial flow regime, causing the flow of the surrounding ice shelves to divert either side of the feature. Ice rumples, however, generally form on less prominent bed highs and result in a predominantly unidirectional flow regime with the upstream ice 40 shelf flowing over the bed anomaly. Ice rises and ice rumples are found all around the perimeter of the Antarctic Ice Sheet, but the mechanisms governing the transition from one flow regime to the other have not yet been investigated and influences of the surrounding ice shelves on the local flow regimes have not yet been quantified. In order to explore these questions, we use the three-dimensional, full Stokes model Elmer/Ice to simulate idealised ice rises and ice rumples under various basal friction scenarios and sea level perturbations. 45 To quantify non-local effects from the surrounding ice shelves, we compare the full Stokes solutions with simpler iceflow approximations (e.g., the shallow ice approximation Hutter (1983); Greve and Blatter (2009)) which do not capture the stress transfer between ice shelf and ice rise. Furthermore, we investigate whether the locality of flow and basal sliding can be determined by examining the mismatch between the Vialov profile and the observed ice thickness. Using sea level perturbations, we explore whether ice rises and ice rumples respond hysteretically and whether multiple steady states exist for 50 a given set of boundary conditions by tracking the grounded area, upstream ice shelf velocity and dome position. Additionally, we investigate under which formation scenarios ice rumples reach a steady state and under which scenarios they are merely a transient feature during ice flow reorganisation.

Methods
Ice rises and ice rumples, and their surrounding ice shelves are investigated in steady state and transient scenarios using the 55 three-dimensional full Stokes numerical model Elmer/Ice .  (2) and (4) indicate cross-sections used for analysis, and (3) is the ice rise dome. Both (2) and (4) are cross-sections through the ice rise dome.
In (e), A marks the highest point of the bed anomaly. The x-direction corresponds with the along-flow direction, the y-direction corresponds with the across-flow direction and the z-direction corresponds with the elevation.

Governing equations
We adopt a coordinate system in which the predominant along-flow direction is aligned with the x-axis, the predominant across-flow direction is aligned with the y-axis and the z-direction marks elevation relative to sea level. The flow of ice is governed by the full Stokes equations, where τ is the deviatoric stress tensor, p is the pressure, ρ i is the ice density and g = gê z is the gravitational acceleration. We assume the ice to be incompressible, and so, the mass conservation equation reduces to The non-linear rheology of ice is modelled using Glen's flow law, which relates the deviatoric stress tensor, τ , to the strain rate 65 tensor,˙ , as where the effective viscosity, η, is Here, n is the Glen's flow law exponent, A is a rheological parameter primarily dependent on ice temperature. Since we assume 70 ice to be isothermal, A is set to a constant value in all simulations. The effective strain rate,˙ e , is calculated from the strain rate tensor,˙ , aṡ e = tr(˙ 2 )/2.

Boundary conditions
The upper surface, z = z s (x, y, t), and lower surface, z = z b (x, y, t), evolve subject to whereȧ s,b are the accumulation / melt rates at the ice-shelf surface (s) and ice shelf base (b), respectively.
The surface accumulation rate,ȧ s = 1.2 m a −1 , reflects the comparatively high rates observed at some ice rises in the Dronning Maud Land area of East Antarctica (Drews et al., 2013). The basal melt rate,ȧ b , beneath the ice shelf is defined as a function of ice thickness, H, based on the parameterisation used in Favier et al. (2016), where ice is grounded, and , where ice is floating,  where α is a tuning parameter and |x − x g | is the closest distance to the grounding line.
A constant flux of Q x=0 = 5.4 × 10 9 m 3 a −1 into the domain is prescribed at the upstream boundary, corresponding to an initial velocity of 300 ma −1 . Ice flows through a fixed calving front where ice is subject to sea pressure. At the lateral boundaries, no friction is applied and the flow velocity is subject to the Dirichlet boundary condition u · n = 0, where n is the 85 normal vector pointing outwards.
Ice in contact with the bed is subject to a non-linear Weertman-type friction law, where τ b is the basal shear stress, C is a constant friction coefficient, u b is the velocity tangential to the bed, and m is the friction law exponent. The position of the grounding line at each time step is determined by solving a contact problem (Durand 90 et al., 2009). The First Floating Elmer/Ice grounding line numerical implementation is used (Gagliardini et al., 2016).

Idealised model domain setup
The evolution of ice rises and ice rumples is simulated in a 60 × 60 km domain (Fig. 1). A bed anomaly is introduced and allows isle-type ice rises and ice rumples to form. The bed takes the form b(x, y) = b 0 + b a , where b 0 is a constant and b a is an anomaly with a flat top, defined as The centre of the bed anomaly is located at (x 0 , y 0 ) and σ controls the horizontal extent. The shape of the bed anomaly is broadly consistent with observations of ice rises, many of which have a plateau-shaped top that is near horizontal (e.g., Derwael Ice Rise (Drews et al., 2015) and ice rises in the Fimbul Ice Shelf (Goel et al., 2020)). All parameters used in the model are summarised in Table 1. 100 The ice thickness is initialised to 300 m throughout the domain, resulting in a geometry that is predominantly floating with a small grounded area at the bed anomaly. To ensure adequate resolution at the grounding line and ice divide, the mesh is refined in the area encompassing the ice rise with a resolution of 350 m (Fig. 2). For this, we use the meshing software Mmg (http://www.mmgtools.org/, version 5.3.10). This is in line with mesh resolution recommendations from other studies (Pattyn et al., 2013;Cornford et al., 2016), but is also the highest mesh resolution that is computationally feasible for the time scales 105 considered here. To account for a possible migration of the ice rise, the radial extent of the area of high resolution is 5 km from the initial grounding line. In the remainder of the domain, a mesh resolution of 2000 m is used. The mesh is vertically extruded resulting in 10 layers spaced equally apart and the horizontal mesh size is kept constant throughout the simulations.

Shallow ice approximation (SIA) comparison
The shallow ice approximation (Hutter, 1983;Greve and Blatter, 2009) describes the flow of ice in the absence of longitudinal 110 and transverse stress gradients and is composed of the deformational velocity (u d ) and basal sliding velocity (u b ) so that the total velocity is u = u d + u b . In SIA, only the vertical shear stress gradients are considered, so that the x-direction and y-direction deformational components of the velocity take the form where h(x, y) is the height of the ice surface relative to a reference horizontal plane. Here, ∇ h denotes the two-dimensional, horizontal gradient operator. We compare the velocity components only at the surface of the ice and also assume that temperature is constant, and so Eq. (10) reduces to The x-direction and y-direction basal sliding components take the form where H is the ice thickness and C b is the basal friction coefficient and relates to the full Stokes basal friction coefficient, C, as follows: 125 where N b = N b e z is the basal normal stress.

Comparison with the Vialov approximation
The Vialov profile (Vialov, 1958) is an analytical solution for an ice sheet profile in the case of a non-slip, flat bed and constant accumulation. The flow in an ice rise is predominantly radial from a point divide and so we use a radial flux condition 130 resulting in an ice geometry profile of the form and Both h 0 and L are calculated from a reference point on the surface of the full Stokes simulation output.
We compare only the lee profile of the ice rises to the Vialov profile as the bed is relatively flat in this area and we assume that small changes in bed topography are negligible. The profiles are compared for a central cross-section from the divide, extending in the along-flow direction into the ice shelf (Label (4) in Fig. 1c).

Design of transient simulations
To allow perturbation simulations to start from a steady state geometry, all simulations are run for 2000 years under constant forcing. Simulations are performed for three different basal friction coefficients C = 3.812 × 10 6 , C = 7.624 × 10 6 and C = 3.812 × 10 8 Pa m −1/3 s 1/3 , which we will refer to as low, intermediate and high friction scenarios, respectively. The In the low and intermediate scenarios, the ice rises transition to ice rumples at some stage during the sea level increase. In the high friction scenario, no such transition occurs after a sea level increase of 80 m. We therefore continue the increase of sea 155 level further at a constant rate of 0.02 m a −1 until the transition occurs. A reversal of the sea level perturbation is performed from a height of 155 m above the initial sea level.

Steady state analysis
After 2000 years of spin-up time, ice rises with a characteristic local flow regime develop in all three scenarios (Fig. 4). From  (Fig. 7).

Ice rise to ice rumple transitions triggered by sea level variation
To understand the response of ice rises and ice rumples with differing basal friction to sea level perturbation, we analyse the grounded area (Figs. 8a,b and 9a), dome migration (Fig. 10) is no abrupt change once a transition from ice rise to ice rumple has occurred (Fig. 8). This is in contrast to the high basal friction scenario, where there is an abrupt change in the upstream ice shelf velocity as a transition from ice rise to ice rumple is approached.   and decrease is performed for the low basal friction scenario starting from the steady states that emerged from the previous sea level perturbation cycle. The response of the grounded area, ice shelf velocity and thickness are calculated as described above (Fig. 8). The hysteresis cycles are now closed and independent of the initial conditions.
When sea level rise is halted in the high basal friction scenario prior to the unstable grounding line retreat (here at a sea level perturbation of 155 m), then the ice rise volume and grounded area also recover asymmetrically resulting in two viable states 225 for a given sea level displacement (Fig. 9).
We hand, during sea level decrease, the hydrostatic grounding lines have a more rapid response.

The influence of basal sliding on the geometry and transient behavior of ice rises
A number of previous studies have argued that basal sliding near ice rise divides is negligible because thermomechanically coupled models often predict ice significantly below freezing point at the ice-bed interface near the summits (Martín et al., 235 2009;Drews et al., 2015;Goel et al., 2020) and because many ice rises exhibit isochronal features called Raymond arches which do not form if basal sliding is dominant (Pettit et al., 2003;Martín et al., 2009). However, low and intermediate scenarios can be relevant in areas where Holocene marine sedimentation results in basal sliding in areas which have regrounded (Pollard et al., 2016). Moreover, differences between observed and simulated Raymond arches under a frozen bed assumption may indicate a delay or suppression of arch growth due to past or present basal sliding (Kingslake et al., 2016).  The greater velocity differences in the lower friction scenarios show that these ice rises are influenced more by the stresses in the surrounding ice shelf. Implications for the presence or absence of a fully local flow regime are twofold: (1) if basal sliding is negligible even in areas close to the grounding zone, then SIA is an appropriate modelling framework, for example, 260 when investigating the surface accumulation history using inverse methods (Callens et al., 2016), and (2) the basal boundary condition determines an ice rise's response to sea level perturbation.  The low and intermediate friction scenarios respond immediately to a rising sea level, with a retreat of the leeward grounding line accompanied by a stossward migration of the dome position. The ice rises progressively thin and eventually transition into ice rumples. There is no significant threshold behaviour between these two states and once the sea level increase is halted, the 265 system converges to a steady state ice rumple with the lee side grounding line located on the retrograde slope at the edge of the basal plateau. The summits are a few tens of meters above the ice shelf surface and the overall geometry is consistent with, for example, the ice rumple located in the Roi Baudouin Ice Shelf (Fig. 13). The minimum overriding velocities of 20 m a −1 are, however, significantly faster than the example observed at the Roi Baudouin Ice Shelf where the ice is effectively stagnant . The smooth transition of the low and intermediate friction ice rises into ice rumples reflects their strong 270 coupling to the surrounding ice shelf, highlighted previously. From a larger scale perspective there are no critical differences between ice rises and ice rumples in those scenarios other than the switch from a local to an overriding flow regime.   (Jezek, 2003).
Conversely, the high friction case only transitions to an ice rumple for sea level perturbations that are greater than what is expected in a glacial-interglacial cycle. In fact, there is no noticeable change in grounded area even for a sea level displacement of 50 m. This stability is in line with, for example, ice promontories at the Ekström Ice Shelf which show a comparatively weak 275 response to the thinning of their surrounding ice shelves (Schannwell et al., 2019). Grounding line retreat rates for higher sea level displacements then remain moderate as long as the leeward side remains grounded on a prograde slope. On a retrograde slope the ice rise becomes unstable and complete ungrounding occurs. We therefore conclude that after a transition from ice rise, there is a threshold basal friction beyond which a steady state ice rumple cannot form.
likely due to a greater grounded area (Fig. 12) and it is worth investigating whether inverse techniques used to predict the basal friction coefficient beneath pinning points produce results which remain valid regardless of the grounded area.
The required sea level perturbation for ungrounding clearly depends on the elevation below sea level of the bed protrusion, but the scenarios shown here with a maximum bed elevation of 80 m below sea level have many real world counterparts (e.g., Kupol Moskovskij, Kupol Coilkovskogo, Leningrad Ice Rise, Djupranen Ice Rise (Goel et al., 2020), Derwael Ice Rise 285 (Drews et al., 2015)). Our study suggests that features with a high basal friction have been and will remain stable local flow features even for comparatively large sea level perturbations. Moreover, it shows that ice rumples with comparatively low surface velocities as in the example provided in Fig. 13, are very unlikely a result of a deglaciated ice rise. An area that requires more investigation is the case of ice rises which do not conform to the plateau-shaped bed topography as prescribed here. The unstable retreat predicted in the high basal friction scenario suggests that ice rises located on retrograde slopes are critically 290 less stable for an equal amount of sea level displacement.

The hysteretic behaviour of ice rises over glacial cycles
In all basal friction scenarios, there are two differing ice rises for a given sea level (Figs. 8,9). These pairs differ in the basal melt rate applied (which is thickness dependent) and in the grounded area. Each pair corresponds to a low and a high buttressing case for which the averaged upstream ice velocity is used as a proxy (Figs. 8c,d,9b and Fig. A1b in the Appendix).

295
The differences in grounded area occupied for the individual pairs is asymmetric. In all cases, the pairs occupy virtually the same region on the obstacle's stoss side, but the extent of grounding on the plateau differs (Fig. 11). The thickness and slopes at the respective grounding lines are comparable, and therefore differences in basal melt (as parameterised in Eq. 7) are small, with differences of only 3.5, 3.0 and 2.4 % in the low, intermediate and high friction scenarios, respectively. The dynamic differences therefore stem mostly from the differing grounded areas that result in a differing form drag (Still et al., 2019) and 300 consequently a differing net resistance to the upstream ice shelf.
A relevant self-stabilising feedback occurs, whereby the ice rise height reduces and the divide migrates stossward during lee side grounding line retreat. This increases the lee side accumulation area in the vicinity of the dome, thereby increasing the ice flux on the lee side and slowing the grounding line retreat. In the same way, sea level decrease results in leeward divide migration, slowing the grounding line advance. Another mechanism that plays a role is the sensitivity of the grounding line to 305 bed shape, with hysteretic behaviour occurring due to the positioning of retrograde and prograde slope segments (Schoof, 2007;Pattyn et al., 2012;Haseloff and Sergienko, 2018;Sergienko and Wingham, 2022). In our study, we also observe grounding line migration patterns linked to the shape of the three-dimensional bed protrusion. Consequently, it matters how the ice rise and ice rumple geometries are initialised to begin with.
Although in our study, we have used a constant surface accumulation, we would expect orographic precipitation to enhance 310 the hysteretic behaviour. In future work it is worth investigating whether effects such as an increased melt rate also produce an hysteretic response in ice rises and ice rumples. Given that the grounded area and basal sliding determine the ice rise evolution, future simulations should include a more informative guess of the basal friction coefficients guided by, for example, seismic studies determining the bed properties (Smith et al., 2015). Inversion of the basal friction parameters from a thermomechanically coupled full Stokes model (Schannwell et al., 2019(Schannwell et al., , 2020 does provide some information in this regard, but also 315 contains lumped uncertainties, e.g., from ice rheology and uncertain boundary conditions. Another process not considered here is changes in the bed protrusion through glacial isostatic adjustment (Kingslake et al., 2018;Wearing and Kingslake, 2019).
The existence of multiple steady states means that the grounding lines of ice rises and ice rumples observed today are dependent on the local ice flow history during the last glacial cycle. Inversely, the dynamics and buttressing effect of ice rises and ice rumples are dependent on the initial geometry prescribed, which is typically unknown. The degree of buttressing is of 320 importance for determining the stability and evolution of the continental grounding line (Favier and Pattyn, 2015;Reese et al., 2018). The representation of ice shelves has been identified as a key cause of continental-scale model spread (Seroussi et al., 2019) and a precise representation of ice rises and ice rumples would reduce spin-up and projection uncertainties.
We have shown that the difference between the simulated grounding line and the hydrostatic equilibrium grounding line is small at each time step. This small error may, however, lead to an error propogation during transient simulation leading to 325 inaccurate grounding line migration if a hydrostatic equilibrium assumption is used.

Conclusions
We examined the effect of basal friction and sea level variation on the evolution of ice rises and ice rumples using idealised simulations including the surrounding ice shelves. In a high basal friction scenario, there is negligible mismatch when comparing simulated full Stokes velocities with SIA velocities, whereas in a low basal friction scenario the mismatch is larger due to 330 a greater influence of stresses from the surrounding ice shelf. The locality of the ice flow and the degree of basal sliding can be diagnosed by examining the (mis-)fit of a Vialov profile to the observed thickness profile. In response to an increasing sea level, a transition from ice rise to ice rumple occurs. Steady state ice rumples form in the low basal friction scenarios whereas the ice rumple in the high friction scenario is ephemeral and ungrounds rapidly. The higher friction ice rise, on the other hand, is largely unresponsive to sea level variations, requiring more than double the sea level rise to trigger the transition compared 335 to the lower friction scenarios. All basal friction scenarios show self-stabilising, hysteretic behaviour, with grounded area and upstream ice shelf buttressing dependent on the evolution history. As a consequence of this behaviour, we identify the need for careful consideration of the grounded area of an ice rise during model initialisation in order for the correct feature to form. Although in our study, we have concentrated only on the response of ice rises to sea level perturbation, further processes such as an increase in basal melt are 340 also likely to result in hysteretic and potentially irreversible behaviour in ice shelf buttressing upstream of ice rises.
Code availability. The code used to run the simulations and the post-processing code can be found at https://doi.org/10.5281/zenodo.