2. Observations of glacier persistence

Central West Greenland (CWG) provides a particularly well-observed laboratory for understanding glacier retreat over bumpy beds. As in most of Greenland, surface melting has been persistently intensifying since the 1970s (or potentially earlier; Trusel and others, Reference Trusel2018). In the late 1990s an influx of warm water from the North Atlantic arrived in glacier fjords in this region (Holland and others, Reference Holland, Thomas, De Young, Ribergaard and Lyberth2008). A compilation of terminus positions recorded by visible satellite imagery (Catania and others, Reference Catania2018) show that many glaciers in CWG retreated between the late 1990s and the early 2000s when ocean temperatures were warm (Fig. S1). However, some glaciers in this region have not retreated during the observational record. Figure 1a shows observations of terminus positions at four such persistent glaciers, Kangerdlugssup Sermerssua (blue), Kangerluarsuup Sermia (green), Sermeq Kujalleq (purple, aka Store Gletscher) and Sermeq Avannarleq (orange). Figure 1b shows the along-flow bed topography at these same glaciers from the BedMachine v3 dataset (Morlighem and others, Reference Morlighem2017). These glaciers have persisted <1 km downstream of bed peaks, indicating the critical importance of peaks in bed topography in potentially delaying or preventing rapid glacier retreat. This persistence likely cannot be explained by a heterogeneity in warm water reaching the termini of these glaciers, since warm water was pervasive below 150 mdepth throughout this region (Holland and others, Reference Holland, Thomas, De Young, Ribergaard and Lyberth2008), and all of these glaciers are connected to the continental shelf through deep bathymetric troughs (Fig. S2; though submarine bathymetry is subject to uncertainty, see Fig. S3). Glaciers in this region that did retreat following the ocean warming event in the late 1990s mostly retreated away from bed peaks (on which they had previously persisted) and some have since ceased retreat upon reaching a new bed peak (Fig. S1; Catania and others, Reference Catania2018). In Figure 1, we show one example of such a glacier, Umiammakku Sermia (red), which began rapidly retreating away from a bed peak $\sim$5 years after the arrival of warm waters in the region, before ceasing retreat at a different bed peak around 2010.

