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

Ice shelf break-up and disintegration events over the past several 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 similar collapse in future. However, the extent to which removal of Larsen C and George VI ice shelves will affect upstream tributary glaciers and add to global sea levels is unknown. Here we apply numerical ice-sheet models of varying complexity to show that 5 the centennial sea-level commitment of Larsen C embayment glaciers following immediate shelf collapse is low (<2.5 mm to 2100, <4.3 mm 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 collapse of 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 instability. Our results demonstrate the varying and relative importance 10 to sea level of the large Antarctic Peninsula ice shelves considered to present a risk of collapse.


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 climate (Mercer, 1978;Morris and Vaughan, 2003).Recent, rapid warming has led to the southward migration of this limit (Vaughan et al., 2003), now threatening the stability of the large Larsen C and George VI ice shelves.
The northernmost remaining ice shelf (Figure 1a), Larsen C, is considered to present the greatest risk of collapse (Jansen et al., 2015).While other mechanisms may contribute to Larsen C ice-shelf stability (Kulessa et al., 2014;Holland et al., 2015;Borstad et al., 2016), the risk of shelf collapse has increased 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).(b) Bedrock elevations below sea level (meters above sea level (m a.s.l)) for the Antarctic Peninsula from BEDMAP2 (Fretwell et al., 2013).
Red inset rectangles delineate location of zoom-in views in Figure 8. Black polygons denote ice-sheet model domains.
Despite the increased research focus on Larsen C Ice Shelf, most of the current mass loss and contribution to sea-level rise from the Antarctic Peninsula originates from large drainage basins feeding 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 two 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 (Figure 1b).These marine-based sectors, which contain a sea-level equivalent of 46.2 mm (25% of the total ice volume in the APIS, Figure 1b), are therefore potentially vulnerable to the marine ice sheet instability mechanism, a tendency of grounding-line retreat to accelerate in the absence of compensating forces (Schoof, 2007;Gudmundsson et al., 2012).
Here we use three state-of-the-art ice-sheet models of varying complexity to compute the upstream glacier response and sealevel rise commitment following collapse of Larsen C and George VI ice shelves.Owing to differences in model setup 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 (Barrand et al., 2013); (ii) the hybrid sheet-shelf model PSU3D (Pollard and DeConto, 2012a), and; (iii) the vertically-integrated sheet-shelf model BISICLES (Cornford et al., 2013).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 Larsen C and George VI ice shelves.

Methods
The ice-sheet models BAS-APISM (Barrand et al., 2013), BISICLES (Cornford et al., 2013), and PSU3D (Pollard and De-Conto, 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 setups from previously published configurations.
2.1 Ice-sheet model description BAS-APISM (Barrand et al., 2013) simulates ice flow by solving the simplest permissible force balance approximation -the linearised shallow-ice approximation (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 simple summation of sea-level rise contributions from individual drainage basins to provide an ice-sheet wide estimate.As the SIA is not valid at the grounding line, only the grounded ice sheet is simulated and grounding-line retreat is parameterised through a statistical model.This model scales the expected retreat of the grounding line in response to ice-shelf collapse to the amount of buttressing at the ice front of each drainage basin (Schannwell et al., 2016).Ice-shelf buttressing was computed from 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 Section 2.5) and immediate ice-shelf collapse in assumed.
PSU3D (Pollard and DeConto, 2012a) simulates ice flow by using a hybrid combination of the scaled SIA and shallow-shelf 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 (Pattyn et al., 2013)) 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 at the grounding line, an additional internal flux boundary condition is employed.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 year 2000 (Le Brocq et al., 2010) throughout the simulations.BISICLES (Cornford et al., 2013) simulates ice flow by solving a vertically integrated stress balance (L1L2 = one-layer longitudinal stress model (Hindmarsh, 2004)) to determine the horizontal velocity.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 resembles more the behaviour of full-Stokes models (Pattyn and Durand, 2013).The equations are solved on an adaptive 2-D grid, allowing for higher resolution in areas of interest such as grounding lines or shear margins, and coarser resolution away from these regions to save computation 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 (Cornford et al., 2016).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 Section 2.3).
Basal sliding sensitivity simulations with BISICLES were also performed with m=1/3 (cubic law) and m=1 (linear law).In addition, a simulation was performed using a Coulomb-limited law (Tsai et al., 2015).This law combines the power law (Equation 1) 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 3 is only valid under the assumption of 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 (Equation 1) 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 is 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;Nick et al., 2010) where Ā is the depth averaged rheological exponent, n=3 is the rheological exponent, d w is the water height in the surface crevasse and ρ 0 is the density of surface liquid.
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 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 sea level to the ice surface.
Water height in surface crevasses (d w in Equation 4) is computed from biased-corrected CMIP5 model projections to 2300 from the model selection presented in Schannwell et al. (2015).The bias-correction and melt computation approach follows 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.3simulations ( van Wessem et al., 2016) such that where T 2m GCM and T 2m RACM O2.3 are the mean DJF near-surface temperatures over the baseline period 1980-2005 from each GCM and RACMO2.3,respectively.The bias calculation (Equation 9) is restricted to the ice-shelf areas in our two model domains (Figure 1).The best performing GCM (lowest bias) for the RCP4.5 (Figure A1, MIROC-ESM) and RCP8.5 (Figure A2, CSIRO) scenarios were then selected as future forcing.
To convert from near-surface temperature to melt, the empirical formula derived by Trusel et al. (2015) was used.This approximates the surface melt available to fill surface crevasses (R).To compute water height in surface crevasses, d w is set to (Pollard et al., 2015).

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 (Barrand et al., 2013).
The motivation for this type of initialisation technique is that the absence of accurate ice thickness datasets leads to the omission of the mechanical model in the cost function employed for the initialisation (Barrand et al., 2013).
BISICLES is initialised by solving an optimisation problem to infer the basal traction coefficient C and the stiffening factor φ (also enhancement factor), 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 nonlinear The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-21Manuscript under review for journal The Cryosphere Discussion started: 6 February 2018 c Author(s) 2018.CC BY 4.0 License.
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 (Cornford et al., 2015).An L-curve analysis was performed to calibrate the Tikhonov parameters and avoid overfitting or overregularisation (Fürst et al., 2015).The selected values are λ C = 10 −1 and λ φ = 10 9 (Figure A3).
In solving this inverse problem, maps of surface elevation and bedrock topography were taken from the BEDMAP2 (Fretwell et al., 2013) dataset, 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 the 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 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 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., ice thickness, bedrock topography, and surface mass balance.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 (Figure 2).

