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 floating ice shelves. An ice rise has an independent flow regime, whereas the flow regime of an ice rumple conforms to that of the ice shelf and merely slows the flow 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 flow, 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 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 flow model with idealised ice rises and ice rumples. 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 stoss side, where ice rise and ice shelf have opposing flow directions, deviations can be significant. Differences are negligible in the absence of basal sliding where the corresponding steady-state ice rise is larger and develops a fully independent flow 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 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 shows not only irreversibility following an equal increase and subsequent decrease in sea level but also that the perturbation history is important when the ice rise or ice rumple geometry is not known. 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.

1 Introduction

Great progress in ice flow modelling has improved the physical representation of dynamical processes at the margins of the Antarctic Ice Sheet, but the transient evolution of the grounding line continues to be challenging, requiring high mesh resolution, small time steps, and advanced model physics (Schoof, 2007; Goldberg et al., 2009; Gudmundsson et al., 2012; Haseloff and Sergienko, 2018; Sergienko and Wingham, 2022). Moreover, the lack of past observational constraints and ice sheet model initialisation inconsistencies result in spin-up simulation geometries which differ from observations (Seroussi et al., 2019) and result in parameter choice uncertainty (Albrecht et al., 2020).

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 influence the migration of the continental grounding line (Favier et al., 2012). 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., 2013, 2015; Schannwell 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) has been recognised in a number of scenarios and shows the key role that ice rises play in large-scale grounding line migration patterns over glacial cycle timescales (Bindschadler et al., 1990, 2005; Barletta 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 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.

To quantify non-local effects from the surrounding ice shelves, we compare the full Stokes solutions with the shallow ice approximation (Hutter, 1983; Greve and Blatter, 2009) and the Vialov profile (Vialov, 1958), 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 full Stokes ice thickness and the Vialov profile, an idealised solution for the ice geometry. Using sea level perturbations, we explore whether ice rises and ice rumples respond hysteretically and whether multiple steady states exist for a given set of boundary conditions by tracking the grounded area, upstream ice shelf velocity, and dome position. Additionally, we investigate the formation scenarios under which ice rumples reach a steady state and the scenarios under which they are merely a transient feature during ice flow reorganisation.

2 Methods

Ice rises and ice rumples and their surrounding ice shelves are investigated in steady-state and transient scenarios using the three-dimensional full Stokes numerical model Elmer/Ice (Gagliardini et al., 2013).

2.1 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,

$$\begin{array}{}\text{(1)}& \mathbf{\nabla }\cdot (\mathit{\tau }-P\mathbf{I})+{\mathit{\rho }}_{\mathrm{i}}\mathit{g}=\mathrm{0},\end{array}$$

where τ is the deviatoric stress tensor, P is the pressure, ρ_i is the ice density, and $g=-g{\widehat{\mathit{e}}}_z$ is the gravitational acceleration. We assume the ice to be incompressible, and so, the mass conservation equation reduces to

$$\begin{array}{}\text{(2)}& \mathbf{\nabla }\cdot \mathit{u}=\mathrm{0}.\end{array}$$

The non-linear rheology of ice is modelled using Glen's flow law, which relates the deviatoric stress tensor, τ, to the strain rate tensor, $\dot{\mathit{ϵ}}$, as

$$\begin{array}{}\text{(3)}& \mathit{\tau }=\mathrm{2}\mathit{\eta }\dot{\mathit{ϵ}},\end{array}$$

where the effective viscosity, η, is

$$\begin{array}{}\text{(4)}& \mathit{\eta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{-\mathrm{1}/n}{\dot{\mathit{ϵ}}}_{\mathrm{e}}^{(\mathrm{1}-n)/n}.\end{array}$$

Here, n is the Glen's flow law exponent, and A is a rheological parameter primarily dependent on ice temperature. Since we assume ice to be isothermal, A is set to a constant value in all simulations. The effective strain rate, ${\dot{\mathit{ϵ}}}_{\mathrm{e}}$, is calculated from the strain rate tensor, $\dot{\mathit{ϵ}}$, as

$$\begin{array}{}\text{(5)}& {\dot{\mathit{ϵ}}}_{\mathrm{e}}=\sqrt{\text{tr}\left({\dot{\mathit{ϵ}}}^{\mathrm{2}}\right)/\mathrm{2}}.\end{array}$$

2.1.1 Boundary conditions

The upper surface, $z=z_{\mathrm{s}}(x,y,t)$, evolves subject to

$$\begin{array}{}\text{(6)}& \left({\displaystyle \frac{\partial }{\partial t}}+\mathit{u}\cdot \mathbf{\nabla }\right)(z-z_{\mathrm{s}})={\dot{a}}_{\mathrm{s}},\end{array}$$

where ${\dot{a}}_{\mathrm{s}}$ is the accumulation rate at the upper ice shelf surface. The lower surface, $z=z_{\mathrm{b}}(x,y,t)$, evolves subject to

$$\begin{array}{}\text{(7)}& \left({\displaystyle \frac{\partial }{\partial t}}+\mathit{u}\cdot \mathbf{\nabla }\right)(z-z_{\mathrm{b}})={\dot{a}}_{\mathrm{b}},\end{array}$$

where ${\dot{a}}_{\mathrm{b}}$ is the melt rate at the ice shelf base. Furthermore, the grounded portion is constrained by the condition

$$\begin{array}{}\text{(8)}& b(x,y)\le z_{\mathrm{b}}(x,y,t)\le z_{\mathrm{s}}(x,y,t),\end{array}$$

where b(x,y) is the bed. The surface accumulation rate, ${\dot{a}}_{\mathrm{s}}$ = 1.2 m a⁻¹, 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, ${\dot{a}}_{\mathrm{b}}$, beneath the ice shelf is defined as a function of ice thickness, H, based on the parameterisation used in Favier et al. (2016),

$$\begin{array}{}\text{(9)}& {\dot{a}}_{\mathrm{b}}=\left\{\begin{array}{ll}\mathrm{0},& \text{where ice is grounded, and}\\ \frac{\mathrm{1}}{\mathrm{50}}{H}^{\mathit{\alpha }}\tanh \left(\frac{|\mathit{x}-{\mathit{x}}_{\mathrm{g}}|}{\mathrm{100}}\right),& \text{where ice is floating},\end{array}\right.\end{array}$$

where α is a tuning parameter, and $|\mathit{x}-{\mathit{x}}_{\mathrm{g}}|$ is the distance to the grounding line. During computation, x represents the position of the current node, and x_g represents the position of the grounding line node closest to the current node. Here, x and x_g are expressed in kilometres, H is expressed in metres, and ${\dot{a}}_{\mathrm{b}}$ has units of metres per year.

A constant flux of $Q{|}_{x=\mathrm{0}}$ = 5.4 × 10⁹ m³ a⁻¹ into the domain is prescribed at the upstream boundary, corresponding to an initial velocity of 300 m a⁻¹. Ice flows through a fixed calving front, where ice is subject to sea pressure. At the lateral boundaries, a free-slip condition is applied, and the flow velocity is subject to the Dirichlet boundary condition $\mathit{u}\cdot \mathit{n}=\mathrm{0}$, where n is the normal vector pointing outwards.

Ice in contact with the bed is subject to a non-linear Weertman-type friction law,

$$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{\mathrm{b}}=-C\mathrm{|}{\mathit{u}}_{\mathrm{b}}{\mathrm{|}}^{m-\mathrm{1}}{\mathit{u}}_{\mathrm{b}},\end{array}$$

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 et al., 2009). The continuous first floating Elmer/Ice grounding line numerical implementation is used (Gagliardini et al., 2016) and was chosen because a discontinuity in basal friction at the grounding line caused undesired numerical artefacts in the ice surface.

Figure 1 The 60 km × 60 km model domain is shown in the case of (a) an ice rise and (b) an ice rumple. A corresponding bird's eye view in (c) and (d) shows the surface velocity magnitude colour-coded and the ice flow direction with arrows. Corresponding along-flow cross-sections are shown in (e) and (f). Sea level is at an elevation of 0 m in the case of the ice rise and 80 m in the case of the ice rumple. In (c), (1), (2), and (3) indicate cross-sections used for analysis, and D is the ice rise dome. Both (2) and (3) 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.

2.2 Idealised model domain set-up

The evolution of ice rises and ice rumples is simulated in a 60 km × 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_{\mathrm{0}}+b_{\mathrm{a}}$, where b₀ is a constant, and b_a is an anomaly with a flat top, defined as

$$\begin{array}{}\text{(11)}& b_{\mathrm{a}}(x,y)=M\exp \left\{{\displaystyle \frac{-\left(\right(x-x_{\mathrm{0}}{)}^{\mathrm{2}}+(y-y_{\mathrm{0}}{)}^{\mathrm{2}}{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{4}}}}\right\}.\end{array}$$

Table 1 List of parameters used in the simulations.

The centre of the bed anomaly is located at (x₀,y₀), σ controls the horizontal extent, and M is the amplitude of the bed anomaly. 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.

Figure 2 Shown is the mesh resolution. In the horizontal, the mesh is unstructured and has a resolution of 350 m in the area surrounding the ice rise or ice rumple. The background resolution (of the surrounding ice shelf) is 2000 m. The mesh is extruded in the vertical with 10 layers. Note that the geometry is exaggerated by a factor of 30 in the vertical direction.

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 (2022). 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 glacial–interglacial timescales 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.

2.3 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 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 $\mathit{u}={\mathit{u}}_{\mathrm{d}}+{\mathit{u}}_{\mathrm{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

$$\begin{array}{}\text{(12)}& {\mathit{u}}_{\mathrm{d}}=-\mathrm{2}\left({\mathit{\rho }}_{\mathrm{i}}g{)}^n\mathbf{\nabla }{z}_{\mathrm{s}}\mathrm{|}\mathbf{\nabla }{z}_{\mathrm{s}}{\mathrm{|}}^{n-\mathrm{1}}\underset{b}{\overset{z}{\int }}A\right(T^{\prime }\left)\right(z_{\mathrm{s}}-\stackrel{\mathrm{‾}}{z}{)}^n\mathrm{d}\stackrel{\mathrm{‾}}{z}.\end{array}$$

We compare the velocity components only at the surface of the ice and also assume that temperature is constant, and so Eq. (12) reduces to

$$\begin{array}{}\text{(13)}& {\mathit{u}}_{\mathrm{d}}(x,y,z_{\mathrm{s}})=-{\displaystyle \frac{\mathrm{2}A({\mathit{\rho }}_{\mathrm{i}}g{)}^n}{n+\mathrm{1}}}(z_{\mathrm{s}}-z_{\mathrm{b}}{)}^{n+\mathrm{1}}\mathrm{|}\mathbf{\nabla }{z}_{\mathrm{s}}{\mathrm{|}}^{n-\mathrm{1}}\mathbf{\nabla }{z}_{\mathrm{s}}.\end{array}$$

The x-direction and y-direction basal sliding components take the form

$$\begin{array}{}\text{(14)}& {\mathit{u}}_{\mathrm{b}}(x,y)=-C_{\mathrm{b}}\left({\mathit{\rho }}_{\mathrm{i}}g\right(z_{\mathrm{s}}-z_{\mathrm{b}}){)}^{p-q}\mathrm{|}\mathbf{\nabla }{z}_{\mathrm{s}}{\mathrm{|}}^{p-\mathrm{1}}\mathbf{\nabla }{z}_{\mathrm{s}},\end{array}$$

where C_b is the basal friction coefficient and relates to the full Stokes basal friction coefficient, C, as follows:

$$\begin{array}{}\text{(15)}& C_{\mathrm{b}}={\displaystyle \frac{N_{\mathrm{b}}}{C^{\mathrm{1}/m}}}\end{array}$$

with

$$\begin{array}{}\text{(16)}& N_{\mathrm{b}}={\mathit{\rho }}_{\mathrm{i}}g(z_{\mathrm{s}}-z_{\mathrm{b}}),\end{array}$$

where N_b=N_be_z is the basal normal stress. In Eq. (14), p and q are chosen for consistency with the non-linear Weertman-type friction law described above.

2.4 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

$$\begin{array}{}\text{(17)}& \mathbf{\nabla }\cdot \mathit{Q}={\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial }{\partial R}}\left(R{\mathit{Q}}_{\mathrm{R}}\right)={\dot{a}}_{\mathrm{s}},\end{array}$$

assuming no azimuthal variance. Here, Q=Q_Re_R denotes the vertically integrated flux at a distance R from the origin (located at the ice rise divide). The resulting ice geometry profile is of the form

$$\begin{array}{}\text{(18)}& h\left(R\right)=h_{\mathrm{0}}{\left[\mathrm{1}-{\left({\displaystyle \frac{R}{L}}\right)}^{\frac{n+\mathrm{1}}{n}}\right]}^{\frac{n}{\mathrm{2}n+\mathrm{2}}},\end{array}$$

where

$$\begin{array}{}\text{(19)}& h_{\mathrm{0}}={\mathrm{2}}^{\frac{n}{\mathrm{2}n+\mathrm{2}}}{\left({\displaystyle \frac{{\dot{a}}_{\mathrm{s}}}{\mathrm{2}{A}_{\mathrm{0}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}n+\mathrm{2}}}{L}^{\frac{n+\mathrm{1}}{n}}\end{array}$$

and

$$\begin{array}{}\text{(20)}& A_{\mathrm{0}}={\displaystyle \frac{\mathrm{2}A({\mathit{\rho }}_{\mathrm{i}}g{)}^n}{n+\mathrm{2}}}.\end{array}$$

L is the horizontal distance from the ice rise divide to the grounding line, and both h₀ 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 (3) in Fig. 1c).

Figure 3 The change in sea level for the transient simulations. The low and intermediate scenarios follow the green curve. A second sea level increase-and-decrease cycle is performed for the low-friction scenario (blue). Sea level is increased to 170 m in the high-friction scenario, and a separate sea level decrease branch is simulated from 155 m (red curve). Sea level is increased and decreased at rates of ± 0.02 m a⁻¹. The crosses indicate points in the low-friction scenario at which a steady-state branch is started with constant sea level in order to compare to the transient simulation.

Figure 4 A cross-section of the ice rise in the along-flow direction for (a) the low-basal-friction, (b) the intermediate-basal-friction, and (c) the high-basal-friction scenarios. The contours show lines of equal velocity (full Stokes) in the x direction, i.e. in the along-flow direction.

2.5 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⁶, C = 7.624 × 10⁶, and C = 3.812 × 10⁸ $\mathrm{Pa}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}$, which we refer to as low-, intermediate-, and high-friction scenarios, respectively. The intermediate-friction scenario has the same basal friction coefficient as that used in MISMIP (Pattyn et al., 2012) and in Favier and Pattyn (2015), where an ice rise is also modelled. The low-basal-friction coefficient is close to the suggested value of 3.16 × 10⁶ $\mathrm{Pa}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}$ in MISMIP+ (Cornford et al., 2020). The high-basal-friction scenario essentially excludes basal sliding, mimicking ice frozen to the bed. For each basal friction coefficient, transient simulations with variable sea level are performed (Fig. 3). In the low- and intermediate-basal-friction scenarios, sea level is increased by 80 m at a rate of 0.02 m a⁻¹ over 4000 years and then stays constant for another 2000 years. Sea level is then decreased at a rate of 0.02 m a⁻¹ back to the initial level followed by a second phase of constant sea level for 2000 years. A second cycle is performed for the low-basal-friction scenario for a comparison with the first cycle. Branches of the low-basal-friction simulation are run to steady state (equilibrium) at discrete intervals while keeping sea level fixed. We run these simulation branches in order to understand how far from steady state the transient simulations are. The choice of sea level perturbation rate is in line with observations, showing periods of sea level rise of up to 0.04 m a⁻¹ during the last deglaciation (Deschamps et al., 2012). Furthermore, we run branches of the simulation beyond the original sea level at the same sea level decrease rate of 0.02 m a⁻¹.

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 in sea level further at a constant rate of 0.02 m a⁻¹ until the transition occurs. A reversal of the sea level perturbation is performed from a height of 155 m above the initial sea level.

3 Results

3.1 Steady-state analysis before sea level perturbation

After 2000 years of spin-up time, ice rises with a characteristic local flow regime develop in all three full Stokes scenarios (Fig. 4). From low to high friction, they vary in maximum thickness (H_max = 213–468 m), grounded area (132–225 km²), and characteristic timescale (t_c = 178–391 a, defined as $t_{\mathrm{c}}=H_{max}/\dot{a_s}$). The characteristic timescale is a metric that gives an indication of the rate of development of Raymond arches (Martín et al., 2009; Goel et al., 2020). The ice divide position in the high-friction scenario has a stossward offset of 0.9 km from the vertical symmetry axis of the bed protrusion. In the intermediate- and low-friction scenarios it is shifted stossward by 2.7 and 3.3 km, respectively. In all three cases, there is substantially more grounding on the stoss side of the bed protrusion than on the lee side.

Topographic and flow divides coincide in all three cases, and ice rise surface velocities are within tens of metres per year. There is negligible basal sliding in the high-friction scenario (with average absolute velocities of roughly 0.5 × 10⁻⁴ m a⁻¹ at the bed–ice interface), whereas basal sliding in the along-flow cross-section (Label (3) in Fig. 1c) accounts for 90 % and 98 % of the local mean horizontal velocities in the intermediate- and low-friction scenarios, respectively. The width of the lateral shear zones, here defined as the lateral distance from the grounding line along a cross-section (Label (2) in Fig. 1c) in which v_x reaches 90 % of v_x at the domain boundary, varies marginally from 10 to 11.3 km from the low- to the high-friction scenarios. Ice fluxes upstream of the protrusion are approximately equal, but mean velocities are 15 % slower, and ice is about 15 % thicker in the high-friction scenario compared with the low-friction scenario.

All ice rises exhibit geometries and flow regimes which are comparable to observations. For example, the high-friction scenario is comparable to Derwael Ice Rise, where previous studies have assumed no basal sliding a priori (e.g. Drews et al., 2015). Basal sliding in the low- and intermediate-basal-friction scenarios means these ice rises are more susceptible to transition into ice rumples when sea level is raised, as shown later.

Figure 5 A bird's eye view of the grounded area corresponding to the steady states at t = 2000 years of the simulations with (a) a low-basal-friction, (b) an intermediate-basal-friction, and (c) a high-basal-friction scenario. In colour, the percentage difference is shown between the calculated SIA surface velocity magnitude and the full Stokes velocity magnitude.

Figure 6 The full Stokes and SIA surface velocities at t = 2000 years in the along-flow direction (a, c, and e) and in the across-flow direction (b, d, and f), as indicated by the cross-sections A–A^′ and B–B^′ through the divide in Fig. 5. Panels (a) and (b) show the low-basal-friction scenario, (c) and (d) show the intermediate-basal-friction scenario, and (e) and (f) show the high-basal-friction scenario.

Figure 7 Cross-sections of the full Stokes simulations at t = 2000 years in the case of the (a) low-, (b) the intermediate-, and (c) the high-basal-friction coefficient. Using a reference point on the ice surface at the grounding line, a Vialov profile is calculated and plotted.

The comparison of full Stokes surface velocities to SIA surface velocities on ice rises illustrates where the local flow assumptions are violated. Figure 5 shows that all three basal friction scenarios have mismatches near the ice rise divides where longitudinal stress gradients are significant. The high-basal-friction scenario shows a good fit otherwise, as do the low- and intermediate-basal-friction scenarios on the lee sides. However, for these cases, surface velocities differ more on the stoss sides of the ice rises (Fig. 6). In the low-friction scenario, absolute deviations increase from 0 %–20 % in the vicinity of the divide (but not at the divide) to over 100 % closer to the grounding line. In the intermediate-friction scenario, deviations are not quite as significant but nonetheless reach deviations of 100 %. In terms of ice thickness, the Vialov approximation captures the high-friction scenario well despite the non-flat bed, while it significantly overestimates the low- and intermediate-basal-friction scenarios in which basal sliding is dominant (Fig. 7).

Figure 8 The response of grounded area and upstream ice shelf velocity to sea level perturbation in the case of low basal friction. Panels (a) and (c) show the evolution for the first sea level increase-and-decrease cycle in blue and red. Panels (b) and (d) show the evolution for the second increase-and-decrease cycle. These curves are also plotted in panels (a) and (c) with dashed red and blue lines for comparison. The crosses represent the results of steady-state branches of the transient simulations at corresponding sea levels. The transition from ice rise to rumple and vice versa is represented by the black dots and a change in colour of the curve.

Figure 9 The response of grounded area and upstream ice shelf velocity to sea level perturbation in the case of high basal friction. In (a), the grounded area is plotted against sea level displacement, and in (b), the average velocity in the x direction in a cross-section upstream of the ice rise (at 20 km from the influx boundary) is plotted against sea level displacement. Red indicates that the system exhibits a characteristic flow regime of an ice rise, and blue indicates that of an ice rumple. The square indicates from where a reversal of the sea level perturbation is simulated.

Figure 10 The response of the dome position to a raising and lowering of sea level in the case of (a) the low-, (b) the intermediate-, and (c) the high-basal-friction coefficients.

3.2 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), lee-side grounding line position, and the upstream ice shelf velocities (Figs. 8c, d and 9b). The upstream ice shelf velocity is defined as the mean velocity of ice in the x direction at x = 20 km, as marked by Label (1) in Fig. 1c. In terms of these metrics, the low- and intermediate-basal-friction scenarios behave distinctly differently than the high-basal-friction scenario. The former transition gradually to ice rumples if sea level is raised past a certain threshold and regrow into ice rises if sea level is reversed. The reversal is not symmetric, and the respective steady-state geometries depend on the history of their evolution (i.e. hysteresis). The high-basal-friction scenario, on the other hand, requires a much larger sea level perturbation to trigger transition into an ice rumple. Once this transition is reached, an ice rumple forms, but the system is unstable, and the ice rumple ungrounds entirely. Details of these differing states are provided in the following.

Before transitioning to an ice rumple, the dome position in the low-friction scenario migrates linearly at a rate of 1.7 m a⁻¹ with increasing sea level (Fig. 10). The dome of the intermediate-friction ice rise migrates first at a rate of 0.8 m a⁻¹ before increasing to a migration rate of 5.7 m a⁻¹ after a sea level increase of 29 m. The dome of the high-basal-friction ice rise exhibits a slow response to sea level displacement during the first 152 m of sea level increase with a divide migration rate of 0.2 m a⁻¹, before increasing to a migration rate of 5.0 m a⁻¹.

After a sea level increase of 20 m in the case of the low-friction case, 30 m in the case of the intermediate-friction case, and 161 m in the high-friction case, the grounding line on the lee side of the ice rise migrates past the highest point of the bed anomaly (marked by A in Fig. 1e) and so is located on a retrograde slope. A transition from ice rise to ice rumple occurs at a further sea level displacement of 30, 16, and 1 m after the grounding line has reached this point in the case of the low-, intermediate-, and high-basal-friction scenarios, respectively.

A steady acceleration of the upstream ice shelf is seen in both the low- and intermediate-basal-friction scenarios, and there 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.

Figure 11 The figures show a cross-section of the ice rises in the along-flow direction for (a) the low-, (b) the intermediate-, and (c) the high-basal-friction scenario. The dotted lines show the geometry of the ice rises before sea level perturbation, and the solid lines show the geometry after a full cycle of sea level increase and decrease.

Figure 12 An along-flow cross-section of the ice rumple at t = 8000 years in the case of (a) the low-basal-friction and (b) the intermediate-basal-friction scenario. The contours show lines of equal velocity in the x direction and are spaced 25 m a⁻¹ apart.

After keeping the sea level constant for 2000 years at a sea level perturbation of 80 m, the low- and intermediate-basal-friction ice rumples evolve to their respective steady states, with minimum velocities of 20 and 38 m a⁻¹ (Fig. 12). Reversal of the sea level perturbation then triggers an asymmetric reversal of the variables of interest described above, with grounded area and upstream ice shelf thickness increasing and upstream velocity decreasing. A transition from ice rumple to ice rise (Figs. 8 and A1) is observed when sea level is 21 and 19 m above the initial sea level in the low- and intermediate-basal-friction scenarios, respectively (as opposed to displacements of 50 and 45 m for low- and intermediate-basal-friction scenarios in the sea-level-increase scenarios, respectively). Once the original sea level is again reached, the ice rises in both the low- and intermediate-basal-friction scenarios are smaller, with a smaller grounded area and a lesser buttressing effect on the upstream ice shelf (Figs. 8, A1, and 11). The upstream ice shelf in the case of the low-basal-friction scenario has a decrease in velocity of 18 m a⁻¹, whereas the ice shelf in the intermediate-basal-friction scenario decreases in velocity by 25 m a⁻¹. A second cycle of sea level increase 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 and ice shelf velocity is calculated as described above and presented in Fig. 8. The hysteresis cycle is now closed, with the final steady state corresponding to the state before the last sea level perturbation cycle.

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), the ice rise volume and grounded area also recover asymmetrically, resulting in two differing states for a given sea level displacement (Fig. 9).

We investigate the migration of the stoss- and lee-side grounding lines of the ice rise and make a comparison with the grounding line position in the case of hydrostatic equilibrium (video in the Supplement). The maximum differences in position are 0.5 km on the stoss side and 0.4 km on the lee side, with mean differences of 0.2 km in both cases. During sea level increase, the hydrostatic grounding line positions have a delayed response in comparison with the Elmer/Ice grounding line. On the other hand, during sea level decrease, the hydrostatic grounding lines have a more rapid response.

4 Discussion

4.1 The influence of basal sliding on the geometry and transient behaviour 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 the freezing point at the ice–bed interface near the summits (Martín et al., 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 simulations show that ice rises can form in scenarios where basal sliding is significant. Surface velocities in the low and intermediate scenarios are within a few metres per years near the crests, similar to the predictions in the high-friction scenario (Fig. 5). In this regard, surface velocities alone are a poor indicator for the presence or absence of basal sliding on ice rises. However, the geometries between the three scenarios differ significantly, and only the high-friction scenario can be adequately approximated with the Vialov profile, whereas the low and intermediate scenarios exhibit significant misfits (Fig. 8). This means that a simple fit with a Vialov profile can serve as a first-order metric for the absence or existence of basal sliding for specific ice rises. This is important, as the degree of basal sliding in the vicinity of the grounding line determines the local ice flow and the ice rise's transient behaviour in response to sea level perturbation. When comparing the grounding line positions of the full Stokes model and the hydrostatic grounding line position, we find that differences are small. However, over the millennial timescales considered here, together with the compounding effect of the small errors in grounding line position at each time step, it is possible that a hydrostatic assumption may result in differing ice rise and ice rumple geometries as well as a differing transition point.

Many ice rises are fully surrounded by ice shelves, and the extent to which isle-type ice rise velocities are affected by longitudinal and shear stresses transferred from the upstream ice shelf is not fully clear. This effect is analysed here using the differences between the non-local full Stokes simulations and the fully local SIA. The flow regime in the high-friction scenario is, to a large extent, independent of the surrounding ice shelf. In the low- and intermediate-basal-friction scenarios, however, the differences between full Stokes and SIA are greater and are especially evident on the stoss side of the ice rise. 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, 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.

Figure 13 An along-flow ground-based radargram (Drews, 2019) showing an ice rumple in the Roi Baudouin Ice Shelf, East Antarctica, is shown in (a). The flow of ice is from left (A) to right (A^′). In (b), the location of the radargram (A–A^′) is shown (Jezek, 2003).

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 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 metres 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⁻¹ are, however, significantly faster than the example observed at the Roi Baudouin Ice Shelf, where the ice is effectively stagnant (Berger et al., 2016). The smooth transition of the low- and intermediate-friction ice rises into ice rumples reflects their strong 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.

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 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.

Interestingly, the low-friction ice rumple exhibits lower minimum velocities than the intermediate-friction ice rumple, most 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 horizontal resolution applied.

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 – 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 less stable for an equal amount of sea level displacement.

4.2 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 and 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 A1b in the Appendix).

There is a difference in the individual pairs, with the grounded area being larger in the sea-level-increase scenario than in the sea-level-decrease scenario. 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. 9) 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 of drag (Still et al., 2019) and consequently a differing net resistance to the upstream ice shelf.

