Dynamic response of Antarctic Peninsula Ice Sheet to potential collapse of Larsen C and George VI ice shelves

. Ice shelf break-up and disintegration events over the past 5 decades have led to speed-up, thinning, and retreat of upstream tributary glaciers and increases to rates of global sea-level rise. The southward progression of these episodes indicates a climatic cause and in turn suggests that the larger Larsen C and George VI ice shelves may undergo a similar collapse in the future. However, the extent to which removal of the Larsen C and George VI ice shelves will affect upstream tributary glaciers and add to global sea levels is un-known. Here we apply numerical ice-sheet models of varying complexity to show that the centennial sea-level commitment of Larsen C embayment glaciers following immediate shelf collapse is low ( < 2 . 5mm to 2100, < 4 . 2mm to 2300). Despite its large size, Larsen C does not provide strong buttressing forces to upstream basins and its collapse does not result in large additional discharge from its tributary glaciers in any of our model scenarios. In contrast, the response of inland glaciers to a collapse of the George VI Ice Shelf may add up to 8 mm to global sea levels by 2100 and 22 mm by 2300 due in part to the mechanism of marine ice sheet


Introduction
The observational history of ice-shelf collapse in the Antarctic Peninsula has led to a proposed northerly limit of ice-shelf viability determined by the −9 • C mean annual isotherm (Mercer, 1978;Morris and Vaughan, 2003). Recent, rapid warming has led to the southward migration of this limit , now threatening the stability of the large Larsen C and George VI ice shelves. The northernmost remaining ice shelf (Fig. 1a), Larsen C, is considered to present the greatest risk of collapse (Jansen et al., 2015). While other mechanisms such as ice-shelf thinning, fracturing, and weakening of shear margins may contribute to Larsen C ice-shelf instability (Kulessa et al., 2014;Borstad et al., 2016), the risk of shelf collapse has increased slightly since summer 2017 when a large iceberg calved off Larsen C. This calving event leaves Larsen C in conditions similar to those present immediately prior to the collapse of Larsen B Ice Shelf in 2002 and may promote instability (Jansen et al., 2015).
Despite the increased research focus on Larsen C Ice Shelf, most of the current mass loss and contribution to sealevel rise from the Antarctic Peninsula originates from large drainage basins feeding the George VI ice shelf, along the English coast, western Palmer Land, in the south-west of the peninsula (McMillan et al., 2014;Martín-Español et al., 2016). Here, outlet glaciers have thinned rapidly in the last 2 decades, contributing ∼ 0.1 mm a −1 to global sea-level rise (Wouters et al., 2015;Hogg et al., 2017). Many of these glaciers are grounded below sea-level with deeply incised bedrock troughs and retrograde-sloping bedrock topography (Fig. 1b). These marine-based sectors, which contain a sealevel equivalent of 46.2 mm (25 % of the total ice volume in the Antarctic Peninsula Ice Sheet, Fig. 1b), are therefore potentially vulnerable to the marine ice sheet instability mecha-  (Fretwell et al., 2013). The colour bar is truncated at 0 m. Red inset rectangles delineate the locations of zoomed-in areas in Figs. 2, 5, 6, and 8. Black polygons denote ice-sheet model domains.
nism, a tendency of grounding-line retreat to accelerate in the absence of compensating forces (Schoof, 2007;Gudmundsson et al., 2012).
Here we use three ice-sheet models of varying complexity to compute the upstream glacier response and sea-level rise commitment following the potential collapse of the Larsen C and George VI ice shelves. Owing to differences in model set-up and physics, this study does not provide a full model intercomparison, but rather presents a multi-model spread sea-level envelope assessment using a range of ice-flow approximations: (i) the linearised shallow-ice approximation (SIA) model BAS-APISM ; (ii) the hybrid sheet-shelf model PSU3D (Pollard and DeConto, 2012a), and (iii) the vertically integrated sheet-shelf model BISICLES . This multi-model approach provides a starting point for regional ice-sheet model forecasts and sea-level impact studies and allows examination of process differences in glacier responses across the drainage basins of the Larsen C and George VI ice shelves.

Methods
The ice-sheet models BAS-APISM , BISICLES , and PSU3D (Pollard and DeConto, 2012a) have been described in detail elsewhere. A summary of model description, parameterisation, and experimental design relevant to this study are presented here, including important changes to model set-ups from previously published configurations. Model domains vary across the models with BAS-APISM and BISICLES including the entire Antarctica Peninsula and PSU3D simulating the Larsen C embayment and George VI embayment separately (red rectangles in Fig. 1).

