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

. We use the BISICLES adaptive mesh ice sheet model to carry out one, two, and three century simulations of the fast-ﬂowing 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, evolves according to variation in meteoric ice accumulation rates and oceanic ice shelf melt rates. Future changes in accumulation and melt rates from no change, by and ocean models driven the and uniform rate that


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 underestimation 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.

Introduction
The present-day West Antarctic Ice Sheet (WAIS) is experiencing 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. Physically 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 Published by Copernicus Publications on behalf of the European Geosciences Union.
migration. Progress has been made in both areas (Goldberg et al., 2012a, b;Pattyn et al., 2013), but the computational expense 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 . Two emission scenarios are included: SRES E1, a mitigation scenario in which emissions are stabilized by 2050 at 500 ppm CO 2 , 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 (Ligtenberg 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 providing 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 complement 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).

Model equations
BISICLES employs a vertically integrated ice flow model based on Schoof and Hindmarsh (2010) which includes longitudinal 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 elevation b and ice thickness h the upper surface elevation s is in which ρ i and ρ w are the densities of ice and ocean water. The ice thickness h and horizontal velocity u satisfy a twodimensional mass transport equation and two dimensional stress-balance equation together with lateral boundary conditions. The terms on the right hand side of Eq.
(2), a and M, 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, and I is the identity tensor. The vertically integrated effective viscosity φhμ is computed by evaluating the integral numerically, with the ice sheet sub-divided into 10 layers, narrowing progressively from 0.16h near the surface to 0.03h near the base. The vertically varying effective viscosity µ(x, y, z) includes a contribution from vertical shear and satisfies where the flow rate exponent n = 3, φ is a stiffening factor (or, equivalently, φ −n is an enhancement factor), and A(T ) depends on the ice temperature T through the Arrhenius law described by Hooke (1981), The Cryosphere, 9, 1579-1600, 2015 www.the-cryosphere.net/9/1579/2015/ where A 0 = 0.093 Pa −3 a −1 , Q/R = 9.48 × 10 3 K, f = 0.53 K k , k = 1.17 and T r = 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 temperature T and the form of A(T ), macroscopic damage, and fabric formation. Finally, the basal traction is determined by a viscous law: 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 excessive grounding line retreat: in previous work on Pine Island Glacier non-linear laws with m < 1 have led to slower rates of retreat (Joughin et al., 2010;Favier et al., 2014). We hold the fields C 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 of C and φ we use (see Sect. 2.4.1) already feature slippery beds and weak shear margins hundreds of kilometers upstream from the grounding line.

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 Embayment, 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. If n is normal to a boundary and t is parallel to it, u · n = 0, t · ∇u · n = 0, ∇h · n = 0. 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 5 Pa m −1 a so that at the interface between V and Q , while at the calving front (which is fixed), we impose the usual conditions on the normal and transverse stress: These boundary conditions, and indeed, Eq. (11) alone for a problem whose entire boundary is a fixed calving front, are sufficient provided that h(x, y, t = 0) is given and that the basal friction coefficient C(x, y) is non-zero in at least part of the ice sheet.

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  www.the-cryosphere.net/9/1579/2015/ 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 resolution x , with 0 ≤ ≤ L and 2 x +1 = x . While restrictive in some senses -all model domains must be rectangular, for example -these meshes have a signal advantage: it is straightforward to generate 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 calculated 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 subkilometer 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 ≤ x ≤ 4000 m, than in the Ronne-Filchner and Ross ice shelf catchments; there we find that a mesh with 625 m ≤ x ≤ 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) condition t < 0.25 x l |u| 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.