A self-stabilising feedback occurs, with divide migration opposing grounding line retreat in a sea-level-increase scenario. The ice rise height reduces, and the divide migrates stossward during lee-side grounding line retreat. Because the divide moves stossward, the area of accumulation adjacent to the divide on the lee side of the ice rise increases. The increased accumulation area promotes an increased flux across the grounding line, opposing grounding line retreat. Analogously, sea level decrease results in leeward divide migration. The resulting reduction in accumulation area adjacent to the divide on the lee side of the ice rise opposes grounding line advance. The existence of negative feedback mechanisms in both the sea-level-increase and sea-level-decrease scenario results in hysteretic behaviour (Figs. 8, 9 and A1).

Another mechanism that plays a role is the sensitivity of the grounding line to 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 the hysteretic behaviour. In future work it is worth investigating whether effects such as an increased melt rate also produce a 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, 2020) does provide some information in this regard but also 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 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 propagation during transient simulation, leading to inaccurate grounding line migration if a hydrostatic equilibrium assumption is used.

5 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 steady-state full Stokes velocities with steady-state SIA velocities, whereas in a low-basal-friction scenario the mismatch is larger due to stronger mechanical coupling to 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 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 importance of perturbation history for the formation of the correct feature. 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 also likely to result in hysteretic and potentially irreversible behaviour in ice shelf buttressing upstream of ice rises.

