Dynamic influence of pinning points on marine ice-sheet stability: a numerical study in Dronning Maud Land, East Antarctica

The East Antarctic ice sheet is likely more stable than its West Antarctic counterpart because its bed is largely lying above sea level. However, the ice sheet in Dronning Maud Land, East Antarctica, contains marine sectors that are in contact with the ocean through overdeepened marine basins interspersed by grounded ice promontories and ice rises, pinning and stabilising the ice shelves. In this paper, we use the ice-sheet model BISICLES to investigate the effect of sub-ice-shelf melting, using a series of scenarios compliant with current values, on the ice-dynamic stability of the outlet glaciers between the Lazarev and Roi Baudouin ice shelves over the next millennium. Overall, the sub-ice-shelf melting substantially impacts the sea-level contribution. Locally, we predict a short-term rapid grounding-line retreat of the overdeepened outlet glacier Hansenbreen, which further induces the transition of the bordering ice promontories into ice rises. Furthermore, our analysis demonstrated that the onset of the marine ice-sheet retreat and subsequent promontory transition into ice rise is controlled by small pinning points, mostly uncharted in pan-Antarctic datasets. Pinning points have a twofold impact on marine ice sheets. They decrease the ice discharge by buttressing effect, and they play a crucial role in initialising marine ice sheets through data assimilation, leading to errors in ice-shelf rheology when omitted. Our results show that unpinning increases the sea-level rise by 10 %, while omitting the same pinning point in data assimilation decreases it by 10 %, but the more striking effect is in the promontory transition time, advanced by two centuries for unpinning and delayed by almost half a millennium when the pinning point is missing in data assimilation. Pinning points exert a subtle influence on ice dynamics at the kilometre scale, which calls for a better knowledge of the Antarctic margins.