Model data requirements
Time-dependent simulations require initial ice thickness data h 0 (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 field T (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 measuring 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 highresolution 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;, 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 fields C(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 . 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 −1 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 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 rates a 0 and M 0 chosen to hold the ice sheet close to equilibrium. Determination of a 0 and M 0especially M 0 -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 relaxation period, we compute a 0 and M 0 from ∇ · [uh], with M 0 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 of a 0 will tend to counter the thinning that is already evident in the ASE, we carry out some additional experiments where a 0 is replaced by an accumulation pattern derived from the RACMO atmosphere model. . Melt rate anomalies integrated over each region. In contrast with the atmosphere models, the ocean models provide similarly 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 Circumpolar Deep Water temperatures.

Initial state
Ice Shelf, and through the Ross Ice Shelf are apparent in both the flow field and the basal traction field C. The flow field itself 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 −1 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 a 0 and M 0 ) is between −5 and 5 m a −1 except at calving fronts. Integrating this thickening rate leads to an annual loss of volume above flotation of 3 km 3 a −1 in the ASE, 7 km 3 a −1 in the Filchner-Ronne ice shelf basin and 7 km 3 a −1 in the Ross ice shelf basin. The synthetic accumulation a 0 used to obtain this thickening rate does include some unrealistic large-amplitude shortwavelength 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, with a 0 between −5 and 5 m a −1 : in particular there is no region of ∼ 100m a −1 ablation needed to counter the flux arriving at the Pine Island Glacier grounding line.

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.

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 a 0 (x, y) and M 0 (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 The Cryosphere, 9, 1579-1600, 2015 www.the-cryosphere.net/9/1579/2015/ 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 T (t) by averaging the projected ocean temperature over volume bounded laterally between the contemporary 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 M(t) = 16 T (t) m a −1 K −1 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 concentrating 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 threefold 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.  Figure 5. Net change in volume above flotation over the course of the combined anomaly experiments. Only the Amundsen Sea Embayment experiences a net loss ( V ) 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 3 km 3 loss by 2100 and 40 × 10 3 km 3 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.

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 dramatic 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 field a 0 (x, y), which we constructed to hold Thwaites Glacier in steady state during the control experiment, includes a spot of unrealistic -5 m a −1 -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 field a 0 (x, y) with the 5 m a −1 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 0 (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 0 (x, y, t).
As none of the melt rate anomalies in the ASE exceed 10 m a −1 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 −1 applied across all floating ice from 1980 onward. 0/U16 used the synthetic accumulation field a 0 (x, y), while 0 /U16 used the HadCM3/E1/RACMO2 1990-2000 mean accumulation field a 0 (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 −1 across the entire extent of floating ice starting from 1980.

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

Combined anomaly experiments
Only the Amundsen Sea Embayment experiences a net loss of volume above flotation ( V ) 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 3 km 3 by 2100 and 3-23 × 10 3 km 3 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 3 km 3 by 2100 and 40 × 10 3 km 3 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. Figure 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 3 km 3 greater loss (− V ) 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 3 km 3 , 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 3 km 3 during the 21st century, and around 20 × 10 3 km 3 by the end of the 22nd century (assuming 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 3 km 3 loss by 2200 and H/A/R/B less than 5 × 10 3 km 3 . 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, V d . For a given region , V d is the difference between the net change in volume above flotation and the cumulative, integrated accumulation anomaly: 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 in V d ≈ 10 × 10 3 km 3 (25 mm SLE) across the whole of West Antarctica by 2100 and around 30 × 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 −1 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 3 km 3 (100 mm SLE) by 2200. Overall, the net V (t) for a simulation with accumulation anomaly a(t) and melt rate anomaly M(t) can be estimated rather precisely from the result, V (t), of a simulation with the same ocean anomaly and a different (or no) accumulation anomaly a (t): 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 accumulation causes low amplitude but widely distributed thickening.

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.
Pine Island Glacier and the ice streams feeding the Crosson and Dotson Ice Shelves experience ∼ 1 km a −1 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 affected by initial conditions, with some simulations showing little retreat and others shedding between 100 and 210 km 3 a −1 volume above flotation over the 21st and 22nd centuries. In calculations with the synthetic mass balance a 0 : 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 −1 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 3 a −1 ) than in the control (∼ 100 km 3 a −1 ).
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 3 a −1 (0.2 mm a −1 SLE) between 2000 and 2100, while the grounding line retreats at a rate ∼ 1 km a −1 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 −1 , and the rate of loss of volume above flotation to 320 km 3 a −1 (0.9 mm a −1 SLE). A similar calculation of accelerating mass loss, with losses of less than 0.25 mm a −1 SLE during the 21st century and up to 1 mm a −1 SLE thereafter was reported by Joughin et al. (2014). 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 3 a −1 over the 21st century, comparable to rates computed in other modeling studies Favier et al., 2014), and present-day observations . 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 3 a −1 . Applying uniform melt rates of 16 m a −1 pushes the glacier over this region earlier, but the mean rate of volume loss still increases, from 110 to 170 km 3 a −1 in the 0 /U16 experiment.

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. Figure 9 shows the corresponding loss of volume above flotation, increasing from 20 km 3 a −1 (5 mm SLE per century) over the 21st century through 60 km 3 a −1 (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 3 a −1 (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 3 a −1 and somewhat lower in the 22nd century: 340 vs. 500 km 3 a −1 .
The Möller and Institute Ice Streams exhibit a centurylong 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 maximum rate of volume loss, approximately 210 km 3 a −1 , 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. , 9, 1579-1600, 2015 www.the-cryosphere.net/9/1579/2015/

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 −1 grounding line retreat in the 0/U16 and U/08 experiments and lose volume above flotation at a rate of around ∼ 200 km 3 a −1 .
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 −1 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 3 a −1 .
There is little retreat apparent in the inactive Kamb ice stream in any of the experiments.

Marie-Byrd Land
Marie-Byrd Land shows little sign of retreat in the H/A/0/F experiments, amounting to around 0.5 × 10 3 km 3 (1.5 mm SLE), despite the elevated melt rates imposed from 2100 onward. Even when uniform 16 m a −1 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 3 km 3 (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 at x min = 625 m, and their beds are inevitably under-resolved by the sparse bedrock data in this region (Fretwell et al., 2013).

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, Institute, Evans, MacAyeal, Bindschadler, Whillans and Mercer 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 Institute 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.

Appendix A: Mesh resolution
Both solution accuracy and computational cost increase as the computational mesh is refined. We need to show that our chosen meshes, with x min = 250 m for the ASE and x min = 625 m elsewhere, are adequate. To this end, we consider just the H/A/0/F simulations, as the dynamic response is similar for every case.
Both grounding line migration and overall mass loss in the ASE model are rather sensitive to mesh resolution close to the grounding line. Figure A1 plots the change in volume above flotation ( V ) over time for meshes with finest resolutions varying from x min = 4000 m to x min = 250 m, while Fig. A1 shows grounding line positions around Smith glacier at the start and end of the simulations, together with the final mesh. At the coarsest resolution there is essentially no loss of volume and the grounding line does not retreat at all, while at finer resolutions there is progressively greater volume loss and grounding line retreat, just as in previous simulations of Pine Island Glacier (Favier et al., 2014;Cornford et al., 2013).
Provided that the finest resolution x min ≤ 1000 m, the volume change V ( x min , t) converges with x min , but we only obtain the theoretical rate of convergence -O( x min ) -when x min ≤ 500 m. If the rate of convergence is O( x min ), we can calculate an estimate of the volume loss for any resolution given the results from two coarser resolutions, through Richardson extrapolation. The error estimate for a given resolution is the estimated volume loss as x min → 0 is and the estimated volume loss for a finer resolution is Figure A1 plots these estimated volume losses alongside the simulated losses: V ( 1 2 x min , t) is a good approximation to V ( 1 2 x min , t) once x min = 500 m. For the meshes with x min = 250 we estimate that volume losses over 200 years fall short through under-resolution, by e( x min , t) = 2.5 × 10 3 km 3 out of 20 × 10 3 km 3 , that is, 7 mm sea level equivalent (SLE) out of 50 mm SLE.
The RISFRIS model is less sensitive to mesh resolution than the ASE model, and we were able to carry out calculations using coarser meshes with x min = 625 m. Figure A2 shows that much of the loss in volume above flotation seen in a simulation with x min = 312.5 m ( V = 32 mm SLE) is also seen in a simulation with x min = 5000 m ( V = 21 mm SLE). At the same time, while the grounding line does retreat more as the resolution is increased, retreat takes place for all resolutions in the Möller, Institute and Foundation ice streams, though only for x min < 2500 m in Evans Glacier. Nonetheless, the simulation with x min = 1250 m loses as much additional volume (4 mm SLE) compared to the x min = 2500 m case as calculations with x min = 2500 m lose on top of the x min = 5000 m case. Convergence with x min is apparent once x min < 2500 m and at the chosen resolution, x min = 625 m, volume loss is underestimated by e( x min , t) = 0.8×10 3 km 3 out of 12×10 3 km 3 (2 mm SLE out of 30 mm SLE).
In summary, under-resolution results in the loss of volume above flotation being underestimated by around 10 % for our chosen meshes, which have

B1 Basal friction and stiffening coefficients
We estimate the basal friction and stiffening coefficients by solving an inverse problem similar to those of MacAyeal (1993), Joughin et al. (2009) and Morlighem et al. (2010). Broadly speaking, we choose smooth fields C(x, y) and φ(x, y) that minimize the mismatch between modeled and observed speeds. A nonlinear conjugate gradient method was employed to seek a minimum of the objective function composed from a misfit function and a Tikhonov penalty function www.the-cryosphere.net/9/1579/2015/ The Cryosphere, 9, 1579-1600, 2015 The coefficient α 2 u (x, y) is related to error estimates for the observed velocity and we set it to 1 where velocity data are available and 0 elsewhere.
Ideally, the penalty function J p would not be required (so that α 2 C , α 2 φ = 0) because it is equivalent to a claim to some prior knowledge of φ and C -specifically, their probability distributions are p(C) ∝ exp(−α −2 C |∇C| 2 ) and p(φ) ∝ exp(−α −2 φ |∇φ| 2 ) -and we clearly do not have such knowledge. In practice, though, we require α 2 C , α 2 φ > 0 for two reasons. First, J does not have a unique minimum with respect to both φ and C: in other words, the inverse problem would be under-determined, because we have one field of data and two fields of unknowns. Second, even if we were only seeking (say) C, the inverse problem would be ill-conditioned, that is, sensitive to small changes in u o . We follow Hansen (1994) and choose α 2 C , α 2 φ such that lower values lead to faster growth in J p than reduction in J m and larger values lead to the converse.
To ensure that C and φ are positive definite, we express them as and minimize Eq. (B1) with respect to p and q. C 0 and φ 0 are initial guesses for the basal traction and stiffening factor, and here we set if α 2 u > 0 10 5 otherwise (B5) and φ 0 = 1. The nonlinear conjugate gradient method requires expressions for the Gâteaux derivatives d q J (q, p; q ) and d p J (q, p; p ). Here we have used similar notation to Arthern and Gudmundsson (2010), for example: We neglect the non-linearity ofμ in the adjoint of the viscous tensor operator in Eq.
(3) -an omission discussed by e.g. Goldberg and Sergienko (2010) and Morlighem et al. (2013) -and approximate the Gâteaux derivatives d q J m (q, p; q ) and d p J m (q, p; p ) by and where the vector field λ is the solution to a linear boundary value problem, with and on the domain boundary λ · n = 0, t · ∇λ · n = 0.
The penalty function (Eq. B3) is easily differentiated to give and The inverse problem requires accurate speed data consistent with the thickness and bedrock data and having nearcomplete spatial coverage. For the RISFRIS and MBL domains we use InSAR observations made between 2007 to 2009 . We could also use this data for the ASE domain, but as that area has accelerated over the last decade we choose measurements made in 1996 (Joughin et al., 2009). The essential difference between the two is acceleration in the region of Pine Island Glacier's ice shelf, which sped up by around 1 km a −1 . Both data sets have regions of missing data: in the region of the grounding line in Joughin et al. (2009) and in isolated spots elsewhere, and there we rely on the Tikhonov regularization terms in C and φ to provide continuous fields.
The misfit function J m defined in Eq. (B2) is reduced rapidly in the first few conjugate gradient iterations in each of the regional models. Figure B1 shows the relative size of J m over 32 iterations for all three models, and in each case, the relative size of J m is reduced by an order of magnitude. Figure B2 shows maps of the mismatch |u o | − |u| at iteration 0 and iteration 24. Table B1 lists the mean values of the pointwise mismatch |u o | − |u| and the observations |u o | for each of the domains and their fast-flowing portion, where |u o | > 300m a −1 . In the ASE region, the ratio of the mismatch to the observations is reduced from 0.50 to 0.16, in the RISFRIS domain from 0.31 to 0.10, and in the MBL domain from 0.62 to 0.16. In the fast-flowing regions (which our choice of J m tends to favor) the ratio drops from 0.37 to 0.04 in the ASE domain, from 0.25 to 0.03 in the RISFRIS domain, and from 0.45 to 0.04 in the MBL domain.
An adjustment in the region of Pine Island Glacier's grounding line was required to prevent sustained thickening of the order of 100 m a −1 . A similar tendency is seen in other models of Pine Island Glacier, and is dealt with elsewhere 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  The Cryosphere, 9, 1579-1600, 2015 www.the-cryosphere.net/9/1579/2015/  Figure B1. Progress of the nonlinear conjugate gradient method. After around 10 iterations the speed misfit J m (Eq. B2) is reduced by an order of magnitude in the ASE, MBL and RISFRIS optimization problems. 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: where R = 60 km, (x P , y P ) is located in the center of the glacier trunk close to the grounding line and α = 0.2. This results in around 500 m a −1 faster flow in the center of Pine Island Glacier's ice shelf, so that the model lies somewhere between the 1996 observed speeds used in the inverse problem, and the 2007-2008 PALSAR observations .

B2 Relaxation, initial accumulation and initial melt rate
Having matched the observed speeds, we find shortwavelength and large amplitude variation in the ice thickening rate (Fig. B3). In the ice shelves, we could attribute a large part of these fluctuations to essentially unknown melt rates, but in the ice streams we assume them to be artifacts of interpolation and other sources of error in the ice sheet geometry Morlighem et al., 2011), or mismatch between the time at which the geometry and velocity were observed. So before carrying out all the targeted experiments we run (relax) the model for a period with a present-day forcing to bring it closer to a steady-state. The relaxation is carried out in two stages. First we set the SMB to the 1990-1999 mean from a high-resolution atmosphere model (RACMO2, with HadCM3 boundary conditions and the E1 emissions), and the sub-shelf melt rate is chosen to keep the ice shelf in steady state. After 50 years, we compute an accumulation rate a 0 required to keep the grounded ice close to steady state, and a melt rate M 0 that will keep the ice shelf close to steady state. We then run the model for 50 years starting again from the original state. The resulting ice sheet is closer to, but not at, equilibrium, and we carry out all projections starting from this state, with SMB and melt rates computed by adding the perturbations described in Sect. 2.5.1 to a 0 and M 0 . Although the melt rate M 0 at the end of the relaxation could be estimated directly from the flux divergence, we will need to adjust it as the grounding line moves. Flowline calculations indicate that elevated melt rates close to the grounding line can result in a dynamic response quite different from the response to elevated melt some distance downstream (Walker et al., 2008;Gagliardini et al., 2010). Observations and ocean models tend to show that melt rates do decay downstream of the grounding line (Jenkins et al., 2006;Payne et al., 2007;Rignot et al., 2013;Dutrieux et al., 2013Dutrieux et al., , 2014, and we expect peak melt rates to move with the grounding line, partly because pressure melting point is higher there, and partly because the steeper underside of the ice shelf leads to more entrainment of heat and salt. We construct a scheme that allows higher melt rates to follow the grounding by decomposing M 0 into grounding line localized and ambient components, as in Gong et al. (2014) and Wright et al. (2014). We set M 0 (x, y, t) = M GL (x, y)ξ(x, y, t) + M A (x, y)(1 − ξ(x, y, t)), where ξ = 1 at the grounding line, ξ = 0 in the ocean, and 1 Figure B2. Difference between observed and model speeds. |u o | − |u| at the start (a) and end (b) of the inverse problem. The initial guess for C and φ results in a flow field that is up to 1 km a −1 too slow in the fast ice streams. After 24 conjugate gradient iterations, the model and observed speeds rarely differ by more than 100 m a −1 .  Figure B3. Thickening rates before relaxation. Spatial and temporal mismatch between the observed velocity and thickness data leads to an ice sheet with large amplitude, short wavelength fluctuations in ∇ · [uh] and hence ∂h ∂t . Pine Island Glacier is particularly affected, with thickening rates ∼ 100 m a −1 far in excess of observed values, and having the wrong sign. the initial grounding line respectively, and then extrapolating over the entire domain. To compute the ambient component, we construct a parabolic equation (B18) Figure B4 shows how M 0 (x, y, t) varies in Pine Island Glacier's ice shelf as the grounding line retreats. Thwaites Glacier thinned and the grounding line began to retreat given the accumulation and melt rates described above, without any additional forcing. This was prevented in the combined anomaly simulations by adding a highly localized (10 km radius) region of 5 m a −1 extra accumulation to a 0 . Recent observations and modeling indicate that the unstable retreat of Thwaites Glacier may have begun early in the 21st century Joughin et al., 2014), and so we also carry out simulations with an accumulation field a 0 that does not include the localized addition.
Our intention in utilizing a synthetic accumulation was to separate projections of future change due to future ocean warming from changes that are already evident in contemporary observations. At least in the ASE, the use of a synthetic accumulation serves to counter observed present-day dynamic thinning, which may lead to future marine ice sheet instability without any additional forcing, or may combine with any future forcing in a non-linear fashion. To address this issue, we also report results computed in that region using the HadCM3/E1/RACMO2 2000-2009 mean accumulation data (which we will label a 0 ) in place of a 0 , in which case the ASE basin loses volume above flotation at a rate of 50 km 3 a −1 in the absence of any anomaly data.