Appendix A: The response of the grounded area and ice shelf velocity to sea level perturbation in the intermediate-basal-friction scenario

Presented in Fig. A1 is the response of the grounded area and upstream velocity to sea level perturbation in the case of the intermediate-basal-friction scenario (C = 7.624 × 10⁶ $\mathrm{Pa}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}$). The transition from ice rise to ice rumple occurs at a sea level displacement of 19 m, and the transition from ice rumple to ice rise occurs at a sea level displacement of 45 m, compared with 21 and 50 m, respectively, in the low-basal-friction scenario. Interestingly, the grounded area of the ice rumple follows a rather linear path in the intermediate-basal-friction scenario compared with the low-basal-friction scenario.

Figure A1 The response of the grounded area and ice shelf velocity to sea level perturbation in the intermediate-friction scenario. In (a), the grounded area is plotted against sea level displacement, and in (b), the average velocity in the x direction in a cross-section upstream of the ice rise (at 20 km from the influx boundary) is plotted against sea level displacement. Red indicates that the system exhibits a characteristic flow regime of an ice rise, and blue indicates that of an ice rumple.

Appendix B: Comparison between the first floating and discontinuous grounding line numerical implementations

In the case of the low-basal-friction scenario, we have run equivalent simulations using a differing grounding line numerical implementation, namely the discontinuous method (Fig. B1). At the grounding line, basal friction is applied if the other two nodes in the element are also grounded, and a free-slip condition is applied if the other two nodes are ungrounded. The first floating numerical implementation, however, assumes a free-slip condition at the grounding line, and a linear reduction in basal friction between it and the upstream node is applied. Although the discontinuous numerical implementation has been shown to have the least dependence on mesh resolution, it can be argued that the first floating is more plausible physically, with effective pressure disappearing at the grounding line (Gagliardini et al., 2016). The simulations show that regardless of the numerical implementation, hysteresis occurs.

Figure B1 Shown is the response of the grounded area in the low-friction case of the first floating (red) and discontinuous (blue) Elmer/Ice numerical grounding line implementations.