Century-scale simulations of the response of the West Antarctic Ice Sheet to a warming climate

(1)

www.the-cryosphere.net/9/1579/2015/

doi:10.5194/tc-9-1579-2015

Century-scale simulations of the response of the West Antarctic Ice Sheet to a warming climate

S. L. Cornford¹, D. F. Martin², A. J. Payne¹, E. G. Ng², A. M. Le Brocq³, R. M. Gladstone¹, T. L. Edwards¹, S. R. Shannon¹, C. Agosta^4,5, M. R. van den Broeke⁶, H. H. Hellmer⁷, G. Krinner⁴, S. R. M. Ligtenberg⁶, R. Timmermann⁷, and D. G. Vaughan⁸

1Centre for Polar Observation and Modelling, School of Geographical Sciences, University of Bristol, Bristol, BS8 1SS, UK

2Computational Research Division, Lawrence Berkeley National Laboratory, Berkeley, California, USA

3Geography, College of Life and Environmental Sciences, University of Exeter, Exeter, EX4 4RJ, UK

4UJF-Grenoble 1/CNRS, Laboratoire de Glaciologie et Géophysique de l’Environnement (LGGE) UMR5183, 38041 Grenoble, France

5Département de Géographie, Université de Liège, Liège, Belgium

6Institute for Marine and Atmospheric Research, Utrecht University, Utrecht, the Netherlands

7Alfred-Wegener-Institut für Polar und Meeresforschung, Bussestrasse 24, 27570 Bremerhaven, Germany

8British Antarctic Survey, Madingley Road, Cambridge, CB3 0ET, UK Correspondence to: S. L. Cornford (s.l.cornford@bristol.ac.uk)

Received: 9 February 2015 – Published in The Cryosphere Discuss.: 23 March 2015 Revised: 14 July 2015 – Accepted: 27 July 2015 – Published: 18 August 2015

Abstract. We use the BISICLES adaptive mesh ice sheet model to carry out one, two, and three century simulations of the fast-flowing ice streams of the West Antarctic Ice Sheet, deploying sub-kilometer resolution around the grounding line since coarser resolution results in substantial underes- timation of the response. Each of the simulations begins with a geometry and velocity close to present-day observations, and evolves according to variation in meteoric ice accumulation rates and oceanic ice shelf melt rates. Future changes in accumulation and melt rates range from no change, through anomalies computed by atmosphere and ocean models driven by the E1 and A1B emissions scenarios, to spatially uniform melt rate anomalies that remove most of the ice shelves over a few centuries. We find that variation in the resulting ice dynamics is dominated by the choice of initial conditions and ice shelf melt rate and mesh resolution, although ice accumulation affects the net change in volume above flotation to a similar degree. Given sufficient melt rates, we compute grounding line retreat over hundreds of kilometers in every major ice stream, but the ocean models do not predict such melt rates outside of the Amundsen Sea Embayment until after 2100. Within the Amundsen Sea Embayment the largest single source of variability is the onset of sustained retreat

in Thwaites Glacier, which can triple the rate of eustatic sea level rise.

1 Introduction

The present-day West Antarctic Ice Sheet (WAIS) is experi- encing an imbalance between the mass it receives as snowfall and that which it loses through discharge to the oceans (Rignot, 2008; Pritchard et al., 2009; Shepherd et al., 2012;

Mouginot et al., 2014; Rignot et al., 2014). In several areas this has led to the persistent loss of ice amounting to a significant contribution to sea-level rise. Continued acceleration of these losses would imply a significant additional global sea-level rise in the coming decades and centuries. Physi- cally based projections of the contribution of the WAIS to sea level rise are hampered by two main factors. The first of these is the lack of a fully coupled climate and ice sheet model, in which the principal forcing on the ice sheet (accumulation at the upper surface and submarine ice shelf melt) is determined within the model. The second is the technical difficulty involved in calculating the flow of ice across the ice sheet’s grounding line and the consequent grounding line

(2)

migration. Progress has been made in both areas (Goldberg et al., 2012a, b; Pattyn et al., 2013), but the computational ex- pense of fully-coupled ice/ocean models at sufficient resolution and for sufficient integration times remains prohibitive.

We approximate full coupling between the ice sheet and the rest of the climate system by imposing combinations of published meteoric accumulation and oceanic melt rate anomaly data on the BISICLES adaptive mesh ice sheet model (Cornford et al., 2013). Two emission scenarios are included: SRES E1, a mitigation scenario in which emissions are stabilized by 2050 at 500 ppm CO2, and A1B, a balanced scenario close to the center of the SRES range.

These were used to drive global climate warming in the UKMO HadCM3 and MPI ECHAM5 global climate models, which have among the highest skill scores in the CMIP3 model group (based on Antarctic SMB, surface air temperature, mean sea level pressure, and height and temperature at 500 hPa, Connolley and Bracegirdle, 2007). The resulting global climate projections provided boundary conditions to two high-resolution atmosphere models: RACMO2 (Ligten- berg et al., 2013) and LMDZ4 (Agosta et al., 2013) and two ocean models: the medium resolution BRIOS (Bremerhaven Regional Ice–Ocean Simulations) (Hellmer et al., 2012) and the higher resolution FESOM (Finite-element Sea ice-Ocean Model) (Timmermann and Hellmer, 2013), ultimately pro- viding seven sets of meteoric accumulation data and eight sets of oceanic melt rate data.

At the same time, we examine the response of the ice sheet model to variability beyond the scope of the atmosphere and ocean models. The climate projections described above tend to agree on the timing and magnitude of future accumulation and melt rate increases, if not the distribution. We com- plement them with some simplified, widespread melt rate increases, as well as projections further into the future, in order to investigate the additional response to more extreme scenarios. The century-scale evolution of the ice sheet model is also sensitive to its present-day state, especially in the Amundsen Sea Embayment, and we evaluate at least a part of this sensitivity – which will prove to be substantial – by varying the initial accumulation rate and hence the initial thinning rate.

In summary, the aim of this paper is to consider the response of the West Antarctic ice streams to process-based and simplified projections of future ocean and atmosphere warming over the 21st and 22nd centuries. We focus on West Antarctica primarily because of constraints on available computational resources; however these areas are also thought to be most vulnerable to future grounding line retreat because of their deep bedrock and changes in oceanic forcing (Hellmer et al., 2012; Pritchard et al., 2012; Ross et al., 2012; Joughin et al., 2014).

2 Methods

2.1 Model equations

BISICLES employs a vertically integrated ice flow model based on Schoof and Hindmarsh (2010) which includes lon- gitudinal and lateral stresses and a simplified treatment of vertical shear stress which is best suited to ice shelves and fast-flowing ice streams. Ice is assumed to be in hydrostatic equilibrium so that given bedrock elevationband ice thick- nesshthe upper surface elevationsis

s=max

h+b,

1− ρ_i ρ_w

h

, (1)

in whichρiandρware the densities of ice and ocean water.

The ice thicknesshand horizontal velocityusatisfy a two- dimensional mass transport equation

∂h

∂t + ∇ ·[uh]=a−M, (2) and two dimensional stress-balance equation