high-resolution dataset), or not correctly resolving PPw (Rignot et al. (2011)'s velocities in combination with ice/bed geometry from Fretwell et al. (2013)). Each point is addressed with transient simulations run over the next millennium with an ensemble of six different sub-ice shelf melting scenarios in combination with two types of sliding laws, because these parameters are poorly constrained by observations. The 36 resulting simulations give a comprehensive overview of the future evolution of the ice sheet and testify of the importance of including even small pinning points in the observational dataset aimed at modelling 5 purpose.

Input data
Each experiment consists of an initialisation by data assimilation and a subsequent set of transient simulations. The former requires surface velocity, ice thickness, bed elevation, and englacial temperatures. The latter requires ice thickness, bed eleva- 10 tion, initial englacial temperatures, an ice stiffening factor and basal friction coefficient field (the latter two computed by the data assimilation), surface mass balance and basal mass balance of the ice shelves.
The computational domain covers an area of about 40,000 km 2 and is illustrated in Figure 1. Two distinct datasets for flowfield and ice/bed geometry were employed. The standard dataset comprises surface velocities from Rignot et al. (2011) and ice/bed geometry from Fretwell et al. (2013). The high-resolution dataset uses the observations of Berger et al. (2016) on the 15 WRG ice shelf, which account for PPw in both surface velocities and ice/bed geometry (the latter called mBedmap2). These two datasets only differ for the WRG ice shelf and are otherwise identical.
Modelling grounding-line advance as a response to ocean-induced perturbation is very sensitive to sub-ice shelf bathymetry, which is roughly interpolated in our studied domain (Le Brocq et al., 2010) and thus largely uncertain. As a consequence, the water column beneath ice shelves is in places very shallow, which can cause spurious ice-shelf re-grounding. In order to 20 make the bathymetry more coherent with both bed elevation at the grounding line and (unpublished) measurements near the ice-shelf front, we lowered the bathymetry beneath the ice shelves in a two-step procedure. First, we excavated a 250 m thick uniform layer 30 km away from the grounding line, ensuring a smooth connection with the grounded area with a 1-D Gaussian function. The second part of the excavation is based on unpublished bathymetric data collected during a 2011 oceanographic survey (K. Leonard, personal communication, 2012), which shows a deep trough more than 850 m deep cutting through the 25 continental shelf between PPw and Derwael Ice Rise (DIR) (Figure 1). This feature may be the relict of past ice sheet erosion from the WRG ice stream when the grounding line was closer to the continental shelf break (Livingstone et al., 2012). We therefore assume the presence of a narrow trough cutting through the bathymetry beneath the ice shelf linked to the deepest section at the grounding line (yellow line in Figure 1). The across flow excavation uses a 1-D Gaussian-shaped function (its half-width is 15 km based on the ice-stream cross-section extent). Both excavations are included in the standard as well as the 30 high-resolution datasets.
The surface mass balance is taken from Arthern et al. (2006) without temporal variation.

4
The Cryosphere Discuss., doi: 10.5194/tc-2016-144, 2016 Manuscript under review for journal The Cryosphere ERG: East Ragnhild Glacier. We also name a group of two pinning points PPhs located at the front of HB, and the pinning point PPw at the front of WRG.
For the ice-shelf basal mass balance, we applied two melt-rates parametrisations M b1 and M b2 , based on Gong et al. (2014) and Beckmann and Goosse (2003), respectively. The former is a scheme that allows the highest melt rates to follow the grounding-line migration, using a combined function of ice thickness and distance to the grounding line, defined as where H is the ice thickness, and G and A are the grounding line and ambient melt rates, respectively. The value of p decreases exponentially with distance to the grounding line, taking the value of 1 at the grounding line and 0 away from it. α 1 and α 2 are tuning parameters. The M b2 parametrisation is based on the difference between the freezing point of water and ocean temperature near the continental shelf break (as developed in Beckmann and Goosse, 2003). The virtual temperature T f at which the ocean water freezes at the depth z b below the ice shelf is defined as where S o is the ocean salinity (set at 34.5 psu from Schmidtko et al. (2014), confirmed by K. Leonard, personal communication, 2012). The melt rates M b2 are prescribed as where ρ w is the density of water, c p0 the specific capacity of the ocean mixed layer, γ T the thermal exchange velocity, T 0 the 10 ocean temperature (set at -1.5°C from Schmidtko et al. (2014) andK. Leonard, personal communication, 2012), L i the latent heat capacity of ice, ρ i the density of ice (see Table 1 for the value of parameters) and F melt a tuning parameter.
Ice temperature data are provided by a three-dimensional thermo-mechanical model (updated from Pattyn, 2010) and are constant in time. in Cornford et al. (2015)). Data assimilation is performed by a control method that solves the adjoint system of equations, as described in Appendix B1 of Cornford et al. (2015). We applied one kilometre resolution at the grounding line during transient simulations (Table 2.1). The relationship between stresses and strain rates is given by the Glen's flow law: where S is the deviatoric stress tensor,ε is the strain rate tensor, η is the effective viscosity (depending on ice temperatures 10 and effective strain rate), and φ is a stiffening factor representing non-thermal viscosity effects, such as crevasse-weakening and ice anisotropy. The basal friction between the grounded ice sheet and the bed is governed by a Weertman-type sliding law (Weertman, 1957): where τ b is the basal traction, C is the friction coefficient, m is the friction exponent and u b is the basal velocity. Initial fields 15 of C and φ were inferred with the control method applied to ice/bed geometry and surface velocities, using the same procedure as in Berger et al. (2016).

Description of the experiments Initialisation
Three sets of initialisations with both linear (m = 1) and nonlinear (m = 1/3) sliding were performed for C, φ (both inferred 20 with the control method), and the initial ice/bed geometry: -B e /S: The control method and the transient simulations use the high-resolution dataset, so that PPw is used for model initialisation and evolution.
-B e /U : This is a variant of B e /S in which transient simulations start from bed elevation and ice thickness without resolving PPw -we use Bedmap2 instead of mBedmap2 -in order to simulate unpinning. Because there is no friction beneath ice shelves, we set the value of the friction coefficient C in case of further ice-shelf regrounding at 500 P a m −1 a. This number causes high basal sliding (comparable to sliding beneath ice streams), which reflects the idea of a sediment-filled bathymetry, and is motivated by the sediment layer observed upstream of the WRG grounding line (Callens et al., 2014).
After model initialisation, the ice-sheet geometry was relaxed for 50 years prior to the transient simulations, in order to 5 decrease the ice-flux divergence due to artefacts of interpolation and other sources of geometry errors (such as in Cornford et al., 2015). During the relaxation, we used mass conservation to compute melt rates beneath ice shelves (assuming steady state), which gives values in line with current observations (Depoorter et al., 2013;Rignot et al., 2013). However, applying such melt rates beneath the HB ice shelf leads to a rapid retreat of the grounding line during the time span of the relaxation.
We solved this issue by applying a positive basal mass balance (i.e., accretion) of 1 m a −1 during the relaxation, which helps