Spin-Up
After initialisation the sheet-shelf models should be in equilibrium, providing ∂h ∂t = 0.However, this condition is not fulfilled, requiring a spin-up or relaxation simulation to reach a steady state for each model to guarantee a steady-state starting position.
To 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 simply (Price et al., 2017)  (Figure 3), even though PSU3D simulations are not as close to steady state as BISICLES simulations at the end of the spin-up period.

Experimental Design
Two sets of experiments were undertaken using the ice-sheet models.In the first simulation (hereafter Experiment 1), immediate ice-shelf collapse was imposed on all three ice-sheet models and combined with a fixed calving front position.This provides 5 7 The Cryosphere Discuss., https://doi.org/10.an envelope of sea-level rise projections from the best available ice-sheet models for the peninsula region and evaluates 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.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 (Benn et al., 2007;Nick et al., 2010).This relation initiates iceberg calving when the combined depth of surface and bottom crevasses reach a threshold percentage of ice thickness (See Section 2.2).Crevasse depth primarily depends on the stress field of the ice shelf, with extensional stresses providing favourable conditions for crevasse opening, though meltwater hydrofracture may also increase calving rates, a process that has been strongly implicated in the 2002 collapse of 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 (Barrand et al., 2013).

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-1.5 mm by 2100 and 0.6-1.6 mm to 2300 (Figure 4a).The sea-level curve rises in the first two decades in response to 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) is restricted to five outlet glaciers in the southern part of the embayment (Figure 5).In contrast, immediate collapse of George VI Ice Shelf perturbs upstream grounded tributaries by up to 0.8 m a −1 , and results in 4-11 mm total sea-level rise (Figure 4b).The more dramatic response of George VI tributary glaciers means that they have not yet reached steady-state by 2100 (Figure 4d), and PSU3D simulations continue to contribute to sea level well beyond this date.Some of this discrepancy between sheet-shelf models is attributed to the fact that PSU3D is not as close to steady-state as BISICLES following initialisation and spin-up (Figures 3,A4).Moreover, ice-sheet 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 similar predicted grounding-line retreat (Table 1).Such a response has been previously attributed to differences in the underlying model physics (L1L2, A-HySSA, (Pattyn et al., 2013)). 9 The Cryosphere Discuss.