∇ ·[φhµ (2¯ ˙+2tr()I)]˙ +τ^b=ρigh∇s, (3) together with lateral boundary conditions. The terms on the right hand side of Eq. (2),a andM, are the meteoric accumulation rate, applied to the upper surface of the entire ice volume, and the oceanic melt rate, applied to the under-side of ice shelves. When Eq. (2) is discretized, oceanic melt is applied only to cells whose center is floating. As for Eq. (3), ˙is the horizontal strain-rate tensor,

˙=1 2

∇u+(∇u)^T

(4) and I is the identity tensor. The vertically integrated effective viscosityφhµ¯ is computed by evaluating the integral

φhµ(x, y)¯ =φ

s

Z

s−h

µ(x, y, z)dz (5)

numerically, with the ice sheet sub-divided into 10 lay- ers, narrowing progressively from 0.16hnear the surface to 0.03hnear the base. The vertically varying effective viscos- ityµ(x, y, z)includes a contribution from vertical shear and satisfies

2µA(T )(4µ²˙²+ |ρig(s−z)∇s|²)^(n−1)/2=1, (6) where the flow rate exponentn=3,φis a stiffening factor (or, equivalently,φ⁻ⁿ is an enhancement factor), andA(T ) depends on the ice temperatureT through the Arrhenius law described by Hooke (1981),

A(T )=A₀exp 3f

[T_r−T]^k − Q RT

, (7)

(3)

where A0=0.093 Pa⁻³a⁻¹, Q/R=9.48×10³K, f = 0.53 K^k,k=1.17 and Tr =273.39 K. The coefficient φ is estimated by solving an inverse problem (see Sect. B1), and it is simply a convenient way to represent several conflated factors: uncertainty in both temperatureT and the form of A(T ), macroscopic damage, and fabric formation.

Finally, the basal traction is determined by a viscous law:

τ^b=