Transient scenarios
Each initialisation is followed by 12 different transient simulations, applying either linear or nonlinear sliding together with 6 different prescribed sub-ice shelf melt rates, M b1 and M b2 , each with 3 different amplitudes -low, medium and high -15 set by tuning the parameters α 1 and F melt (see Table 1). The naming convention adopted for transient simulations and the corresponding parameters are given in Table 2.
The sum of medium melt rates over the ice shelves yields values that are comparable to current values Depoorter et al., 2013, andM. Depoorter, personal communication, 2016). The sum of low and high melt rates represent approximately half and twice the sum of medium melt rates, respectively. Initial melt rates M b1 and M b2 of medium amplitude   3 Results

Data assimilation
The L-curve analysis performed by Berger et al. (2016) to optimise regularisation still holds for our extended domain and nonlinear sliding, even though it was originally applied to a smaller domain and linear sliding.
The root mean square error between modelled and observed velocities after data assimilation is ≈ 14 m a −1 for B e /S and 5 B e /U initialisations, and ≈ 13 m a −1 for RF/S initialisation, and is independent of the applied sliding law. Such mismatches are similar to what was already computed by control methods applied to the Antarctic ice sheet (e.g. Fürst et al., 2015;Cornford et al., 2015). The largest mismatch is found at the calving front and at the ice rises and promontories. We also find a large mismatch upstream of the TB/HB promontory (Figure 4). We attribute it to the poor consistency between the high observed 9 The Cryosphere Discuss., doi: 10.5194/tc-2016-144, 2016 Manuscript under review for journal The Cryosphere  surface slope and thickness combined with low surface velocities (Supplementary Figure 1), as high driving stresses should induce high velocities. The control method cannot deal with such a non-physical combination for a steady-state ice sheet: it decreases the friction during the first iterations, and further attempts to catch up with the consequent mismatch through ice stiffening during the following iterations.
A significant difference between the two datasets appears in the vicinity of PPw (Figure 4), where the mismatch is lower 5 when using the high-resolution dataset. There, omitting PPw in the control method leads to an excessive ice stiffening ( Figure 5 in Berger et al., 2016).
The central parts of ice shelves are comparatively more viscous, except within rifting areas, where the viscosity can be few orders of magnitudes smaller. The friction coefficient is comparatively small beneath the ice streams of WRG, HB and TB, and few orders of magnitude higher where ice velocity is small, such as in between ice streams and beneath ice promontories and 10 rises. We show these results in Figure 4 with linear sliding.

Initial speed up after unpinning
Unpinning (for B e /U initialisation) induces an instantaneous acceleration of the WRG ice shelf by up to 300 m a −1 at the former location of PPw. After 50 a, the acceleration has propagated over almost the entire ice shelf up to the grounding line, but unpinning does not affect the nearby ice shelves of HB and East Ragnhild Glacier ( Figure 5). The central flowline of the WRG ice stream migrates westward and relocates at an almost equal distance from the HB/WRG promontory and DIR within a few 5 years. The velocities at the ice-shelf front are ≈ 20% larger than for B e /S initialisation. Overall, the comparatively faster ice shelf induces a less advanced grounding line at the end of simulations (about 10 km). The velocity increase near the HB/WRG promontory leads to thinning of its eastward side, making its saddle area afloat and turning it into an ice rise more rapidly than for B e /S and RF/S initialisations. 10 The grounding line migrates similarly for medium melt-rates experiments with linear sliding (shown in Figure 7) and nonlinear sliding. Here we present the common successive steps of all scenarios regarding grounding-line migration and ice dynamics ( Figure 6 and Supplementary Movie).