Fig. 1. Observational evidence of terminus and terminus persistence at bed peaks in Central West Greenland (CWG). (a) Terminus positions ($x$-axis) over time ($y$-axis) at five CWG glaciers derived from satellite-based sensors (Catania and others, Reference Catania2018). (b) Along-flow bed topography at CWG glaciers in panel (a), with the $x$-axis is the along-flow distance relative to recent (2016) terminus position, with $x = 0$ representing the present position of the glacier termini and gray shading indicating where there is currently grounded glacier ice ($x< 0$). Nearest bed peaks upstream of the current terminus denoted by a filled circle in each case. For glaciers with strong cross-fjord variations in topography (Kangerluarsuup, Kujalleq), the deepest bed topography across the fjord is used; for the others, mean topography across the fjord is used. Bed topography error range is plotted in Figure S3, with typical error for proglacial fjords of $< 10$ m and typical error for near-terminus subglacial topography of 10–100 m. Bathymetry from BedMachine data compilation (panel b) (Morlighem and others, Reference Morlighem2017). (c) Location of CWG glaciers in panels (a) and (b) on polar stereographic north projection (EPSG:3413). (https://epsg.io/3413).

Geological evidence from regions of past glacier retreat further demonstrates the importance of bed peaks in the response of glaciers to climate change. The bathymetry of the Ross Sea, Antarctica is composed of smooth, flat troughs separating large plateaus. Amid this smooth bathymetry, localized recessional moraines and grounding zone sediment wedges record locations where the deglacial retreat of glaciers in the Ross Sea embayment paused for prolonged time periods (Simkins and others, Reference Simkins2017; Greenwood and others, Reference Greenwood, Simkins, Halberstadt, Prothro and Anderson2018). Figure 2 focuses on three particular locations in the Ross Sea where high-resolution multibeam bathymetric observations show pervasive grounding zone wedges connecting, parallel to, and on top of isolated seamounts on otherwise flat topography. Grounding zone wedges are formed when the grounding line of a glacier persists at a location for a sufficiently long time period that a significant sediment deposit forms at that location (Dowdeswell and Fugelli, Reference Dowdeswell and Fugelli2012). These grounding zone wedges are not present in surrounding portions of the seafloor in the Ross Sea, indicating that bed peaks (which are generated by non-glaciological processes) at both of these locations in the Ross Sea exerted an important control on glacier retreat. Such grounding line persistence at isolated bed peaks also occurs in 3D numerical ice flow simulations of the retreat of West Antarctic ice streams through the Ross Sea following the last glacial period (e.g. Golledge and others, Reference Golledge2014). Other marine geophysical surveys of the seafloor in regions of past glacier retreat also reveal widespread geological evidence for prolonged periods of terminus persistence at bed peaks over a wide range of time periods and local conditions (Stoker and others, Reference Stoker, Bradwell, Howe, Wilkinson and McIntyre2009; Todd and Shaw, Reference Todd and Shaw2012; Greenwood and others, Reference Greenwood, Clason, Nyberg, Jakobsson and Holmlund2017).

Fig. 2. Bathymetry (in meters below sea level, m b.s.l.) of the southwestern Ross Sea (inset) with linear to sub-linear grounding zone wedges (i.e.paleo-grounding lines) concentrated on and between isolated volcanic seamounts. (a) Paleo-grounding lines positioned (indicated by brown lines) between clustered, flat-topped seamounts. (b) Pinning of a paleo-grounding line (brown line) on and around a seamount. (c) Paleo-grounding line (indicated by a pointer) pinning on and between seamounts that are separated by several kilometers. Multibeam echo sounding bathymetry was collected on cruise NBP15-02A (Simkins and others, Reference Simkins2017; Greenwood and others, Reference Greenwood, Simkins, Halberstadt, Prothro and Anderson2018).

3. Model description

To simulate a typical marine-terminating glacier near a bed peak, we use a one-dimensional flowline model (in the x-direction) of a marine-terminating glacier (with the same equations and numerical methods as Schoof, Reference Schoof2007a, and many other studies). We simulate the glacier velocity ($u$), thickness ($h$) and terminus position ($x_t$). Velocity is determined from the shallow stream/shelf approximation (SSA) of the momentum balance

(1)$$\eqalign{{\partial\over \partial x} \left(2\bar{A}^{-1/n} h \left\vert {\partial u\over \partial x} \right\vert ^{1/n - 1} {\partial u\over \partial x} \right)- C\vert u\vert ^{m-1} u - \rho_{\rm i} gh {\partial \over \partial x} \left(h-b \right) = 0,\;} $$

where $\bar {A}$ is the depth-averaged flow law rate factor, $C$ is the sliding law coefficient, $m$ is the sliding law exponent, $\rho _{\rm i}$ is the density of ice, $g$ is gravity, $b$ is the bed topography. Ice thickness is determined from the mass conservation equation

(2)$${\partial h\over \partial t} + {\partial\over \partial x} \left(uh \right) = a,\; $$

where $a$ is the surface mass balance (SMB, net annual snowfall and surface melt). In the cases we consider in this study, ice flow is dominated by sliding over a moderately slippery bed (parameters listed in Table 1), and the role of lateral shear stresses are not considered (in order to aid comparison with established theory). At the ice divide ($x = 0$), we set $u = 0$ and ${\partial ( h-b) \over \partial x} = 0$. We assume that the terminus is just at flotation (i.e. an unbuttressed grounding line), so we prescribe

(3)$$h( x_{\rm t}) = -{\rho_{\rm w}\over \rho_{\rm i}} b( x_t) ,\; $$

where $\rho _{\rm w}$ is the density of seawater and bed topography is negative when below sea level. Ice velocity at the terminus is calculated (not prescribed) from the stress balance, assuming that the driving stress is fully supported by extensional stresses

(4)$${\partial\over \partial x} \left(2\bar{A}^{-1/n} h \left\vert {\partial u\over \partial x} \right\vert ^{1/n - 1} {\partial u\over \partial x} \right)- \rho_{\rm i} \left(1-{\rho_{\rm i}\over \rho_{\rm w}} \right)gh {\partial h\over \partial x} = 0 .$$

Table 1. Parameter values for steady-state and transient retreat simulations (unless otherwise specified in text)

The momentum and mass conservation equations are solved using an implicit backward Euler method with a 1 year time step, and the constraint that the terminus is at flotation. The entire model is formulated on a moving mesh that is specified to be coarse (1–3 km) in the glacier interior and fine (50–100 m) in the grounding zone region, which is specified as the final 3% of the grounded glacier domain. The moving mesh is constructed so that the model domain always encompasses the entirety of the grounded glacier. This refined moving mesh numerical modeling approach has been shown to be a highly accurate method for simulating marine-terminating glaciers (Vieli and Payne, Reference Vieli and Payne2005).

The four idealized bed topographies we consider in this study all have a single sharp bed peak, but with different reverse bed slopes just upstream of the peak, and otherwise the same forward-sloping bed (i.e. shallowing toward the interior). We specify the bed topography through a piecewise-linear function that changes slope at the bed peak and trough. We design such idealized bed topographies to resemble the height and slopes around bed peaks in Greenland (Fig. 1b) and Antarctica (Fig. 2). More extensive high-resolution observations of subglacial bed topography indicate pervasive regions of Greenland and Antarctica where bed slope changes sign on horizontal length scales of less than a kilometer (Morlighem and others, Reference Morlighem2017, Reference Morlighem2020). We consider other versions of the bed peak topography in simulations plotted in the supplementary materials.

In the next two sections, we consider steady-state glacier configurations and transient simulations of glacier evolution. Glacier steady-states are determined by numerically solving for glacier states, with rates of change that are zero to within machine precision. To confirm stability of these steady-states, we perturb the SMB by 0.001 m a$^{-1}$ and solve the transient glacier evolution, ensuring it does not evolve to a very different steady-state. Transient simulations are initialized from a steady-state and perturbed through a step reduction in SMB. In the supplementary materials, we consider various other forms of forcing (e.g. ocean melt and a linear trend in SMB) and find qualitatively similar results to the step reduction in SMB.

The type of model described here has been widely implemented, and shown to compare well to real glaciers in Greenland and Antarctica which flow primarily through sliding and float at their terminus (e.g. Enderlin and others, Reference Enderlin, Howat and Vieli2013; Jamieson and others, Reference Jamieson2014; Huybers and others, Reference Huybers, Roe and Conway2017). We have also replicated the substance of the results described hereafter in simulations with a range of horizontal resolutions (Fig. S4) and high-resolution simulations of the Elmer/Ice Full-Stokes numerical ice-flow model (see Figs S5, S6), indicating that the resolution or SSA approximation does not affect the substance of the conclusions in this study.

4. Enhanced glacier stability at bed peaks

We first consider how bed peaks affect glacier stability over a range of climate forcing. Simulations show (Fig. 3b) that over a wide range of SMB, glacier termini persist indefinitely (i.e. reside at a stable steady-state) near bed peaks. We find that the sharp bed peaks we consider in this study, which entail an instantaneous transition (in space) from a forward-sloping bed to a reverse-sloped bed, lead to glacier stability over a wider range of SMB than what is predicted in prevailing theories of terminus stability (dotted line in Fig. 3b, reproduced from theory of Schoof, Reference Schoof2012). The steeper the reverse sloped bed upstream of the bed peak, the wider the range of SMB over which the glacier will remain stable. For the steepest reverse slope (shown in blue), the glacier remains stable a short distance downstream of the topographic high for a significant range of SMB from 0.5 to 1.0 m a$^{-1}$.

Fig. 3. Simulated stable terminus positions in the vicinity of a bed peak. (a) Four idealized bed topographies with differing bed slope just upstream of bed peak. (b) Bifurcation diagrams showing steady-state terminus positions over a range of surface mass balance. For each value of SMB, an initial guess on either side of the bed peak is used to determine if more than one steady-state exists. Lines are plotted from simulations of steady-state glacier state at 0.001 m a$^{-1}$ increments of SMB (seaward of the bed peak simulated steady-state terminus positions are 1–10 m apart). The dotted line is the stable terminus positions calculated by solving $ax_{\rm g} = Q_{\rm g}( x_{\rm g})$, with analytical approximation for ice flux from Schoof (Reference Schoof2007b) (from topographies in panel a), derived by neglecting the effect of local slope. (c) The fractional difference between the ice flux at the terminus predicted from our numerical solution, $Q_{{\rm num}}$, and the ice flux that would be predicted on neglect of the effect of local slope, $Q_{\rm g}( x_{\rm g})$, as a function of distance from the bed peak (normalized by terminus ice thickness).

Ice must flow uphill to reach a terminus located at a bed peak. This slows the flow of ice approaching the bed peak and decreases the ice flux through a terminus retreating toward a bed peak. Figure 3c shows how ice flow near the terminus is changed by local bed slope by plotting the deviation of simulated steady-state terminus ice flux from the ice flux predicted by theories assuming negligible bed slope near the terminus (i.e. Schoof, Reference Schoof2007b). As the terminus approaches within $\sim 10$ ice thicknesses of the bed peak, the ice flux decreases much more rapidly than is predicted under the assumption of negligible bed slope near the terminus. Indeed the magnitude of this reduction in ice flux near the bed peak (10–50% in these examples) is comparable to the effect of substantial ice shelf buttressing on grounding line ice flux (Reese and others, Reference Reese, Winkelmann and Gudmundsson2018; Mitcham and others, Reference Mitcham, Gudmundsson and Bamber2022). The cause of this rapid decline in ice flux is a lowered surface slope and driving stress and hence lowered velocity on the ice flowing up the bed peak, which then lowers terminus ice velocity, thickness and ice flux through longitudinal viscous stresses. Figure 4 shows surface elevation profiles for three steady-state simulations with termini just downstream of the bed peak, demonstrating how the surface slope flattens where ice flows up to the bed peak. Seen another way, ice that slows as it flows up the bed peak thickens, increasing the gradient in height above buoyancy near the terminus, which reduces the potential ability for changes in SMB or ocean melt at the terminus to affect terminus position.

Fig. 4. Simulated near-terminus surface elevations for three stable glacier configurations near the bed peak for the $b_x = 0.004$ bed topography (red lines in Fig. 3). Gray shading indicates region of reverse-sloping bed.

The enhanced stability of simulated termini near bed peaks, as compared to prior theory, explains counterintuitive aspects of observations. There is a wide range of external forcing over which a terminus will persist at a bed peak, explaining why so many glacier termini are observed at bed peaks on bumpy bed topography. Indeed, repeating these steady-state simulations for corrugated bed topography (a repeated series of peaks and troughs) indicates stable glacier terminus positions exist almost exclusively at bed peaks for retreating glaciers (Fig. S7, though advancing glaciers may have stable configurations away from bed peaks). The reduced glacier sensitivity to climatic changes at bed peaks also explains why many glaciers are observed to persist, seemingly on the precipice of instability, even while experiencing substantial fluctuations in local climate (Tinto and Bell, Reference Tinto and Bell2011; Hillenbrand and others, Reference Hillenbrand2017; Catania and others, Reference Catania2018). Such enhanced stability of glaciers at bed peaks is in contrast to the prevailing idea that glaciers at bed peaks are necessarily ‘vulnerable’ to even small changes in climate due to their spatial proximity to reverse-sloping beds over which the marine ice-sheet instability occurs (Ross and others, Reference Ross2012; Morlighem and others, Reference Morlighem2020). As we have shown here, glaciers can be close to a bed peak, but are only ‘vulnerable’ in the sense in that they may initiate a rapid retreat in response in a large change in climate forcing.

The enhanced range of stability near points of destabilization (in state space) is a hallmark of a ‘crossing-sliding bifurcation’ (di Bernardo and others, Reference di Bernardo, Budd, Champneys and Kowalczyk2008). The system behavior in the vicinity of such bifurcations is different from the canonical ‘saddle-node bifurcation’, which has previously been identified as the route through which grounding lines lose stability (Mulder and others, Reference Mulder, Baars, Wubs and Dijkstra2018; Pegler, Reference Pegler2018). In a saddle-node bifurcation, the loss of stability occurs suddenly in state space, without a change in the sensitivity of the stable branch (i.e. lines in Fig. 3b) to parameter changes upon approach to the bifurcation point. In a crossing-sliding bifurcation, the stable branch instead ‘slides’ along the bifurcation manifold with parameter variation (i.e. the system state remains close to the bifurcation without crossing it). Crossing-sliding bifurcations typically arise in system with non-smooth variations in system properties coinciding with bifurcations. In the case of a glacier retreating toward a sharp bed peak, as the system approaches the bifurcation point due to parameter variation (i.e. SMB or ocean melt) the stable glacier state (Fig. 3b) becomes much less sensitive to parameter variation, before eventually crossing the bed peak and initiating a large change in glacier state. This distinction in the type of bifurcation is important because it leads to much larger jumps in the system state (i.e. ice volume loss) upon crossing the bifurcation. As the forcing changes (i.e. SMB decreases), the onset of rapid ice loss is delayed, leading to a higher rate of ice loss if and when the terminus crosses the bed peak. In Figure 3b, this amounts to the difference between a 25 km retreat in terminus position for a relatively smooth bump (i.e. the dotted line), compared to a retreat of 40–100 km for sharper peaks (purple, yellow, red, blue lines).

5. Distinguishing glacier stability from transient persistence

In transient simulations of terminus retreat over sharp bed peaks (those plotted in Fig. 3a), a glacier is initialized at a steady state with its terminus just downstream of a bed peak, and is then subjected to a 40% step reduction in SMB uniformly over the glacier catchment. Figure 5a shows that some of the simulated glaciers retreat up to, then transiently persist just downstream of the bed peak for a period of time spanning decades to centuries (yellow and red lines), before eventually crossing the bed peak and rapidly retreating over the reverse-sloping bed. There are also cases where there is merely a brief slowdown in the rate of retreat at the bed peak (purple line), and other cases where the persistence continues indefinitely (blue line). We define such an indefinitely persistent case as stable in the same mathematical sense that we do in the previous section, where a system state persists forever with no change in forcing. Similar behaviors of transient and indefinite persistence of glaciers at bed peaks also occur in equivalent full-Stokes simulations of glacier retreat over bed peaks (Figs S5, S6) and in SSA simulations of glacier retreat with different types of forcing and smoothed bed peaks (Figs S8–S10). These transient simulations show that even for SMB values that do not correspond to a stable glacier configuration (plotted in Fig. 3b), there may still be prolonged periods of transient terminus persistence. Longer periods of transient persistence lead to more rapid subsequent retreat, which continues even as the terminus encounters forward-sloping bed topography.

Fig. 5. Simulated terminus retreat in the vicinity of a bed peak. (a) Evolution of a terminus from steady state, in response to an instantaneous 40% reduction in surface mass balance over the glacier catchment (1.1–0.66 m a$^{-1}$), for a variety of upstream bed slopes. Terminus position ($y$-axis) is relative to bed peak location as in Figure 3a. (b) Thinning rate 50 km upstream of terminus in transient simulations. (c) Ice velocity 50 km upstream of terminus in transient simulations.

Glaciers that persist at bed peaks continue to lose mass through thinning upstream of the terminus following changes in climate forcing, as seen in observations of recent thinning upstream of the terminus at persistent glaciers throughout Greenland (Kjeldsen asnd others, Reference Kjeldsen2015; Felikson and others, Reference Felikson2017; Mouginot and others, Reference Mouginot2019; Shepherd and others, Reference Shepherd2020). At persistent glaciers in CWG, this thinning is mostly being driven by negative SMB anomalies, which are largely offset by dynamic thickening bringing ice from upstream portions of glacier catchments (Felikson and others, Reference Felikson2017). Ultimately, this upstream-intensified thinning leads a decrease in ice surface slope and upstream slowing, which is captured in observations of persistent CWG glaciers (Joughin and others, Reference Joughin, Smith, Howat, Scambos and Moon2010) and our simulations (Figs 5b-c). Though such thinning is less than that occurring at retreating glaciers through dynamic thinning, it nonetheless shows that persistence of a glacier terminus is not necessarily indicative of a glacier in mass equilibrium.

We can compare our simulated glaciers which stabilize at bed peaks to those which merely pause at bed peaks to ascertain whether observations of persistent glaciers may provide evidence of their eventual fate. We find that, regardless of whether they ultimately remain at the bed peak or retreat from it, the glaciers we simulate have upstream thinning rates within millimeters per year of each other (Fig. 5b and Fig. S11), and ice velocities within meters per year of each other (Fig. 5c and Fig. S12) while they persist at a bed peak. It would thus be exceedingly difficult to observationally distinguish glaciers that are merely paused from those that have stabilized indefinitely at bed peaks. Other studies have also found that, in realistic simulations of the future retreat of glaciers away from bed peaks, small uncertainties in the observed glacier state, bed topography or the climate forcing produce large uncertainties in the timing of the onset of rapid glacier retreat which is then amplified by the divergence of retreat predictions due to marine ice-sheet instability (Gladstone and others, Reference Gladstone2012; Robel and others, Reference Robel, Seroussi and Roe2019). Ultimately, the delicate balance between advection and thinning at persistent glaciers makes it exceedingly difficult to project retreat of glaciers over bumpy bed topography, and further emphasizes the need for more accurate observational constraints on glacier state and rate of change, bed topography and local climate change.