(−C|u|^m−1u if _ρ^ρⁱ

wh >−b

0 otherwise , (8)

with m=1. Like φ, C will be determined by solving an inverse problem, as described in Sect. B1. Our choice of a linear viscous law may well bias our results toward exces- sive grounding line retreat: in previous work on Pine Island Glacier non-linear laws withm <1 have led to slower rates of retreat (Joughin et al., 2010; Favier et al., 2014).

We hold the fieldsC andφconstant throughout our simulations. That is not to say that these fields ought not change over the course of one or two centuries; for example regions of damage (low φ) might well propagate with the grounding line as englacial stresses grow in regions previously dominated by the balance between gravitational and basal shear stress. Rather, we lack models of sufficient skill for the present, and anticipate incorporating progress in damage models (Borstad et al., 2012) and hydrology models (Werder et al., 2013) in future calculations. We note, however, that the maps ofCandφwe use (see Sect. 2.4.1) already feature slippery beds and weak shear margins hundreds of kilometers upstream from the grounding line.

2.2 Model domains and boundary conditions

We carried out calculations on three rectangular domains, shown in Fig. 1. The largest of these (RISFRIS) covers the Ross and Filchner–Ronne ice shelves and their tributary ice streams, while two smaller domains cover the Amundsen Sea Embayment (ASE) and Marie-Byrd Land (MBL). Each of the rectangular domains is split into an active region _V, where ice is permitted to flow, and a quiescent region _Q where ice is taken to be stationary. For example, in the RIS- FRIS domain,_V covers the present-day drainage basins of the Ross and Filchner–Ronne ice shelves, and the ice shelves themselves, while Q covers the Amundsen Sea Embay- ment, Marie-Byrd Land, and part of the Antarctic Peninsula.

Likewise, inside the ASE domain V spans the drainage basin of Pine Island, Thwaites, Smith, Pope, and Kohler glaciers. This construction assumes that the ice divides will not stray from their current configuration, and so limits us to simulations over a few centuries.

Reflection boundary conditions were applied at the edge of each domain. Ifnis normal to a boundary andtis parallel to it,

u·n=0, t· ∇u·n=0, ∇h·n=0. (9)

MBL

1000 km RISFRIS

1 ASE

2 4

5 6

8 3

7

1. Pine Island Glacier 2. Thwaites Glacier

3. Crosson and Dotson ice shelves 4. Evans Ice Stream

5. Carson inlet and Rutford ice streams 6. Moller and Institute ice streams 7. Foundation Ice Stream 8. Siple Coast

Figure 1. West Antarctica divided into three computational do- mains. Simulations are carried out in three rectangular model domains: RISFRIS, ASE and MBL. Each of these has an active region

_V bounded by the dashed contours and the calving front (black), while the remaining area_Qis made quiescent. Integration of, say, volume above flotation is carried out only over the active regions.

Sea level rise results are given separately for the Filchner–Ronne Ice Shelf (FRIS) and Ross Ice Shelf (RIS) regions, but the simulations are carried out on a domain joining both regions together (RISFRIS). The Amundsen Sea Embayment (ASE) and Marie-Byrd Land (MBL) are simulated separately.

In practice, these boundary conditions are unimportant because of the presence of quiescent regions and calving fronts inside the domain. In the quiescent regions, we set the basal traction coefficient to a large value,C=10⁵Pa m⁻¹a so that at the interface between_Vand_Q,

u≈0 and uh≈0, (10)

while at the calving front (which is fixed), we impose the usual conditions on the normal and transverse stress:

n·[φhµ (2¯ ˙+2tr()I)]˙ =1 2ρig

1− ρ_i

ρ_w

h²n. (11) These boundary conditions, and indeed, Eq. (11) alone for a problem whose entire boundary is a fixed calving front, are sufficient provided thath(x, y, t=0) is given and that the basal friction coefficientC(x, y)is non-zero in at least part of the ice sheet.

2.3 Adaptive mesh refinement

Fine horizontal resolution, or other careful treatment, is held to be crucial when simulating grounding line migration (Vieli and Payne, 2005; Durand et al., 2009). Indeed, the

(4)

1000 km

200 km

< 1 5 30 166 910

> 5000 u (^ma⁻¹)

(a)

1000 km

200 km

< −1000 −600 −200 200 600

> 1000 uo−u (^ma⁻¹)

(b)

1000 km

200 km

<10⁰ 10¹ 10² 10³ 10⁴

>10⁵ C(Pa m⁻¹a)

(c)

1000 km

200 km

<10⁴ 10⁵ 10⁶ 10⁷ 10⁸

>10⁹

φµ (^{Pa a})

(d)

1000 km

200 km

< −10 −6 −2 2 6

> 10 a0−M0(^{m a}⁻¹)

(e)

1000 km

200 km

< −10 −6 −2 2 6

> 10

∂h∂t(^{m a}⁻¹)

(f)

1

Figure 2. Model initial state. Panels show (a) the ice flow speed|u|, (b) the difference between observed and model speed|u_o| − |u|, (c) the basal traction coefficientc, (d) the vertically averaged effective viscosityφµ, (e) the synthetic mass balance¯ a₀(x, y)−M₀(x, y, t=0), and (f), the thickening rate ^∂h_∂t, all at the start of the prognostic calculations. Data from all three domains are shown, with a magnified inset showing the ASE data in more detail.

(5)

BISICLES ice sheet model was designed primarily with this in mind, and discretizes the stress and mass balance Eqs. (2) and (3) on block-structured meshes built from rectangular subsets of uniform grids with resolution1x^`, with 0≤`≤L and 21x^`+1=1x^`. While restrictive in some senses – all model domains must be rectangular, for example – these meshes have a signal advantage: it is straightforward to gen- erate new meshes as the ice sheet evolves, and to transfer the previous time-step’s ice thickness data to the new mesh in a conservative fashion. It is also relatively easy to study convergence with mesh resolution by running the same experiment for successive values of L and verifying that the differences between, say, the volume above flotation calcu- lated in each case converge to zero at the expected rate. We regard such a convergence study as a pre-requisite for any ice dynamics simulation, since there is no general proof that any particular mesh resolution is adequate. We include the relevant results in appendix A, where we show that sub- kilometer resolution around the grounding line is necessary and adequate in all of West Antarctica, but that finer resolution is needed in the Amundsen Sea Embayment, where we employ a mesh with 250 m≤1x^`≤4000 m, than in the Ronne-Filchner and Ross ice shelf catchments; there we find that a mesh with 625 m≤1x^`≤5000 m is sufficient.

There is a relationship between horizontal mesh spacing, ice velocity and time-step. Since the advection scheme chosen to evolve Eq. (2) is explicit, we re-compute the time-step periodically so that a Courant-Friedrichs-Lewy (CFL) con- dition1t <0.25^1x_|u|^l is satisfied everywhere in the domain.

In practice, this leads to as many as 128 time-steps per year in the ASE domain, where the mesh is finest and the flow is fastest.

2.4 Model data requirements

Time-dependent simulations require initial ice thickness data h₀(x, y) as well as accumulation rates a(x, y, t ) and melt rates M(x, y, t )for Eq. (2), together with a bedrock elevation map b(x, y), a basal friction coefficient field C(x, y), a temperature fieldT (x, y, z)and a stiffening factorφ (x, y) to solve Eq. (3). Bedrock elevation and initial ice thickness data for the RISFRIS and MBL domains were taken from the ALBMAP DEM (Le Brocq et al., 2010). ALBMAP is provided at a lower resolution (5 km) than the more recent Bedmap2 (1 km) (Fretwell et al., 2013), but the model is less sensitive to finer-scale variations in bedrock (Sun et al., 2014) and in any case the distance between the flight-lines measur- ing ice thickness for both DEMs is typically not finer than 5 km. A custom map of bedrock elevation and ice thickness set on a 1 km grid was used for the ASE domain: it is close to the Bedmap2 data, and was used before for studies of Pine Island Glacier (Favier et al., 2014). It was prepared in a similar manner to ALBMAP, but includes extra data from high- resolution airborne radar (Vaughan et al., 2006) and submarine surveys (Jenkins et al., 2010). It also includes a pinning

point at the tip of Thwaites Glacier’s slower flowing eastern ice shelf, a feature that is clearly visible in the velocity data (Joughin et al., 2009; Rignot et al., 2011), that corresponds to peak one of the two described in Tinto and Bell (2011), but is absent in the bathymetry data. We raised the bathymetry by 120 m to ground the ice in that region. Ice temperature data are provided by a three-dimensional thermo-mechanical model (Pattyn, 2010) and is held fixed in time.

The basal friction and stiffening coefficients are chosen by solving an inverse problem similar to those of MacAyeal (1993), Joughin et al. (2009) and Morlighem et al. (2010).

A detailed description is given in Appendix B1, but in summary, we construct smooth fieldsC(x, y)and φ (x, y) that minimize the mismatch between the modeled speed and the published InSAR observations. For the RISFRIS and MBL domains we use InSAR observations made between 2007 to 2009 (Rignot et al., 2011). For the ASE domain, which has accelerated over the last decade, we use measurements made in 1996 (Joughin et al., 2009). However, the observed speeds are not entirely compatible with the thickness and bedrock data. Notably, computing the flux divergence from the thickness and velocity data results in a region of 100 m a⁻¹thickening across Pine Island Glacier’s grounding line. Others have noted this strong thickening, and address it by imposing a large synthetic mass balance (Joughin et al., 2010), by constraining the ice viscosity and accepting a worse match to the observed velocity (Favier et al., 2014), or by modifying the bed to give acceptable thickening rates while matching the observed velocity field (Rignot et al., 2014; Nias et al., 2015). Here, we soften the ice around the grounding line by reducing the stiffening factor relative to the field in the inverse problem.

The accumulation and melt rates are computed by adding future climate anomalies, described in Sect. 2.5, to initial accumulation and melt ratesa₀andM₀chosen to hold the ice sheet close to equilibrium. Determination ofa₀ and M₀ – especiallyM₀ – is somewhat involved, and is described in more detail in Appendix B2. Essentially, we evolve the ice sheet geometry for 50 years while holding the ice shelf constant in order to dampen short-wavelength, large-amplitude fluctuations in the flux divergence. At the end of this relax- ation period, we computea0andM0from∇ ·[uh], withM0

parametrized as a function of time and space so that peak melt rates follow the grounding line as it migrates. Since a synthetic mass balance along the lines ofa0 will tend to counter the thinning that is already evident in the ASE, we carry out some additional experiments wherea0is replaced by an accumulation pattern derived from the RACMO atmosphere model.

2.4.1 Initial state

Figure 2 illustrates the state of the three regional models at the end of the initialization procedure. The large ice streams flowing into the Amundsen Sea, through the Filchner–Ronne

(6)

t (a)

∆a(ma−1, 32 year moving average)

2000 2050 2100 2150 2000 2050 2100 2150

−0.2 0 0.2

Amundsen Sea Embayment West Antarctica

−0.1

−0.05 0 0.05 0.1

−0.2 0 0.2

Marie−Byrd Land

−0.05 0 0.05

Ross Ice Shelf

−0.05 0 0.05

Filchner−Ronne Ice Shelf

GCM / Scenario HadCM3 / A1B HadCM3 / E1 ECHAM5 / A1B ECHAM5 / E1 Atmosphere model

RACMO LMDZ4

Figure 3. Accumulation anomalies integrated over each region. The atmosphere models provide enhanced accumulation only for the A1B emission scenarios.

t (a)

∆M(ma−1)

2000 2050 2100 2150 2000 2050 2100 2150

−10 0 10 20 30

0 2 4

0 10

Marie−Byrd Land

−2 0 2

Ross Ice Shelf

−2 0 2

Filchner−Ronne Ice Shelf Ocean model

FESOM BRIOS

GCM / Scenario HadCM3 / A1B HadCM3 / E1 ECHAM5 / A1B ECHAM5 / E1

Figure 4. Melt rate anomalies integrated over each region. In con- trast with the atmosphere models, the ocean models provide simi- larly growing melt rates in both A1B and E1 scenarios. Note that the Amundsen Sea Embayment and Marie-Byrd Land melt rate anomalies are not given directly by the ocean models, which do not resolve the smaller ice shelves, but are characterized in terms of nearby Cir- cumpolar Deep Water temperatures.

Ice Shelf, and through the Ross Ice Shelf are apparent in both the flow field and the basal traction fieldC. The flow field it- self is close to the observed speed, with the largest mismatch due to the softening in Pine Island Glacier described above, and within 200 m a⁻¹of the observations elsewhere. Regions of fast flow are fringed by weak shear margins that are due to both the shear thinning action of Glen’s flow law and localized low stiffening factorφ. Note that the observed velocity field cannot be matched withφ=1, whatever the basal traction. For example the shear margins in Pine Island Glacier and the division between the western and eastern portions of Thwaites Glacier’s ice shelf requireφ∼0.1.

The thickening rate (given a0 and M0) is between −5 and 5 m a⁻¹except at calving fronts. Integrating this thickening rate leads to an annual loss of volume above flotation of 3 km³a⁻¹ in the ASE, 7 km³a⁻¹ in the Filchner–Ronne ice shelf basin and 7 km³a⁻¹ in the Ross ice shelf basin.

The synthetic accumulation a0 used to obtain this thickening rate does include some unrealistic large-amplitude short- wavelength features, with the largest values in mountainous regions with steep slopes: the ring-shaped features in the ASE surround isolated peaks, for example. Strong ablation is limited to the ice shelves, witha₀between−5 and 5 m a⁻¹: in particular there is no region of∼100m a⁻¹ablation needed to counter the flux arriving at the Pine Island Glacier grounding line.

2.5 Prognostic experiments

Twenty-two simulations were performed for one or more of the three model domains. Each simulation makes use of the same initial geometry, basal traction coefficient, and stiffness coefficient, but differs from the others in terms of the meteoric accumulation and oceanic melt rates imposed. Each experiment is named after these forcing data, and falls into one of three groups: two control calculations, which are subject to a constant climate, fourteen experiments forced by combinations of time-dependent climate model data, and six melt rate anomaly experiments, which are subject to constant-in-time accumulation. The experiments are summarized in Table 1 and described in detail below.

2.5.1 Combined anomaly experiments

Future climate forcings were derived from the atmosphere and ocean models by computing space- and time-dependent anomalies with respect to the 1980–1989 mean, and adding them to a0(x, y) andM0(x, y, t ). By combining the seven atmosphere projections with the eight ocean projections, we have fourteen experiments, as shown in Table 1. These are named after the anomalies: for example, the experiment named H/A/R/F combines the HadCM3/A1B/RACMO2 accumulation anomalies with the HadCM3/A1B/FESOM melt rate anomalies. Given the fourteen forcing combinations, the ice sheet model was evolved, starting from its initial state in

(7)

Table 1. Accumulation data, oceanic melt rate data, and duration for all 22 experiments. Each experiment is named for its climate forcing data, which are specified by combinations of GCM, emission scenario, high-resolution atmosphere model and high-resolution ocean model.

Synthetic Mass balance anomalies

Experiment accumulation GCM Emissions Atmosphere Ocean Final year

control a₀(x, y) – – – – 2200

control⁰⁰ a₀⁰⁰(x, y) – – – – 2200

H/A/R/F a₀(x, y) HadCM3 A1B RACMO2 FESOM 2200

H/A/R/B a₀(x, y) HadCM3 A1B RACMO2 BRIOS 2200

H/A/L/F a₀(x, y) HadCM3 A1B LMDZ4 FESOM 2200

H/A/L/B a₀(x, y) HadCM3 A1B LMDZ4 BRIOS 2200

H/E/R/F a₀(x, y) HadCM3 E1 RACMO2 FESOM 2150

H/E/R/B a₀(x, y) HadCM3 E1 RACMO2 BRIOS 2200

H/E/L/F a₀(x, y) HadCM3 E1 LMDZ4 FESOM 2150

H/E/L/B a₀(x, y) HadCM3 E1 LMDZ4 BRIOS 2200

E/A/R/F a₀(x, y) ECHAM5 A1B RACMO2 FESOM 2100

E/A/R/B a₀(x, y) ECHAM5 A1B RACMO2 BRIOS 2100

E/E/R/F a₀(x, y) ECHAM5 E1 RACMO2 FESOM 2100

E/E/R/B a₀(x, y) ECHAM5 E1 RACMO2 BRIOS 2100

E/E/L/F a₀(x, y) ECHAM5 E1 LMDZ4 FESOM 2100

E/E/L/B a₀(x, y) ECHAM5 E1 LMDZ4 BRIOS 2100

H/A/0/F a₀(x, y) HadCM3 A1B – FESOM 2300

H/A/0⁰/F a₀⁰(x, y) HadCM3 A1B – FESOM 2200

H/A/0⁰⁰/F a₀⁰⁰(x, y) HadCM3 A1B – FESOM 2200

0/U16 a₀(x, y) – – – 16 m a⁻¹ 2200

0/U8 a₀(x, y) – – – 8 m a⁻¹ 2200

0⁰⁰/U16 a₀⁰⁰(x, y) – – – 16 m a⁻¹ 2200

1980 to at least 2100 and on to 2150 or 2200 if the forcing data were available. The HadCM3/A1B/FESOM ocean data, both sets of HadCM3/A1B/BRIOS ocean data and all of the HadCM3/A1B atmosphere projections were sufficient to run the ice sheet model until 2200, the HadCM3/E1/FESOM data run to 2150, and the ECHAM5 data to 2100.

Neither ocean model produced substantial melt rate increases in the ASE or MBL domains, presumably because they are not able to resolve the small ice shelves along those coasts. We computed melt rates in those regions from projections of nearby ocean temperatures. The melt rates and consequent thinning experienced by small ice shelves, such as Pine Island Glacier, is thought to be forced by changes in the temperature of near-coast water masses (Jacobs et al., 2011; Pritchard et al., 2012). We compute a local ocean temperature anomaly 1T (t )by averaging the projected ocean temperature over volume bounded laterally between the con- temporary ice front, the sector boundaries, and the 1000 m bathymetric contour and vertically between depths of 200 and 800 m, on the grounds that water contributing to melting must be deep enough to interact with the base of an ice shelf but shallow enough to cross the continental shelf break.

Finally, a melt rate anomaly 1M(t )=161T (t )m a⁻¹K⁻¹ was chosen to be at the upper end of the range of observa-

tional and modeling studies (Holland et al., 2008; Rignot, 2002).

The accumulation and melt rate anomalies, plotted in Figs. 3 and 4, have a notable feature. The A1B atmosphere models project increased accumulation during the 21st century, and a further increase during the 22nd century, over and above the E1 models. Although the two atmosphere models distribute snowfall differently, with RACMO2 concen- trating its increased accumulation over the Amundsen Sea Embayment and Filchner–Ronne Ice Shelf drainage basins and LMDZ4 heaping mass over Marie-Byrd Land and the Ross Ice Shelf drainage basin, both models project a three- fold increase for A1B over E1. At the same time, the A1B and E1 ocean models both provide enhanced melt rates from 2100, with the most obvious difference between trends being the choice of FESOM or BRIOS. Even before carrying out any simulations, we can expect to see similar dynamic thinning in the two emissions scenarios, which, coupled with the extra accumulation in A1B, means that we expect to simulate more sea level rise for E1 (mitigation) emissions than A1B (business as usual) emissions.

(8)

t (a)

∆V(103 km3 =2.5 mm s.l.e)

2000 2050 2100 2150 2000 2050 2100 2150

● ● ● ● ● ● ● ● ● ●

−30

−20

−10 0

● ● ● ● ● ● ● ● ● ●

−30

−20

−10 0

● ● ● ● ● ● ● ● ● ●

−10 0 10

Marie−Byrd Land

● ● ● ● ● ● ● ● ● ●

−10 0 10

Ross Ice Shelf

● ● ● ● ● ● ● ● ● ●

−10 0 10

●

Experiment control H/A/R/F H/A/R/B H/A/L/F H/A/L/B H/E/R/F H/E/R/B H/E/L/F

H/E/L/B E/A/R/F E/A/R/B E/E/R/F E/E/R/B E/E/L/F E/E/L/B

Figure 5. Net change in volume above flotation over the course of the combined anomaly experiments. Only the Amundsen Sea Em- bayment experiences a net loss (1V) in all of the combined experiments. Nonetheless, the result is a net loss over West Antarctica as a whole. Note that Thwaites glacier does not retreat in the combined anomaly experiments (which use the synthetic accumulation), and the ASE could contribute an extra 9×10³km³ loss by 2100 and 40×10³km³by 2200. A magnified version of this figure, covering the period 1980-2100, is included in the supplement, as is a .csv file of the data.

2.5.2 Melt rate anomaly experiments

The climate-forced experiments outlined above present a rather limited view of future change. Since both the ocean and atmosphere models project similar futures, they cannot provide much information about the response to earlier or more widely distributed ice shelf thinning. At the same time, the assumption that the ice sheet was in steady state at the end of the 20th century does not allow us to examine changes that may already be under way. A number of experiments with melt rate anomalies but no accumulation anomalies were carried out to address these limitations.

Parts of the ice sheet model might be on the brink of dra- matic change in 2200, so we ran a longer set of calculations, starting in 1980 and running until 2300, based on the HadCM3/A1B/FESOM melt rates. This experiment has only the synthetic accumulation field, so we label it H/A/0/F. The HadCM3/A1B/FESOM melt rates were chosen because they run up to 2200, produce enhanced melting in all of the basins, and rise constantly from 2050 onward to give the largest melt rates at the end of the 22nd century. From 1980 to 2200 we applied the melt rate anomalies as before, and applied the 2200 melt rate anomaly for the remainder of the simulation.

Imposing a synthetic accumulation field to hold the ice sheets close to steady state given the present-day geometry and velocity is a questionable choice in the Amundsen Sea Embayment, where observations over the last decades show extensive thinning. Furthermore, the synthetic mass balance fielda₀(x, y), which we constructed to hold Thwaites Glacier in steady state during the control experiment, includes a spot of unrealistic – 5 m a⁻¹– accumulation close to the Thwaites Glacier grounding line (see Fig. 2). In light of these issues, we carried out three additional simulations. None of these experience accumulation anomalies but the first, H/A/0⁰/F, has a synthetic accumulation fielda₀⁰(x, y)with the 5 m a⁻¹accumulation spot removed, and the second and third, H/A/0⁰⁰/F and control⁰⁰ do not make use of a synthetic accumulation field at all, but employ the HadCM3/E1/RACMO2 1990–1999 temporal mean, which we will call a₀⁰⁰(x, y), from 1980 onward. H/A/0⁰/F and H/A/0⁰⁰/F are subject (like H/A/0/F) to the HadCM3/A1B/FESOM melt rate anomaly data, while the control⁰⁰ experiment maintains the same melt rate parametrization as the control experiment, that is M₀(x, y, t ).

As none of the melt rate anomalies in the ASE exceed 10 m a⁻¹until after 2050, we also examined the model’s response to earlier ocean warming. Two simulations, 0/U16 and 0⁰⁰/U16, were performed, with melt rate anomalies of 16 m a⁻¹ applied across all floating ice from 1980 onward.

0/U16 used the synthetic accumulation fielda0(x, y), while 0⁰⁰/U16 used the HadCM3/E1/RACMO2 1990–2000 mean accumulation fielda⁰⁰₀(x, y).

Both FESOM and BRIOS ocean models produce similar melt rate anomalies, with enhanced melt rates concentrated around Berkner Island in the Filchner–Ronne ice shelf, and around Roosevelt Island in the Ross Ice Shelf. Those similar patterns are due to the physics of ocean circulation in the two models, but it makes sense to consider ice sheet sensitivity to melt rates that cover a greater extent. At the same time, as in the ASE, melt rates begin to grow around 2100 in all of the drainage basins, so we need to consider our sensitivity to earlier warming. With those aims in mind, we carried out a pair of uniform melt rate experiments (0/U8,0/U16) in the RISFRIS domain, applying melt rate anomalies of 8 and 16 m a⁻¹across the entire extent of floating ice starting from 1980.

3 Results and discussion

Combining melt rate and accumulation anomalies leads to essentially the same patterns of dynamic thinning and grounding line retreat as melt rate anomalies alone. For example, the H/A/R/F and H/A/L/F simulations exhibit similar grounding line retreat to the H/A/0/F results. With that in mind, we will discuss the variation in volume above flotation between the combined anomaly experiments in Sect. 3.1 before dis-

(9)

t (a)

∆V−Α∆a(103 km3 =2.5 mm s.l.e)

2000 2050 2100 2150 2000 2050 2100 2150

● ● ● ● ● ● ● ● ● ●

−30

−20

−10 0

● ● ● ● ● ● ● ● ● ●

−30

−20

−10 0

● ● ● ● ● ● ● ● ● ●

−10 0 10

Marie−Byrd Land

● ● ● ● ● ● ● ● ● ●

−10 0 10

Ross Ice Shelf

● ● ● ● ● ● ● ● ● ●

−10 0 10

●

Experiment control H/A/0/F H/A/R/F H/A/R/B H/A/L/F H/A/L/B H/E/R/F H/E/R/B

H/E/L/F H/E/L/B E/A/R/F E/A/R/B E/E/R/F E/E/R/B E/E/L/F E/E/L/B

Figure 6. Dynamic change in volume above flotation over the course of the combined anomaly experiments. Dynamic losses (computed by subtracting the accumulation anomaly from the net loss) occur in all regions, and are nearly independent of the accumulation anomalies, so that the net change in (for example) the H/A/R/F simulation is much the same as the sum of the volume change computed for the ocean-forced H/A/0/F simulation and the HadCM3/A1B/RACMO2 accumulation anomaly. Note that Thwaites glacier does not retreat in the combined anomaly experiments (which use the synthetic accumulation), and the ASE could contribute an extra 9×10³km³loss by 2100 and 40×10³km³by 2200. A magnified version of this figure, covering the period 1980- 2100, is included in the supplement, as is a .csv file of the data.

cussing grounding line migration in the context of the melt rate anomaly experiments in Sect. 3.2.

3.1 Combined anomaly experiments

Only the Amundsen Sea Embayment experiences a net loss of volume above flotation (1V) in all of the combined anomaly experiments (Fig. 5). Both the Ross Ice Shelf drainage basin and Marie-Byrd Land show a positive imbalance, while the Filchner–Ronne Ice Shelf drainage basin remains close to balance. Adding all four trends together, West Antarctica sees a net loss of 0–8×10³km³ by 2100 and 3–23×10³km³ by 2200. Note that Thwaites Glacier does not retreat in these combined calculations, as they all apply the synthetic accumulation. When Thwaites glacier does retreat the ASE model loses an extra 9×10³km³by 2100 and 40×10³km³ by 2200, based on the difference between the H/A/0/F and H/A/0⁰⁰/F melt rate anomaly experiments (see Sect. 3.2.1).

The differences between responses to the ocean models are quantitative rather than qualitative, with the higher BRIOS melt rates leading to faster retreat along the same paths. Fig- ure 5 shows the volume above flotation trends for the combined anomaly experiments where, with everything else held equal, the BRIOS simulations exhibit a 10×10³km³greater loss (−1V) by the end of the 22nd century than the FE- SOM simulations. Around half of this difference is concentrated in the Amundsen Sea Embayment, with the remainder divided between the Filchner–Ronne Ice Shelf region and Marie-Byrd Land.

The difference between RACMO2 and LMDZ4 simulations with a given ocean model and the A1B scenario is as large as the difference between both ocean models across emission scenarios. Although the H/A/R/F and H/A/L/F grounding line retreat is essentially the same, the decrease in volume above flotation over West Antarctica as a whole differs by 10×10³km³, with the majority of that difference accounted for by the larger LMDZ4 accumulation over the drainage basin of the Ross Ice Shelf. The H/A/R/B and H/A/L/B experiments differ in the same way. Variation between the atmosphere models for the E1 scenario is smaller, with all of the E1 models having similar mass loss trends.

The choice of GCM and emission scenario leads to the largest variation between the combined anomaly experiments. Melt rates grow over time in both A1B and E1 scenarios, but accumulation grows much less in the E1 scenario.

The four HadCM3 E1 experiments produce a net volume loss between 6 and 10×10³km³during the 21st century, and around 20×10³km³by the end of the 22nd century (assum- ing that the E1/FESOM experiments follow the same trend from 2150 onward). Of the A1B simulations, only H/A/R/B results in a similar trend, with H/A/R/B and H/A/L/F giving rise to around 10×10³km³loss by 2200 and H/A/R/B less than 5×10³km³. The ECHAM5 E1 and A1B simulations are generally closer to balance, but do not run beyond 2100 when the majority of HadCM3 imbalance occurs.

Despite the wide variation in accumulation anomalies, we see little interaction between atmosphere model and ice dynamics. Figure 6 shows the loss of volume above flotation due to ice dynamics,1Vd. For a given region,1Vdis the difference between the net change in volume above flotation and the cumulative, integrated accumulation anomaly:

1V_d(t )=1V (t )−

t⁰=t

Z

t⁰=1980

Z

1addt⁰. (12)

In each case the difference between curves is dominated by the difference in ocean anomaly. The H/A/R/F and H/A/L/F curves lie close to one another (and to the H/A/0/F curve), each resulting in1V_d≈10×10³km³ (25 mm SLE) across the whole of West Antarctica by 2100 and around 30× 10³km³(75 mm SLE) by 2200. The H/A/R/B and H/A/L/B trends are also close to one another, but lead to rather more

(10)

Figure 7. Grounding line migration and bedrock elevation in the melt rate anomaly (no accumulation anomaly) experiments. Bedrock (b) contours are drawn every 400 from 1200 m below to 1200 m above sea level Pine Island Glacier and the ice streams feeding the Dotson and Crosson Ice Shelves retreat throughout the 21st, 22nd, and 23rd century CE in all simulations apart from the control. Thwaites Glacier retreats over 200 km during the 21st and 22nd centuries when subjected to the 1990s HadCM3/A1B/RACMO2 accumulation (H/A/0⁰⁰/F,0⁰⁰/U16), and still retreats by more than 100 km in the control⁰⁰experiment when no melt rate anomaly is applied, but its retreat is delayed when subjected to a synthetic accumulation field (H/A/0/F,0/U16). In the Filchner–Ronne Ice shelf region, the Möller, Institute, and Evans Ice Streams begin to retreat during the 22nd century when forced by the HadCM3/A1B/FESOM (H/A/0/F) melt rates, and their retreat accelerates in the following century, but all three ice streams, along with the Foundation and Rutford Ice Streams and Carlson inlet, retreat during the 21st century if uniform 16 m a⁻¹melt rates are applied (0/U16). Along the Siple coast, Whillans Ice Stream, and to a lesser extent Mercer Ice Stream, retreat over the 21st and 22nd century in both the H/A/0/F and control simulations, but their grounding lines have merely swept over a lightly grounded area between the model initial state and the present-day state. The MacAyeal and Bindschadler Ice Streams are driven to retreat 100 km during the 22nd century by the H/A/0/F melt rates but all four streams retreat more than 200 km in the 0/U16 experiment. An alternative version of this figure, intended for larger displays, is included in the Supplement.

excess discharge – around 40×10³km³(100 mm SLE) by 2200. Overall, the net1V (t )for a simulation with accumulation anomaly1a(t )and melt rate anomaly1M(t )can be estimated rather precisely from the result,1V⁰(t ), of a simulation with the same ocean anomaly and a different (or no) accumulation anomaly1a⁰(t ):

1V (t )≈1V⁰(t )+

t⁰=t

Z

t⁰=1980

Z

1a−1a⁰

ddt⁰. (13)

This result is valid for the Amundsen Sea Embayment and the Filchner–Ronne and Ross ice shelf drainage basins, and only breaks down in Marie-Byrd Land, which contributes little to the projections. We account for it by noting that, in these century scale simulations, increased melt rates lead to large amplitude but localized thinning, whereas increased

(11)

Time, t (a)

∆V(103 km3 )

2000 2050 2100 2150 2200 2250

−80

−60

−40

−20 0

−200

−150

−99

−50 0

∆V(mm s.l.e)

H/A/0/F 0/U16 control H/A/0'/F H/A/0''/F 0''/U16 control''

Figure 8. Change in volume above flotation (1V (t )) in the Amundsen Sea Embayment during the melt rate anomaly experiments. The ASE discharges an excess volume between 5×10³and 20×10³km³ by 2100, and between 20×10³ and 60×10³km³ by 2200. The difference is dominated by the onset of retreat in Thwaites Glacier and Pine Island Glacier. Pine Island Glacier begins its retreat around 2000 in all simulations, apart from the control and control⁰⁰experiment, as do the glaciers feeding Dotson Ice Shelf and Crosson Ice Shelf. Thwaites Glacier, on the other hand, begins to retreat immediately in the H/A/0⁰⁰/F, control⁰⁰and 0⁰⁰/U16 experiments, in around 2100 in the H/A/0⁰/F and 0⁰/U16, and after 2200 in the H/A/0⁰/F and O⁰/U16 experiments.

accumulation causes low amplitude but widely distributed thickening.

3.2 Melt rate anomaly experiments

The melt rate anomaly experiments (that is, the experiments with no accumulation anomalies) all exhibit grounding line retreat in excess of the control simulation. Fig- ure 7 depicts this retreat for the H/A/0/F experiment, along- side RISFRIS and ASE uniform melt rate simulations (0/U16) and ASE simulations with no synthetic accumulation (control⁰⁰,H/A/0⁰⁰/F and 0⁰⁰/U16). Provided that melt rates are sufficient, deep bedded glaciers flowing into the Filchner–Ronne and Ross Ice Shelves see their grounding lines retreat by as much as 100 km in a century, as do Pine Is- land and Thwaites Glaciers. However, while the ASE retreats during the 21st, 22nd, and 23rd century in both the FESOM and uniform melt experiments, the RISFRIS glaciers do not show significant retreat until the 22nd century when driven by FESOM melt rates (which do not start to grow until the late 21st century). Animations showing several of the melt rate anomaly experiments are included in the supplement.

3.2.1 Amundsen Sea Embayment

The Amundsen Sea Embayment thins throughout the melt rate anomaly simulations, losing 5–14×10³km³ (15–

40 mm SLE) volume above flotation between 2000 and 2100 and 20–70×10³km³ (50–190 mm SLE) by 2200 (Fig. 8).

Pine Island Glacier and the ice streams feeding the Crosson and Dotson Ice Shelves experience ∼1 km a⁻¹ grounding line retreat from the present day onward in all of the experiments apart from control and control⁰⁰, while Thwaites Glacier sees its retreat delayed in some simulations. The major distinction between simulations is the onset of retreat in Thwaites Glacier: experiments in which its grounding line begins to retreat around 2000 lose volume at more than twice the rate of those in which retreat begins around 2200.

Projections of retreat in Thwaites Glacier are strongly af- fected by initial conditions, with some simulations showing little retreat and others shedding between 100 and 210 km³a⁻¹volume above flotation over the 21st and 22nd centuries. In calculations with the synthetic mass balance a0: H/A/0/F, 0/U16, and the control experiment 0/U0), the grounding line remains close to the present-day position until after 2200, despite the near complete removal of its ice shelf by that time in H/A/0/F and 0/U16. It seems that even with the isolated promontory in the slowly flowing eastern section, the ice shelf exerts little back-pressure on the ice stream and so its loss is not felt strongly. The majority of the grounding line retreat only begins after 2200, triggered by the retreat of the small stream which diverts from Thwaites glacier to flow into the south-western corner of Pine Island Glacier and from then on the grounding line retreats at a rate of 1 km a⁻¹until 2300. In contrast, the glacier begins to retreat around 2100 in the H/A/0⁰/F experiment, and around 2000 in the H/A/0⁰⁰/F, 0⁰⁰/U16 and even the control⁰⁰experiments. In these last three simulations, marine ice sheet instability is already acting at the beginning of the simulation, and at no point does the ice shelf provide enough buttressing to prevent it. Provided that retreat is initiated, the higher melt rates applied in H/A/0⁰⁰/F, 0⁰⁰/U16 do result in faster retreat (∼200 km³a⁻¹) than in the control⁰⁰(∼100 km³a⁻¹).