Ice-sheet model description
BAS-APISM  simulates ice flow by solving the simplest permissible force balance approximation -the linearised SIA. Owing to the linearisation, the model is less sensitive to ice thickness errors than traditional SIA-based models. The linear nature of the model equations permits a simple summation of sea-level rise contributions from individual drainage basins to provide an ice-sheet-wide estimate. As the SIA is not valid for floating ice shelves (Hutter, 1983), only the grounded ice sheet is simulated and grounding-line retreat is parameterised through a statistical model. In response to ice-shelf collapse, this model scales the expected retreat of the grounding line to the amount of buttressing at the ice front of each drainage basin (Schannwell et al., 2016). Ice-shelf buttressing was computed from the output of an ice-sheet model inversion (Arthern et al., 2015). As BAS-APISM cannot simulate grounding-line advance, ice-shelf flow, or ice-shelf buttressing, this model is only employed in Experiment 1 (immediate ice-shelf collapse) where ice-shelf flow is not explicitly simulated (See Sect. 2.5) and immediate ice-shelf collapse is assumed.
PSU3D (Pollard and DeConto, 2012a) simulates ice flow by using a hybrid combination of the scaled SIA and shallowshelf approximation (SSA) equations. The SSA is valid for ice shelves and ice streams characterised by low basal drag. This type of ice-sheet model (A-HySSA, asymptotic hybrid SIA-SSA model;  provides the required physics to simulate the ice sheet -ice shelf system, including explicit tracking of the position of the grounding line. To make the model less sensitive to grid resolution, an additional internal flux boundary condition is employed at the grounding line. The model set-up used here is similar to Pollard et al. (2015), but cliff failure and bedrock deformation are not included. PSU3D solves the time-varying 3-D temperature equation, but surface air temperature forcing is held constant at the year 2000 (Le Brocq et al., 2010) throughout the simulations.
BISICLES  simulates ice flow by solving a vertically integrated stress balance (L1L2 = onelayer longitudinal stress model; Hindmarsh, 2004) to determine the horizontal velocity. The ice rheology is given by Glen's flow law Here S is the deviatoric stress tensor, η is the effective viscosity,˙ is the strain-rate tensor, and φ is the stiffening factor that accounts for ice damage, anisotropy, and temperature uncertainties . This type of stress balance is similar to the SSA, but includes vertical shearing in the effective viscosity calculation, resulting in softer ice at the grounding line in comparison to traditional SSA models and more resembles the behaviour of full Stokes models . The equations are solved on an adaptive 2-D grid, allowing higher resolutions in areas of interest such as grounding lines or shear margins and coarser resolution away from these regions to save computational time. A subgrid interpolation scheme for basal drag near the grounding line was employed to improve the accuracy of the grounding-line position at each time step . In all BISICLES simulations ice temperature data are provided by a three-dimensional thermo-mechanical model (Pattyn, 2010) and is held fixed in time.
Basal traction in PSU3D and BISICLES is determined by a viscous law where m = 0.5 (quadratic law), τ b is the basal traction, u is the horizontal velocity, ρ i and ρ w are ice and ocean densities, b is the bedrock elevation, h is ice thickness, and C is the basal friction parameter inferred by solving an inverse problem (See Sect. 2.3). Due to the linearisation of the evolution equations in BAS-APISM, there is no need to specify whether or not basal sliding is occurring. All rates are determined by the ice flux, which is directly derived from the data . Basal sliding sensitivity simulations with BISICLES were also performed with m = 1/3 (cubic law) and m = 1 (lin-ear law). In addition, a simulation was performed using a Coulomb-limited law (Tsai et al., 2015). This law combines the power law (Eq. 2) with the Coulomb friction law by ensuring that basal traction cannot exceed the Coulomb friction that is proportional to the effective pressure N e : where the first term in the parentheses is the Coulomb friction law with a = 0.5, m = 0.5 and the effective pressure N e is where g is the acceleration of gravity and h f is the flotation thickness. Equation (4) is only valid under the assumption of a full connection between the basal hydrology and the ocean.
Since the Coulomb law implies that basal drag approaches zero towards the grounding line, this type of basal sliding law ensures a smooth transition from grounded to floating ice, unlike the traditional power law (Eq. 2) which implies that basal drag is highest near the grounding line (Tsai et al., 2015).

Calving
In simulations where the calving front is not fixed, e.g. where ice-shelf flow and retreat are explicitly simulated, calving depends on the depths of surface (d s ) and basal crevasses (d b ), relative to total ice thickness. Crevasse depths are computed by Benn et al. (2007) and Nick et al. (2010).
where A is the depth-averaged rheological coefficient, n = 3 is the rheological exponent, d w is the water height in the surface crevasse, and ρ 0 is the density of surface water. The parameter˙ is the longitudinal strain rate approximated in PSU3D through the isotropic ice divergence: BISICLES essentially implements the same criterion, but computes crevasse depths from membrane stresses such that Eqs. (5)-(6) become where τ is the deviatoric stress tensor, tr() is the trace operator, and h ab is the thickness above flotation (Sun et al., 2017).
In PSU3D ice is calved off when the combined ice thickness of the surface and bottom crevasses reach at least 75 % of the column ice thickness (Pollard et al., 2015), whereas in BISICLES icebergs calve when the sum of the surface and bottom crevasses reaches the distance from the ice surface to the waterline. Water height in surface crevasses (d w in Eq. 5) is computed from biased-corrected CMIP5 model projections up to 2300 from the model selection presented in Schannwell et al. (2015). The bias correction and melt computation approach follow Trusel et al. (2015). In brief, December-January-February (DJF) near-surface temperatures from the CMIP5 historical simulations were compared to high-resolution (5.5 km) RACMO2.3 simulations (van Wessem et al., 2016) such that where T 2m GCM and T 2m RACMO2.3 are the mean DJF near-surface temperatures over the baseline period 1980-2005 from GCM and RACMO2.3, respectively. The bias calculation (Eq. 10) is restricted to the ice-shelf areas in our two model domains ( Fig. 1). The best-performing GCM (lowest bias) for the RCP4.5 (Fig. A1, MIROC-ESM) and RCP8.5 (Fig. A2, CSIRO) scenarios were then selected as future forcing.
To convert near-surface temperature to melt, the empirical formula derived by Trusel et al. (2015) was used. This formula scales surface melt (R in Eq. 11) exponentially with mean DJF near-surface temperatures and approximates the surface melt available to fill surface crevasses. To compute water height in surface crevasses, d w is set to (Pollard et al., 2015) d w = 100R 2 . (11)

Model initialisation
BAS-APISM employs a combined altimetric and velocity initialisation scheme, permitting a steady-state starting condition after initialisation under the assumption that the current ice sheet configuration is close to steady state . This is accomplished through the computation of balance fluxes. The motivation for this type of initialisation technique is that the absence of accurate ice thickness data sets leads to the omission of the mechanical model in the cost function employed for the initialisation . BISICLES is initialised by solving an optimisation problem to infer the basal traction coefficient C and the stiffening factor φ (also enhancement factor, Eq. 1), by matching modelled velocities with observed velocities (Rignot et al., 2011). This type of initialisation is well known and widely employed in ice sheet modelling (MacAyeal, 1992;Cornford et al., 2015). A non-linear conjugate gradient method was employed to seek a minimum of the objective function where J m is the misfit between observed and modelled velocities and J p is a Tikhonov penalty function described by where λ C and λ φ are the Tikhonov parameters and J reg C and J reg φ represent the spatial gradients of C and φ integrated over the domain . An L-curve analysis was performed to calibrate the Tikhonov parameters and avoid overfitting or over-regularisation (Fürst et al., 2015). The selected values are λ C = 10 −1 and λ φ = 10 9 (Fig. A3).
To solve this inverse problem, maps of surface elevation and bedrock topography were taken from the BEDMAP2 (Fretwell et al., 2013) data set, and a steady state 3-D temperature field was used from a higher-order model (Pattyn, 2010). It is only necessary to find solutions with a single sliding law, as coefficients can be computed from one another to give the same basal traction τ b ; e.g. the coefficients for the cases m = m 1 and m = m 2 must satisfy C 2 |u| m 2 = C 1 |u| m 1 . We chose m = 1 for the inversion simulation.
PSU3D utilises a different algorithm to infer the basal traction coefficient. Instead of matching velocities, the algorithm implemented in PSU3D seeks to minimise the misfit between local surface elevation observations and modelled local surface elevations (Pollard and DeConto, 2012b). To achieve this, the ice-sheet model is run forward in time, and basal traction coefficients are periodically compared and adjusted according to the local surface elevation error. This iterative process is continued until the modelled surface elevation converges to the best fit with observed surface elevation (Pollard and DeConto, 2012b). Note that this simpler algorithm does not infer a stiffening factor φ for ice shelves. Input maps needed for the inversion algorithm are from ALBMAP (Le Brocq et al., 2010), e.g. on ice thickness and bedrock topography. For all PSU3D simulations, basal traction fields are interpolated onto the respective model grid from a 5 km Antarctica inversion simulation. The coarser resolution leads to some interpolation artefacts in the basal traction coefficient fields (Fig. 2).

Spin-up
Following initialisation, the sheet-shelf models should aim to be as close to the steady-state initial conditions provided by observations, as long as the ice sheet itself is in steady state, such that ∂h ∂t = 0. However, owing to data inconsistencies and in part to a violation of this steady-state assumption, this condition is not fulfilled, requiring a spin-up or relaxation simulation to reach a steady state for each model. To tease out the sea-level rise contributions from ice-shelf removal and facilitate comparison across all three ice-sheet models, the employed spin-up approach aims to keep the ice sheet geometry as close as possible to the initial geometry. This is necessary because BAS-APISM provides a stable starting condition after initialisation. To ensure a minimal change in ice-sheet geometry, we compute a synthetic mass balance (MB) which is where FC is the negative of the modelled thickness field change when the model is run forward a single time step. This synthetic mass balance is applied in all spin-up and perturbation simulations. All simulations are then run forward in time for 50 years with only this forcing applied. To reach steady state, the volume above flotation change with time should be near zero ( ∂V ∂t ∼ 0) at the end of the spin-up. All of our simulations fulfil this criterion (Fig. 3), even though PSU3D simulations are not as close to steady state as BISI-CLES simulations at the end of the spin-up period.

Experimental design
Two sets of experiments were undertaken using the icesheet models (Table 1). In Experiment 1, immediate iceshelf collapse was imposed on all three ice-sheet models and combined with a fixed calving front position. This provides an envelope of sea-level rise projections for the peninsula region and allows for an evaluation of the importance of each shelf to the tributary glaciers upstream. Simulations with the sheet-shelf models (PSU3D and BISICLES) were carried out at different horizontal resolutions to investigate the grid dependence on the sea-level rise projections and to select the best compromise between computational demand and appropriate grid resolution (Table 1). In the second simulation (Experiment 2), the two sheet-shelf models (PSU3D and BISICLES) were run at 1 km resolution to simulate ice-shelf retreat and collapse and subsequent tidewater glacier retreat using a physically based calving relation . Volume change ∂V ∂t spin-up plot for BISICLES (solid lines) and PSU3D (dashed lines) at different horizontal resolutions. (Benn et al., 2007;Nick et al., 2010). This relation initiates iceberg calving when the combined depth of surface and bottom crevasses reaches a threshold percentage of ice thickness (See Sect. 2.2). Crevasse depth primarily depends on the stress field of the ice shelf, with extensional stresses providing favourable conditions for crevasse opening. Meltwater hydrofracture may also increase calving rates in a process that has been strongly implicated in the 2002 collapse of the Larsen B Ice Shelf (Scambos et al., 2003). In all simulations of Experiment 2 ice-shelf thickness is allowed to evolve freely. This more realistic experiment permits the evaluation of a more gradual loss of buttressing to the upstream glaciers and assesses the effect of a dynamic calving front. In all simulations, perturbations to the surface mass balance are ignored as these are expected to be small in comparison to ice dynamic changes resulting from shelf loss . Moreover, ocean melting is set to zero in the perturbation experiments unless stated otherwise.

Experiment 1: immediate ice-shelf collapse
Projections of sea-level rise from Larsen C embayment glaciers following immediate shelf collapse (Experiment 1) are small, ranging from 0.5 to 1.5 mm by 2100 and from 0.6 to 1.6 mm by 2300 (Fig. 4a). The sea-level curve rises in the first 2 decades in response to a loss of backstress provided by the shelf, then decelerates with tributary glaciers, adjusting to the new configuration ∼ 25 years after collapse. Grounding-line retreat of > 5 km and extensive dynamic thinning (> 0.6 m a −1 , propagating ∼ 75 km inland) are restricted to five outlet glaciers in the southern part of the embayment (Fig. 5). In contrast, immediate collapse of George VI Ice Shelf perturbs upstream grounded tributaries by up to 0.8 m a −1 averaged over 300 years and results in a total sea-level rise of 4-11 mm by 2300 (Fig. 4b). The more dramatic response from George VI tributary glaciers means that they have not yet reached steady state by 2100 (Fig. 4d), and PSU3D simulations continue to contribute to sea level well beyond this date. This discrepancy between the sheet-shelf models may be attributed to a combination of differences in initialisation, inferred basal traction fields, and PSU3D not being as close to steady state as BISICLES following initialisation and spin-up (Figs. 3 and A4). Moreover, icesheet thinning in response to the collapse event propagates further upstream in PSU3D and is more widespread than in BISICLES, leading to higher rates of mass loss despite similarly predicted grounding-line retreat (Tables A1 and A2). This response has been previously attributed to differences in the underlying model physics (L1L2, A-HySSA). Using synthetic geometries, A-HySSA models have been shown to be more sensitive to grounding-line advance as well as retreat. These differences are most likely caused by neglected vertical shearing terms in the pure membrane ice-sheet models .
While there is a notable grid dependence in BISICLES projections (Fig. A6), this is much reduced in the PSU3D projections, supporting the findings of previous modelling studies (Pollard and DeConto, 2012a) that ice-sheet models with the implementation of an internal flux boundary condition are less sensitive to grid resolution. The required firstorder convergence  of the sea-level rise projections in the BISICLES simulations is met for simulations at 1 and 0.5 km resolution. To facilitate comparison between the two sheet-shelf models, a 1 km grid is employed for Experiment 2.
At 1 km resolution, the combined sea-level rise by 2300 from glaciers in both embayments ranges from 5.6 to 10.4 mm sea-level equivalent, with > 60 % of the total provided by George VI outlet glaciers. BAS-APISM projects a similar total to 2300 (11.5 mm), though a poor match in simulated spatial patterning of dynamic thinning (Figs. 5 and 6) shows that the simplified model physics and statistical approach to grounding-line retreat do not perform satisfactorily in some areas.
Across all three ice-sheet models, and in both Larsen C and George VI embayment domains, ice-shelf collapse does not result in widespread and extensive grounding-line retreat (Figs. 5 and 6). This was expected for Larsen C outlet glaciers due to a combination of prograde-sloping bedrock topography and the moderate backstress currently provided by the shelf (Fürst et al., 2016). George VI Ice Shelf, however, provides both strong buttressing (Fürst et al., 2016) and mostly marine-based outlet glaciers on retrograde sloping bedrock topography (Fig. 1b), conditions expected to be favourable for marine ice sheet instability. Despite this, grounding-line retreat of George VI outlet glaciers is limited to a few locations and < 15 km in length (Fig. 6). These findings suggest that stabilising forces such as basal and lateral drag may provide enough resistance for the ice sheet in western Palmer Land to remain in a stable configuration following the initial response to ice-shelf collapse. This is supported by earlier modelling studies with idealised geometries, showing that the magnitude of grounding-line retreat is a function of the retrograde sloping channel width (Gudmundsson et al., 2012;Gudmundsson, 2013). The smaller the channel width, the less retreat was simulated (Gudmundsson et al., 2012). Considering the small size of the drainage basins in the peninsula region with channel widths < 30 km, the remaining lateral buttressing from shear margins likely impedes any runaway grounding-line retreat.

Experiment 2: gradual ice-shelf retreat
When ice-shelf frontal changes are explicitly simulated (Experiment 2) with sheet-shelf models (PSU3D, BISICLES) using a stress-field-dependent calving law, sea-level rise projections span a much larger range. With forcing from the representative concentration pathway (RCP) 8.5 high-emission scenario, Larsen C and George VI embayment basins com-  bined provide up to 23 mm sea-level equivalent ice loss by 2300 (Fig. 7b), with 95 % of this total coming from George VI tributary glaciers. The contribution to the sealevel budget from Larsen C embayment glaciers is small (< 1.5 mm) and remains equivalent to Experiment 1. The sea-level commitment from Larsen C glaciers is modest as complete shelf collapse is not forecast until 2150 in RCP8.5 (Fig. 7c), and only 45-60 % of the shelf area is lost by 2300 in RCP4.5 ("business-as-usual" scenario). This leads to limited grounding-line retreat, and dynamic thinning is restricted to five outlet glaciers in the southern part of the embayment (Fig. 8a). The larger grounded area loss simulated with PSU3D for Larsen C (Fig. 7c) is not a response to the loss of buttressing force, rather it is due to a more seaward advanced initial grounding-line position introduced in the model spin-up phase, the effect of which is small on sealevel projections (0.28 mm sea-level equivalent).
Although projections for Larsen C Ice Shelf glaciers agree reasonably well in absolute numbers across both sheet-shelf models despite differences in their underlying physics, projections diverge for simulations of George VI Ice Shelf glaciers (Fig. 7b). Experiment 2 (dynamic calving front) sim-  ulations with PSU3D provide very similar sea-level projections to Experiment 1 runs for George VI (immediate iceshelf collapse: 6.8-7.1 mm from RCP4.5 and 8.5, respectively). In contrast, BISICLES projects little sea-level rise under RCP4.5 for George VI, as the amount of available meltwater for hydrofracturing is insufficient to initiate iceshelf collapse or retreat due to the different implementations of the calving law (See Sect. 2.2). Under RCP8.5, break-up of the George VI Ice Shelf occurs approximately in 2100 (Fig. 7d), resulting in widespread grounding-line retreat and sea-level rise of 22 mm by 2300. Both sheet-shelf models project similar sea-level rise up to 2150 in Experiment 2 for the George VI domain. In the BISICLES RCP8.5 simulation, however, following the collapse of the shelf, calving fronts and grounding lines retreat further back into the marine-based sectors (Fig. 8d). After ice-shelf collapse, the grounding line and calving front for all drainage basins are almost in identical locations. Increasing rates of calving permit the grounding line to retreat much further inland in the RCP8.5 BISICLES simulation for George VI. As this enhanced grounding-line retreat is only present in Experiment 2, it suggests that this retreat is most likely due to a combi-The Cryosphere, 12, 2307-2326, 2018 www.the-cryosphere.net/12/2307/2018/ nation of the dynamic calving front and the marine ice-sheet instability mechanism. Even with a dynamic calving front, enhanced grounding-line retreat for George VI is not triggered before some time after ice-shelf collapse (> 15 years, Tables A3 and A4), indicating that fast grounding-line retreat is not triggered until the calving front and the grounding line reach a retrograde sloping bedrock topography. As a result of widespread grounding-line retreat for George VI in the RCP8.5 scenario (Fig. A5), extensive dynamic thinning occurs (> 1 m a −1 ), extending up to 100 km inland in the southern parts of the embayment (Fig. 8d). PSU3D simulations do not show enhanced grounding-line retreat in this sector. The discrepancy in sea-level rise projections between Experiment 1 and Experiment 2 is a result of the different applied perturbations. In Experiment 1, the entire ice shelf is removed at the start of the simulation before a fixed calving front is employed. In contrast, Experiment 2 with the crevasse calving law has much more potential to vary. Our simulations show that either very little ice can calve (Fig. 7, RCP4.5 scenario) or given enough surface water the entire shelf can collapse and emerging new floating areas that were formerly grounded keep on calving. So, unlike Experiment 1 where collapse is only enforced once, repeated/continuing collapse of the shelf can occur in Experiment 2 (Fig. 7, RCP8.5 BISICLES simulation).
We attribute the good agreement across both models for Larsen C to the fact that the area of the marine-based sectors is limited in this domain (2.1 mm contained in marinebased sectors) due to the very mountainous bedrock topography constraining potential grounding-line retreat. This is supported by all simulations across all ice-sheet models as even under a wide range of different forcings the Larsen C embayment does not contribute more than 4.2 mm by 2300. The greater potential to initiate grounding-line retreat is presented by the George VI Ice Shelf, where much of the ice sheet is marine based with retrograde sloping bedrock topography (Fig. 1b). As this large grounding-line retreat is only initiated in the BISICLES simulation, large differences in sea-level rise projections occur. The most likely explanation for this differing behaviour is due to the difference in the inferred basal traction coefficient fields that affect each model's response to ice-shelf removal. PSU3D predicts much higher-friction bedrock conditions in the George VI embayment than BISICLES (Fig. 2). These high-friction bedrock conditions result in little acceleration of the major outlet glaciers following ice-shelf break-up. This in turn means that the calving law applied only to floating ice cells cannot drive the initial retreat into the marine-based sectors as the outlet glaciers do not thin sufficiently to form floating ice tongues. In contrast in the RCP8.5 BISICLES simulation for George VI, speed-up in response to ice-shelf break-up leads to enhanced dynamic thinning of the main outlet glaciers. This thinning in conjunction with the calving law drives the calving front into the marine-based sectors where further retreat is initiated by a combination of the marine ice-sheet instabil-ity and the meltwater-driven calving law, resulting in much higher simulated sea-level rise projections.

Uncertainty assessment
A range of sensitivity experiments were undertaken to assess the robustness of our model simulations to additional forcings. To assess the impact of an additional ocean forcing, a pair of basal melt anomalies were applied to areas of fully floating ice in addition to the freely evolving calving front forcing (Experiment 2). In a first, moderate simulation, the anomaly was set to the current thinning signal of the respective ice shelf  for the duration of the forecast period (0.5 m a −1 for Larsen C, 1.1 m a −1 for George VI). In a second, extreme scenario, the same initial anomaly was applied, then increasing linearly to 3 times the current thinning signal by 2100, remaining at this magnitude to 2300. In each case, sea-level projections with these additional forcings are within 0.2 mm sea-level equivalent of simulations without additional forcings: in other words, it is ice-shelf break-up in combination with the calving criteria that dominates our results. As the basal boundary condition remains poorly constrained in ice-sheet models, yet our model projections show a strong dependence on this condition, Experiment 1 (immediate shelf collapse) was repeated with BISICLES using a range of basal sliding laws. Each of the traditionally employed power laws result in similar sealevel rise projections to 2300 (1.4-1.6 mm for Larsen C embayment glaciers and 4-6 mm for George VI glaciers; Fig. 4). Projections to 2300 increase by a factor of 2 for simulations using a Coulomb-limited sliding law (Tsai et al., 2015), resulting in ∼ 3 mm from Larsen C glaciers and ∼ 12 mm for George VI glaciers. This type of basal sliding law reduces the basal drag in a mobile ∼ 1 km layer which forms immediately upstream of the grounding line, resulting in greater discharge throughout the simulations (Tables A1 and A2).
The importance of better-constrained boundary conditions in the peninsula region (bedrock topography and ice thickness) is highlighted by a discrepancy between sea-level rise projections for Larsen C embayment basins using different data products. Although total ice volume and ice volume below sea-level differences between ALBMAP and BEDMAP2 products are small (< 15 %), a more recent higher-resolution data set (Huss and Farinotti, 2014) provides an increase of ∼ 100 % in ice volume below sea level. When incorporated into ice-sheet model simulations, this altered bedrock topography results in larger grounding-line retreat rates for some basins occupying deeper bedrock troughs. A consequence is a sea-level rise projection for Experiment 1 (immediate iceshelf collapse) with the reference sliding law (Eq. 2) that increases by a factor of ∼ 3 (4.2 mm) for the Huss and Farinotti (2014) bedrock topography data set, underlining the significance of an accurate boundary data set for sea-level rise projections. In addition, our experiments show that, for simulations of grounding-line motion in response to ice-shelf break-up, sheet-shelf models are necessary. The simple model BAS-APISM fails to reproduce the results of the sheet-shelf models due to the simplified physics. Even across sheet-shelf models differences in model physics, model initialisation, calving law implementation and other numerics (e.g. meshing) can lead to substantially different projections under the same forcing (Fig. A5). Sea-level rise projections are the most sensitive to sliding law and bedrock geometry. The peninsula is not the only region where these parameters greatly affect decadal to centennial sea-level rise projections as similar conclusions were drawn from models of outlet glaciers in the Amundsen Sea embayment (Nias et al., 2018). The wide range of sea-level rise responses to different forcing parameters underlines the need for perturbed ensembles to explore key parameter uncertainties (e.g. basal sliding law) for sea-level rise projections in greater detail for the peninsula region. Owing to the increase in computational power, these types of ensemble projections have become feasible on regional (e.g. Nias et al., 2016) and continental scales (e.g. DeConto and Pollard, 2016).

Comparison with Larsen B Ice Shelf collapse response
To further assess the impact of ice-shelf break-up, five drainage basins from the Larsen C embayment (LarI-LarV, Fig. 8) and George VI embayment (GeoI-GeoV, Fig. 8) were selected for additional analysis. This provides a comparison to real-world examples of the magnitudes and pattern of glacier response to ice-shelf collapse. For Experiment 1 (immediate ice-shelf collapse), the speed-up following iceshelf removal is short lived (∼ 15 years) for both models. Maximum speed-up of ∼ 300 % is possible, though the mean maximum speed-up is ∼ 50 % (Tables A1 to A4). These values are smaller than those observed following Larsen B collapse with a maximum of 8-fold speed-up (Rignot et al., 2004). This may be due to the different areas selected for the speed-up calculation. Both rates of ice discharge (mass loss) and grounding-line retreat are highest immediately following shelf collapse. For 65 % of the selected 10 drainage basins, more than 50 % of the total modelled grounding-line retreat takes place within 15 years of ice-shelf collapse. Maximum mass loss rates for Larsen C (1.6-5.1 Gt a −1 ) are smaller than observations for a similar time period for Larsen B (8.0 Gt a −1 ) (Scambos et al., 2014). Responses of individual drainage basins for Larsen C and George VI are highly variable with grounding-line retreat ranging from 3.7 to 26.8 km for Larsen C and from 3.1 to 12.2 km for George VI (Tables A1 and A2). This high spatial variability in response across the selected basins indicates that the importance of ice-shelf buttressing is also highly spatially variable. Most of the grounding-line retreat, in particular for Larsen C, occurs in areas of bedrock channels (Fig. 5). Since these deep bedrock channels are smaller for Larsen C, this leads to smaller mass loss than for George VI, even though maximum grounding-line retreat numbers are larger. But groundingline retreat is spread across a wider area of the drainage basin front at George VI (Fig. 6). For Experiment 2, maximum grounding-line retreat (4 to 33.4 km, Tables A3 and A4) is of similar magnitude to Experiment 1 (3.1 to 26.8 km, Tables A1 and A2), including the spatial variability across the selected basins for PSU3D in both embayments. However, for BISICLES this only holds true for the Larsen C embayment where mass loss averaged across the five basins remains at 0.1 Gt a −1 for both experiments. For the George VI domain the spatial variability across the basins is strongly reduced with maximum grounding-line retreat now ranging from 10 to 25.5 km (Table A4). When averaged over the five basins, grounding-line retreat increases from 6.4 km for Experiment 1 to 21.3 km in Experiment 2 for George VI. The retreat in this experiment spreads over the entire width of the drainage basin www.the-cryosphere.net/12/2307/2018/ The Cryosphere, 12, 2307-2326, 2018 front (Fig. 8d), resulting in an increase in mass loss over the 300 years from 0.5 Gt a −1 in Experiment 1 to 3.0 Gt a −1 in Experiment 2. This is in agreement with computed sea-level rise projections (Fig. 7). Significant speed-up is absent in the years following ice-shelf removal across all basins due to the more gradual loss of buttressing in Experiment 2 (compared to the complete ice-shelf removal in Experiment 1). This results in a less dramatic dynamic response than in Experiment 1, with the exception of several basins of the George VI Ice Shelf where retreat rates can lead to large mass losses. The gradual loss of buttressing simulated by Experiment 2 leads to grounding-line retreat and mass loss response occurring > 15 years after ice-shelf removal.

Conclusions
The most important contributor to the global sea-level budget up to 2300 from Antarctic Peninsula ice-shelf -ice sheet dynamics is glaciers in western Palmer Land feeding George VI Ice Shelf. Our envelope of sea-level rise projections ranges from 4-12 mm sea-level equivalent water by 2300 in Experiment 1 to 6-22 mm sea-level equivalent water by 2300 in Experiment 2 for George VI. As the highest projection represents only 55 % of the grounded ice below sea level in this region (Fretwell et al., 2013), there may yet be even more ice at risk to dynamic mass loss. This is in contrast to the Larsen C embayment, where the majority of the projections result in a loss of more than 55 % of grounded ice below sea-level. The highest projections remove all grounded ice below sea level as well as ice elsewhere in the basin. The sea-level rise induced by ice-shelf removal in the Larsen C embayment is therefore largely limited by the small area of marine-based sectors. All projections are relatively insensitive to increased ocean forcing yet are highly sensitive to changes in the basal boundary condition and to the boundary data set, highlighting the need for improved bed topography data and a more rigorous uncertainty analysis. While Larsen C Ice Shelf's recent calving event may increase its vulnerability to ice-shelf instability, our simulations under a wide range of future forcing scenarios show that the sea-level commitment of Larsen C embayment glaciers following shelf collapse or retreat is limited to less than 4.2 mm by 2300 (0.6-4.2 mm for Experiment 1; 0.4-1.5 mm for Experiment 2). Individual drainage basin analysis indicates a wide range of responses to ice-shelf removal, but overall ice-flow speed and mass changes are expected to be of similar magnitude to those observed following the 2002 collapse of the Larsen B Ice Shelf.
Data availability. Data sets are available upon request from the corresponding author.      Table A1. Maximum grounding-line retreat (dGL, km), mass change rate (dM/dt, Gt a −1 ) and speed-up (dU/dt 15 ) for selected sample basins in Larsen C embayment for Experiment 1 (immediate collapse). Subscript indicates that numbers are averaged over 15 years after ice-shelf collapse, with the exception of dGL 15 , which is the sum over this period. For columns without subscript the numbers are for the entire 300-year simulation period. For speed-up, in addition to the average number over 15 years (dU/dt 15 ), the maximum speed-up in this time period is provided (dU/dt max ). Speed up was calculated for a region within 5 km of the current grounding line. As BAS-APISM does not simulate the ice shelf, speed-up calculations were not carried out.