Main steps of grounding-line migration
The HB ice shelf/sheet system is by far the most dynamic of the three glaciers. During the first 100 a, its grounding line retreats relatively slowly and the pinning points PPhs (Figure 1) detach from the ice-shelf base. The subsequent unpinning of 15 PPhs is followed by an acceleration of the grounding-line retreat over the deepest part of the bed, along with a speed up of ice increasing from ≈ 20% to 100% in a hundred year or so. During these rapid changes, two sudden jumps (the second being less The B e /U initialisation produces faster retreat of grounding lines than the B e /S initialisation, which produces faster retreat than the RF/S initialisation. In particular, the saddle of the HB/WRG promontory gets afloat the most rapidly. The grounding lines of TB and WRG re-advance over up to tens of km for low-melt scenarios.

4 Discussion
Most of the continental shelf beneath the WAIS is deeply depressed, making the ice sheet prone to widespread MISI (Ritz et al., 2015). With respect to the shelf depression, the EAIS is potentially more stable, but its volume of ice is ten times larger than its western counterpart. It is therefore crucial to investigate a potential unstable retreat of grounding lines that may further affect the ice-sheet stability. Here, our simulations systematically show an unstable retreat of HB over the next few hundreds years 10 regardless of the applied sub-ice shelf melt rates, sliding laws and initialisations (Figure 7 and Supplementary Movie). Half of the simulations also predict the retreat of the neighbouring glacier TB for melt rates comparable to current observations.
In total, the contribution of the studied area to sea level rise is 25 ± 10 mm for the next millennium, which needs to be put in perspective with the comparatively small domain (representing about 1% of the Antarctic ice sheet) and the possible nonlinear effects due to future oceanic forcing that are neglected in this study. 15 After a few hundred years, the HB grounding line is quickly retreating at 1 km a −1 , and the ice-shelf velocities reach 600 m a −1 when the grounding line retreats over the most depressed part of the bed (Figure 6). The retreat is only slightly modulated by the type and amplitude of melt rates, indicating that it is mostly driven by a MISI. However, none of the simula-12 The Cryosphere Discuss., doi: 10.5194/tc-2016-144, 2016 Manuscript under review for journal The Cryosphere  Figure 1). This valley is also narrow and starts tens of kilometres upstream of the current grounding line, while the depression beneath the HB grounded ice is wider and starts closer to the grounding line. This accords with the ideal simulations of Gudmundsson et al. (2012), who showed that a wider trough upstream of a grounding line reduces the buttressing exerted by the ice shelf, which enhances the grounding line retreat rate.