Both mass loss and grounding line migration in Thwaites Glacier accelerate during the second century of retreat. For the H/A/0⁰⁰/F simulation, volume above flotation decreases at a mean rate of 75 km³a⁻¹ (0.2 mm a⁻¹SLE) between 2000 and 2100, while the grounding line retreats at a rate

∼1 km a⁻¹ across a region featuring a broad area less than 800 m below sea level and a narrow trough between 800 and 1200 m below sea level. Over the following century, the grounding line crosses a widening region of deeper bedrock (>1200 m below sea level), so that the greater rate of flow associated with thicker ice at the grounding line is integrated over a broadening front. The average rate of grounding line retreat grows to ∼2 km a⁻¹, and the rate of loss of volume above flotation to 320 km³a⁻¹(0.9 mm a⁻¹SLE).

A similar calculation of accelerating mass loss, with losses of less than 0.25 mm a⁻¹SLE during the 21st century and up to 1 mm a⁻¹SLE thereafter was reported by Joughin et al.

(2014).

(12)

Time, t (a)

∆V(103 km3 )

2000 2050 2100 2150 2200 2250

−80

−60

−40

−20 0

−200

−150

−99

−50 0

∆V(mm s.l.e)

H/A/0/F 0/U16 0/U8 control

Filchner−Ronne Ice Shelf Ross Ice Shelf