Quad=Quadratic
While there is a notable grid dependence in BISICLES projections (Figure A4), this is absent in the PSU3D projections supporting the findings of previous modelling studies (Pollard and DeConto, 2012a;Pattyn et al., 2013) that ice-sheet models with the implementation of an internal flux boundary condition are less sensitive to grid resolution.The required first order convergence (Cornford et al., 2016) of the sea-level rise projections in the BISICLES simulations is met for simulations at 1 km 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 sealevel 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 (Figures 5,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 (Figures 5,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 (Figure 1b), conditions expected to be favourable for marine ice sheet instability.Despite this, grounding-line retreat of 10 George VI outlet glaciers is limited to a few locations and <15 km in length (Figure 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.

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 dependendent 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 combined provide up to 23 mm sea-level equivalent ice loss by 2300 (Figure 7b), with 95% of this total coming from George VI tributary glaciers.

5
The contribution to the sea-level budget from Larsen C embayment glaciers is small (<1.5 mm) and remains equivalent to the fixed calving front simulation, immediate shelf-collapse scenario (Experiment 1).The sea-level commitment from Larsen C glaciers is modest as complete shelf collapse is not forecast until 2150 in RCP8.5 (Figure 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 just five outlet glaciers in the southern part of the embayment (Figure 8a).The larger grounded area loss simulated 10 with PSU3D for Larsen C (Figure 7c) is not a response to 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 on sea-level projections is small (0.28 mm sea-level equivalent). 12 The Cryosphere Discuss.Although projections for Larsen C Ice Shelf glaciers agree 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 (Figure 7b).Experiment 2 (dynamic calving front) simulations with PSU3D provide very similar sea-level projections to Experiment 1 runs (immediate ice-shelf collapse; 6.8-7.1 mm from RCP4.5 and 8.5, respectively).In contrast, BISICLES projects no sea-level rise under RCP4.5, as the amount of meltwater available for hydrofracturing is insufficient to initiate ice-shelf collapse or retreat due to the different implementation of the calving law (See Section 2.2).Under RCP8.5, break-up of George VI Ice Shelf occurs at approximately 2100 (Figure 7d), resulting in widespread grounding-line retreat and sea-level rise of 22 mm by 2300. 13 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-21Manuscript under review for journal The Cryosphere Discussion started: 6 February 2018 c Author(s) 2018.CC BY 4.0 License.
Both sheet-shelf models project similar sea-level rise by the mid 22nd century in Experiment 2. In the BISICLES simulation, however, grounding lines are forced further back into the marine-based sectors (Figure 8d).This enhanced retreat is absent from Experiment 1 because a fixed calving front is used.In Experiment 2 however, increased rates of calving permit the grounding line to retreat much further inland.As this enhanced grounding-line retreat is only present in Experiment 2, it suggests that this retreat is most likely due to a combination of the dynamic calving front and the marine ice-sheet instability mechanism.
Even with a dynamic calving front, enhanced grounding-line retreat is not triggered before some time after ice-shelf collapse (>15 years, Table 2), indicating that fast grounding-line retreat is not triggered before the calving front reaches a retrograde sloping bedrock topography.As a result of widespread grounding-line retreat (Figure A5), extensive dynamic thinning occurs (>1 m a −1 ), extending up to 100 km inland in the southern parts of the embayment (Figure 8d).PSU3D simulations do not show enhanced grounding-line retreat in this sector.We attribute this to the high levels of basal drag inferred by the inversion procedure (Figure 2).This higher basal drag prevents outlet glaciers from speeding up enough to cause dynamic thinning such that the glacier fronts reach flotation, thus limiting the extent of deep basal crevasses, which form only in ice close to or at flotation.As a consequence PSU3D never retreats far enough into the marine-based sectors where the marine ice sheet instability could cause further retreat.

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 (Paolo et al., 2015) 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 subject to 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 sea-level rise projections to 2300 (1.4-1.6 mm for Larsen C embayment glaciers, and 4-6 mm for George VI glaciers, respectively; Figure 4).Projections to 2300 increase by a factor of two for simulations using a Coulomblimited 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 (Table 2).
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 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-21Manuscript under review for journal The Cryosphere Discussion started: 6 February 2018 c Author(s) 2018.CC BY 4.0 License.
(<15%), a more recent higher-resolution dataset (Huss and Farinotti, 2014) provides an increase of ∼100% in ice volume below sea level.When incorporated into ice-sheet model simulations, this boundary condition results in larger grounding-line retreat rates for some basins occupying deeper bedrock troughs.A consequence of this is a sea-level rise projection for Experiment 1 (immediate ice-shelf collapse) with the reference sliding law (Equation 1) that increases by a factor of ∼3 (4.2 mm) for the Huss and Farinotti (2014) boundary input dataset, underlining the significance of accurate boundary dataset for sea-level rise projections.

Comparison with Larsen B Ice Shelf collapse response
To further assess the impact of ice-shelf break-up, five drainage basins from each ice-shelf embayment 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 ice-shelf removal is short lived (∼15 yrs) for both models.Maximum speed-up of ∼300% is possible, though the mean maximum speed-up is ∼50% (Table 2).
These values are smaller than those observed following Larsen B collapse with a maximum of 8fold 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 greatest 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).
For Experiment 2, total maximum grounding-line retreat rates remain similar for PSU3D in comparison to Experiment 1, while they almost triple for George VI in the BISICLES simulation (21.3 km in Experiment 2, 6.4 km in Experiment 1) in agreement with our computed sea-level rise projections (Figure 7).Due to the gradual loss of buttressing, significant speed up is absent in the years following ice-shelf removal across all basins.The overall dynamic response is less vigorous than in Experiment 1 with the exception of the George VI embayment in the BISICLES simulation where larger grounding-line retreat rates lead to high mass loss rates.However, unlike in Experiment 1, most of grounding-line retreat and mass loss occur some time after ice-shelf collapse (>15 years after ice-shelf removal).12 mm sea-level equivalent water in Experiment 1, to 6-22 mm sea-level equivalent water 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.These projections are relatively insensitive to increased ocean forcing,

Figure 1 .
Figure 1.(a) Location map of the Antarctic Peninsula including locations of Larsen C and George VI ice shelves and localities mentioned.
where F C is the negative of the modelled thickness field change when the model is run forward a single time step.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 6 The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-21Manuscript under review for journal The Cryosphere Discussion started: 6 February 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 2 .
Figure 2. Inferred basal traction fields C for the Larsen C (a,b) and George VI model domains (c,d).Black line denotes modelled drainage basins.

Figure 4 .
Figure 4. Upper panel (a,b) shows SLR projections from Experiment 1 (immediate shelf collapse) for BISICLES (solid lines), PSU3D (dashed lines), and BAS-APISM (dotted line).Lower panel (c,d) shows the derivative (rate of change) of the corresponding SLR projections in the upper panel (a,b).Grey shading displays uncertainty associated with SLR projections from BAS-APISM.Note different y-axis scales.

Figure 7 .
Figure 7. SLR projections from Experiment 2 (dynamic calving front) for Larsen C (a) and George VI ice shelves (b) with corresponding area loss of grounded ice (c,d) and ice shelf area loss (e,f).MIROC and CSIRO denote selected global climate model forcing.

Figure 8 .
Figure 8. Dynamic thinning (dh/dt pattern from Experiment 2 (dynamic calving front) averaged over the simulation period 2000-2300 in the RCP8.5 scenario for Larsen C (a,c) and George VI embayments (b,d).Black line denotes modelled drainage basins.LarI-V and GeoI-V indicate drainage basins selected for analysis.
most important contributor to the global sea-level budget to 2300 from Antarctic Peninsula ice-shelf-ice sheet dynamics are glaciers in western Palmer Land feeding George VI Ice Shelf.Our envelope of sea-level rise projections ranges from 4-The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-21Manuscript under review for journal The Cryosphere Discussion started: 6 February 2018 c Author(s) 2018.CC BY 4.0 License.
yet are highly sensitive to changes in the basal boundary condition and the choice of 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 of 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 are limited to less than 4.3 mm by 2300 (0.6-4.3 mm for Experiment 1; 0.4-1.5 mm for Experiment 2).Individual drainage basin analysis indicates a wide range of responses in response 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 Larsen B Ice Shelf.Appendix A The Cryosphere Discuss., https://doi.org/10.5194/tc-2018-21Manuscript under review for journal The Cryosphere Discussion started: 6 February 2018 c Author(s) 2018.CC BY 4.0 License.

Figure A1 .
Figure A1.Near-surface temperature bias for the baseline period 1980-2005 in GCMs for RCP4.5 projections.Dashed black line indicates multi-model mean (-3.1±2.0°C).The selected forcing is highlighted by the red box.

Figure A2 .
Figure A2.Near-surface temperature bias for the baseline period 1980-2005 in GCMs for RCP8.5 projections.Dashed black line indicates multi-model mean (-2.8±1.7°C).The selected forcing is highlighted by the red box.

Figure A3 .
Figure A3.BISICLES L-curve analysis to select Tikhonov parameters λ φ and λC : (a) 3-D scatter plot of the model-data misfit Jm as a function of the regularisation terms J reg C and J reg φ .(b) 2-D cross section for variable λ φ and λC fixed at 10 −1 Pa −2 m 6 a −4 .(c) Reverse case with constant λ φ at 10 9 m 4 a −2 and λC varying.The units of Jm and J reg C are m 4 a −2 and Pa 2 m −2 a 2 , respectively.J reg φ is unitless.Selected

Figure A4 .
Figure A4.Observed ice thickness distribtution (BEDMAP2(Fretwell et al., 2013)) for Larsen C (a) and George VI (d) embayments and modelled ice thickness distribution after spin-up for Larsen C (b,c) and George VI (e,f) embayments.Black line denotes modelled drainage basins.

Figure A7 .
Figure A7.Comparison of modelled grounding-line positions using BISICLES with different basal sliding laws for Larsen C (a) and George VI embayments (b).Black line denotes modelled drainage basins.
Figure 3. ∂V ∂t spin-up plot for BISICLES (solid lines) and PSU3D (dashed lines) at different horizontal resolutions.