5
During the unstable retreat of HB, the ice-shelf thickness is halved compared to the initial thickness. Meanwhile, the thickness of the WRG ice shelf remains almost constant in time near the east side of the HB/WRG promontory. The consequence is an increase of the ice flux coming from the promontory's saddle and going towards the HB ice shelf, reducing the width of the saddle from its western side and eventually making the HB/WRG promontory an ice rise when its saddle becomes afloat. The retreat of TB depends on the melt-rates type and amplitude. All the low amplitude and the M b2 medium amplitude melt rates 10 lead to an advance of its grounding line, while the other scenarios lead to a retreat. However, this contrasting behaviour only slightly modulates the time span by which the saddle of the TB/HB promontory gets afloat, for which the substantial thinning of the HB ice shelf is the major driver.
Current grounding lines fringed and buttressed by ice promontories (such as HB) are relatively stable in the studied area, even resting on upsloping bed (also shown by Gudmundsson et al., 2012, for synthetic numerical experiments). However, 15 13 The Cryosphere Discuss., doi: 10.5194/tc-2016-144, 2016 Manuscript under review for journal The Cryosphere Published: 16 June 2016 c Author(s) 2016. CC-BY 3.0 License.
small amounts of sub-ice shelf melting clearly induce rapid grounding-line retreat and collapse of the promontories into ice rises. This unstable behaviour is corroborated by Favier and Pattyn (2015), showing that promontories are transient features of grounding-line retreat, when they are characterized by an overdeepening upstream of the pinning area.
Most low and several medium melt rates scenarios lead to an advance of the WRG grounding line upstream of DIR (Figure 7), even though we excavated the area below the ice shelf. Because the bathymetry of ice-shelf cavities is poorly constrained, 5 advancing grounding lines must be cautiously interpreted. However, the related effect on sea level calls for a better knowledge of bathymetry beneath ice shelves. Unpinning of the WRG ice shelf mildly affects the global contribution to sea level, which is rather similar to the experiments using B e /S initialisations (Figure 8). However, the decrease of buttressing stemming from unpinning thins the WRG ice shelf and accelerates the retreat of the HB/WRG promontory's saddle from its eastern side: the saddle gets afloat a few hundred years 10 earlier. Such a large difference in timing compared to the differences in sea-level contribution indicates a large sensitivity of promontories deglaciation to a loss of buttressing, similarly to the unstable retreat pointed out in Favier and Pattyn (2015). The Besides the MISI-driven consequences on sea level, sub-ice shelf melting is the other main driver of the retreat. Different behaviours emerge from the two types of melt-rate parametrisations. During the first few hundreds of years, sea-level contri-5 bution is more or less a linear function of melt-rates amplitude. The form of M b1 induces high melt rates at the grounding line when it retreats over the deep trough beneath HB. The contribution to sea level is then a function of pure melting and dynamic thinning, inducing peaks of sea level contribution after about 150 a. In the case of M b2 melt rates, this peak is replaced by a milder bump in sea level contribution (Figure 7) since the pure-melting contribution is lower. After 500 a, the retreat of the HB grounding line is less rapid and the contribution to sea level is then mostly due to melting, and to a lesser extent due to dynamic 10 thinning. Since the M b1 melt rates induce more melting at large depth and almost no melting closer to the surface compared to the M b2 of similar amplitudes, the M b1 melt rates become lower compared to the M b2 melt rates.
Compared to linear sliding, nonlinear sliding (with m = 1/3) enhances basal sliding when ice velocity increases. The acceleration of HB during its unstable retreat consequently yields higher velocities and faster retreat rates of the grounding line for the nonlinear case, hence leading to a higher contribution to sea level (Figure 8). over, any further local change in the boundary condition between the pinning points and the ice shelf, including the extremebut possible -event of unpinning (for instance induced by a substantial thinning of ice shelves; see Paolo et al., 2015) cannot be simulated by the model if the pinning point is omitted in the first place.
Since the early 2000s, uncertainties of ice-sheet modelling outputs have been reduced by substantial numerical improvements, enabling to grasp more accurately key processes such as grounding-line migration (Pattyn and Durand, 2013). This 25 improvement was also made possible by the increasing computational power. We are now able to simulate the behaviour of the WAIS using higher order models at a high spatial resolution in the relevant areas for a wide range of scenarios over the next centuries (Cornford et al., 2015), which was not feasible a few years ago. Nevertheless, the lack of knowledge of essential parameters still affects simulations of the Antarctic ice sheet behaviour, hence preventing further decrease of uncertainties in sea level predictions. Sub-ice shelf melting is a major driver of ice-sheet retreat and sea level contribution (Figure 8). Even though 30 forcing the ice sheet with parametrised melt rates (such as in this study) gives qualitative and informative insights on future sea level contribution, the lack of knowledge of the cavity beneath ice shelves prevents the use of more advanced assessment based on ocean modelling (such as in Hattermann et al., 2014). Moreover, the ill-constrained shape of the ice-shelf cavity dictates how and if the grounding line advances, which also biases future sea level predictions. Here, we demonstrates that sea level predictions and timing of deglaciation can be substantially affected by the type of sliding law, a too shallow bathymetry and the 35

15
The Cryosphere Discuss., doi: 10.5194/tc-2016-144, 2016 Manuscript under review for journal The Cryosphere Published: 16 June 2016 c Author(s) 2016. CC-BY 3.0 License. absence of small pinning points, which all affect ice-sheet initialisation. Also, the exact representation of pinning points (ice rumples, rises and promontories) in the observational datasets, even if they are small, is key for more accurate predictions of future sea level and timing of ice-sheet retreat. Therefore, improving these predictions by the use of ice-sheet modelling relies on future improvements of our knowledge of the bathymetry beneath ice shelves and (small) pinning points.

5
We use the ice-sheet model BISICLES to evaluate the contribution of the outlet glaciers between the Lazarev and Roi Baudouin ice shelves in East Antarctica to future sea level rise, with two different input datasets including or excluding an existing small pinning point (PPw) at the calving front. We also investigate the influence of various sub ice-shelf melt rates parametrisation and two types of weertman-like sliding law (linear and nonlinear). Our results show the likely future unstable retreat of the outlet glacier Hansenbreen (HB) within the next 150 a, which is driven by the marine ice sheet instability (MISI), while the 10 other outlet glaciers are relatively stable over the next millennium. Where the ice sheet is stable (no MISI), sub-ice shelf melting strongly controls sea level contribution. Nonlinear sliding increases the sea level contribution by 20% but does not affect the timing of deglaciation compared to linear sliding. Surprisingly, unpinning (removing PPw after ice-sheet initialisation) hardly impacts the sea level. However, it affects the timing of ice-sheet retreat in the most sensitive parts, such as the HB/WRG promontory which collapses into an ice rise 200 a in advance. On the other hand, omitting PPw during the initialisation of the