Figure 9. Change in volume above flotation (1V (t )) in the Filchner–Ronne and Ross ice shelf regions during the melt rate anomaly experiments. The FESOM-forced simulations (H/A/0/F) lose little volume before 2100, in contrast to the uniform melt rate calculations 0/U8 and 0/U16. All three simulations feature a period where volume is lost from the Filchner–Ronne region at a rate of more than 200 km³a⁻¹– starting immediately for the uniform melt rate experiments but delayed until 2200 for H/A/0/F – which corresponds to the retreat of the Möller and Institute Ice Streams across the Robin Sub-glacial Basin and the flotation of Bungenstock ice rise, which abates after a century. The Ross region, on the other hand, tends to see its rate of mass loss increase throughout the simulation.

Variation between the remaining ASE projections, dominated by Pine Island Glacier, is due to both ocean forcing and initial conditions. Neither the control nor the control⁰⁰simulations exhibit as much grounding line retreat in Pine Island Glacier as those with enhanced melt rates, with the control grounding line holding its initial position and the control⁰⁰ grounding line retreating by around 50 km over 200 years. In contrast, the FESOM-forced simulations all see the grounding line retreat by 60 km in the 21st century and a further 160 km in the 22nd century. The H/A/0⁰⁰F simulation, which employs the RACMO2 surface mass balance rather than a synthetic mass balance, loses volume above flotation at a rate rising from 70 km³a⁻¹over the 21st century, comparable to rates computed in other modeling studies (Seroussi et al., 2014; Favier et al., 2014), and present-day observations (Mouginot et al., 2014). Grounding line retreat slows toward the end of the 21st century around a bedrock rise and stronger bed 60 km upstream from the present-day position (Joughin et al., 2010), after which increase in melt rates from around 2100 drives the grounding line over this stabilizing region and into the deeper beds upstream. From then on, volume above flotation losses increase to 150 km³a⁻¹. Applying uniform melt rates of 16 m a⁻¹pushes the glacier over this region earlier, but the mean rate of volume loss still increases, from 110 to 170 km³a⁻¹in the 0⁰⁰/U16 experiment.

3.2.2 Filchner–Ronne Ice Shelf

The Weddell Sea ice streams feeding the Filchner–Ronne Ice Shelf speed up and lose mass in response to the FESOM melt rates from the mid 21st century, thinning faster in the 22nd century and faster still in the 23rd. The H/A/0/F experiment sees 22nd-century grounding line retreat in the Evans, Möller, Institute and Foundation Ice Streams, with the remaining ice streams starting to retreat only after 2200. Fig- ure 9 shows the corresponding loss of volume above flotation, increasing from 20 km³a⁻¹ (5 mm SLE per century) over the 21st century through 60 km³a⁻¹ (15 mm SLE per century) over the 22nd – a rate comparable to Pine Island Glacier in the 21st century. Although the melt rate anomaly is held at its late-22nd century values over the 23rd century, the rate of retreat continues to increase, reaching 210 km³a⁻¹ (50 mm SLE per century).

All of the ice streams respond immediately to the higher melt rates imposed in the uniform melt rate experiments.

Grounding line retreat is most apparent in the Evans, Möller and Institute Ice Streams, but even the narrow Rutford Ice Stream and Carlson Inlet see grounding line migration over 50 km or more by 2200. Rates of retreat in the Evans, Möller and Institute Ice streams are comparable to rates of retreat, given the same forcing, in Thwaites and Pine Island Glacier, while volume loss in the two regions is similar over the 21st century: 220 vs. 180 km³a⁻¹ and somewhat lower in the 22nd century: 340 vs. 500 km³a⁻¹.

The Möller and Institute Ice Streams exhibit a century- long period of accelerated retreat in all of our simulations, though its onset varies considerably. The H/A/0/F experiment shows some retreat from 2100–2200 in response to the FESOM projection of increased melt rates under the Filchner–Ronne ice shelf, but from 2200 onward the grounding lines of the two ice streams merge and then migrate across the bulk of the deep-bedded Robin Sub-glacial Basin in a single century to reach down-sloping beds by 2300 CE (Fig. 7). At the same time, the Bungenstock ice rise is isolated and then floats. Compared to that, retreat begins right away in the 0/U16 experiment, reaches the down-sloping bed by 2100 CE, and retreats little more after that, while the 0/U8 calculation produces a similar period of retreat starting in 2050. Figure 9 shows the change in volume above flotation for each of these simulations: all three sustain a max- imum rate of volume loss, approximately 210 km³a⁻¹, for the hundred-year period corresponding to this retreat. These streams, it would seem, are close to marine ice sheet instability as seen in Wright et al. (2014) and can be forced into unstable retreat by physically plausible (FESOM) melt rates.

We compute faster retreat here than Wright et al. (2014) simply because the melt rates are greater.

(13)

3.2.3 Ross Ice Shelf and Siple Coast

The Siple coast ice streams feeding the Ross Ice Shelf lose up to 50 mm SLE by 2100 and 150 mm SLE by 2200 in response to increased melt rates, but respond less to melt rates derived from FESOM. Figure 9 shows the mass loss trend for the control, H/A/0/F, 0/U16 and U/08 experiments, and grounding lines are shown for the first three in Fig. 7. The MacAyeal, Bindschadler, Mercer and Whillans ice streams all flow over retrograde beds and exhibit∼1 km a⁻¹grounding line retreat in the 0/U16 and U/08 experiments and lose volume above flotation at a rate of around ∼200 km³a⁻¹. Although the Mercer and Whillans ice streams do retreat in response to the H/A/0/F experiment, they retreat in a similar fashion in the control experiment. Both the H/A/0/F and control grounding lines sweep across an area which is lightly grounded in the model’s initial state but just floating in (for example) Bedmap2 (Fretwell et al., 2013). The MacAyeal and Bindschadler ice streams do start to retreat around 2150, when the FESOM melt rates grow to 10 m a⁻¹in the region of Roosevelt Island; that retreat is sustained, and accelerates, until the end of the experiment in 2300. The accompanying loss of volume above flotation amounts to about 60 km³a⁻¹. There is little retreat apparent in the inactive Kamb ice stream in any of the experiments.

3.2.4 Marie-Byrd Land

Marie-Byrd Land shows little sign of retreat in the H/A/0/F experiments, amounting to around 0.5×10³km³ (1.5 mm SLE), despite the elevated melt rates imposed from 2100 onward. Even when uniform 16 m a⁻¹ melt rates are imposed across the ice shelves, the grounding line retreats by only a few kilometers over the 200 years, with an accompanying loss of volume above flotation of 1.6×10³km³(5 mm SLE).

It is apparent in Fig. 7 that the present-day grounding line in this region runs perpendicular to a down-sloping bed, with much of the bed above sea level, so there is no possibility of marine ice sheet instability. As a result, we might not expect to see century-scale retreat. On the other hand, all of the glaciers in this region are rather narrow, which leaves the possibility that they are under-resolved even at1xmin=625 m, and their beds are inevitably under-resolved by the sparse bedrock data in this region (Fretwell et al., 2013).

4 Conclusions

Our most extreme simulation of widespread dynamic thinning in West Antarctica’s fast-flowing ice streams results in 200 mm of eustatic sea level rise by 2100 and 475 mm by 2200. Pine Island and Thwaites Glaciers see their grounding lines retreat by hundreds of kilometers, as do the Möller, In- stitute, Evans, MacAyeal, Bindschadler, Whillans and Mer- cer ice streams and to a lesser extent Carlson Inlet and the Rutford Ice Stream. All of these ice streams flow along beds

that deepen inland, and so can be subject to marine ice sheet instability. Some of the ice streams appear to be on the edge of critical change; for example Pine Island Glacier and the Möller, Institute and Evans ice streams remain close to their present-day configurations unless melt rates are increased.

Our model of Thwaites Glacier, on the other hand, depends strongly on its initial state: either it remains steady for up to 200 years after its ice shelf has all but disappeared, or it retreats rapidly, raising sea level by at least 25 mm each century, even if its ice shelf remains in place.

Wholesale retreat occurs only if enhanced oceanic melt rates are imposed across all the ice shelves, but neither the FESOM nor the BRIOS ocean circulation models project substantial warming beneath the Filchner–Ronne or Ross ice shelves until after 2050. Simulations based upon these more realistic projections also result in significant dynamic losses in the Amundsen Sea Embayment: up to 50 mm SLE by 2100 and 150 mm SLE by 2200 provided that Thwaites Glacier retreats. On the other hand, there is little retreat in the Filchner–

Ronne or Siple Coast ice streams until after 2100, and only around 30 mm SLE of thinning by 2200. The Möller and In- stitute ice streams do exhibit accelerated retreat immediately after 2200, increasing their contribution from 20 mm SLE in 2200 to 75 mm SLE by 2300.

Both the RACMO2 and LMDZ4 atmosphere models project increasing snowfall given the A1B emissions scenario, which partly offsets any dynamic thinning. We found that the effect of increasing accumulation could be separated from the effect of increasing melt rates, with the ice sheet model responding to the two perturbations independently.

Up to 20 mm SLE of extra accumulation by 2100 and as much as 75 mm by 2200 is dispersed across West Antarctica, sufficient to balance the FESOM- or BRIOS-driven contribution of Pine Island Glacier and the Möller and Institute ice streams – but not Thwaites Glacier – over the same period. The E1 projections do not show increased accumulation, and in those cases the dynamic thinning, which varies much less between emissions scenarios, is offset by no more than 20 mm SLE by 2200.