Future Antarctic bed topography and its implications for ice sheet dynamics

. The Antarctic bedrock is evolving as the solid Earth responds to the past and ongoing evolution of the ice sheet. A recently improved ice loading history suggests that the Antarctic Ice Sheet (AIS) has generally been losing its mass since the Last Glacial Maximum. In a sustained warm-ing climate, the AIS is predicted to retreat at a greater pace, primarily via melting beneath the ice shelves. We employ the glacial isostatic adjustment (GIA) capability of the Ice Sheet System Model (ISSM) to combine these past and future ice loadings and provide the new solid Earth computations for the AIS. We ﬁnd that past loading is relatively less important than future loading for the evolution of the future bed topography. Our computations predict that the West Antarctic Ice Sheet (WAIS) may uplift by a few meters and a few tens of meters at years AD 2100 and 2500, respectively, and that the East Antarctic Ice Sheet is likely to remain unchanged or subside minimally except around the Amery Ice Shelf.


Introduction
Projecting the evolution of the Antarctic Ice Sheet (AIS) into the next few centuries relies on simulating a complex and non-linear coupled Earth system. A recent survey of experts by Bamber and Aspinall (2013) reveals that projections for AIS contribution to the rate of sea level rise at the year AD 2100 are generally rather moderate (∼ 1.7 mm yr −1 ) and that the upper end of the spectrum of projections would be about 7 times this value, mainly owing to intensification of dynamics of the West Antarctic Ice Sheet (WAIS). However, projections beyond AD 2100 are much more uncertain (Bindschadler et al., 2013) and are mainly limited by the poor knowledge of the physics involved in the grounding line (GL) migration and ice shelf melting and calving (e.g., Vaughan and Arthern, 2007;Walker et al., 2013). Furthermore, there is strong evidence that over the past four million years, during times of increased global atmospheric temperatures by 2-3 • C, the marine WAIS collapsed, and possibly some portions of the larger East Antarctic Ice Sheet (EAIS) as well (Naish et al., 2009;Raymo et al., 2011;Cook et al., 2013).
If we are to improve our abilities to assess the risk of the catastrophic consequences of partial collapse of AIS marine-based ice currently locked to the Antarctic continent, steps must be taken to assess fully the role of solid Earth deformation over tens to thousands of years, during which time gravitational viscoelastic flow of the underlying mantle may act to change the stability of the AIS. Past assessments have been that isostatic uplift following the ice sheet retreat promoted stability of the WAIS near the end of the last deglaciation during the mid-Holocene (Thomas, 1976). More recently, the promotion of stability by glacial isostatic adjustment (GIA) has been shown with increasingly sophisticated computations (Gomez et al., 2010(Gomez et al., , 2013. It is also now recognized that past AIS recession and ice flow direction are plausibly explained by strong interaction of ice loading with the solid Earth (Siegert et al., 2013). It is this ice-sheet/solid-Earth (IS/SE) interaction that we now explore in this paper for the AIS as a whole, using the vertical surface motions of bedrock GIA (Ivins and James, 1999) and the 3-D thermomechanical ice flow (Pattyn, 2003) capability of the Ice Sheet System Model (ISSM) (Larour et al., 2012b).

IS/SE feedback: GL migration
Perhaps the most important IS/SE feedback is associated with the GL migration (e.g., Lingle and Clark, 1985). For equilibrium sea level, any change in vertical bedrock elevation also changes the local sea depth. This perturbs buoyant forces in the regional ocean water and may promote the migration of the GL. Uplift of the seabed occurs, for example, in response to thinning of the inland ice sheet, and it causes local sea depth to decrease. If the GL is in initial equilibrium but sea depth decreases due to bed uplift, the GL tends to advance towards the continental shelf (e.g., Weertman, 1974;Schoof, 2007). The conceptual illustration of this negative feedback is shown in Fig. 1a. Lingle and Clark (1985) explored the effect of GIA-related seabed uplift on GL migration for Ice Stream E, now known as the MacAyeal Ice Stream, in the WAIS. The modeled GIA uplift, caused by thinning of the ice stream catchment area, delays the onset of GL retreat, thus reducing the rate of retreat during Holocene sea level rise. Notably, it was argued that the regional advance of Ross Sea GL may have occurred over the past three millennia. The gravitational attraction effects on local sea level developed later by Gomez et al. (2010) act to amplify this negative feedback during ice sheet retreat, since the diminished mass behind the GL has less mutual gravitational attraction with adjacent sea water, thus causing local sea level to drop.
The pace and magnitude of GL migration are also dictated by the bedrock slope. If sea depth decreases (due to bed uplift and sea level drop associated with the thinning of inland ice) at the initially stable GL on a reverse bed slope, for example, the GL tends to advance further on a relatively flat bed. Therefore the GLs associated with the Ross and Ronne ice shelves in the WAIS having relatively flat reverse bed slopes are sensitive to this effect (Conway et al., 1999). Capturing the dynamics of such GL sensitivity demands proper understanding of interactions between the ice, the ocean and the solid Earth, and is indeed the key to successful modeling of the WAIS (e.g., Vaughan and Arthern, 2007;Katz and Worster, 2010).

Additional IS/SE feedbacks
There are other important feedback mechanisms that the solid Earth deformation provides to ice sheets. The GIA uplift can be important in providing basal resistance to ice flow and buttressing the ice sheet by raising bedrock pinning points (e.g., Favier et al., 2012;Siegert et al., 2013). This concept is captured in Fig. 1b. The GIA-induced changes in surface elevation and regional slope of the ice sheet may affect the gravitational driving stress, as well as some processes at the ice-atmosphere interface (e.g., surface mass balance). These perturb the momentum balance and affect the englacial velocity field (e.g., Cuffey and Paterson, 2010;Winkelmann et al., 2012), which in turn may impact the ice thermodynamics via changes in strain heating. Spatially varying bed uplift also affects the hydraulic potential field and hence the subglacial hydrology and the sliding rate of the ice sheet. There is additionally the complication of bulge migration, a broad-scale phenomenon involving bending of the elastic lithosphere and mantle lateral flow. Due to the lateral motion of this topographic bulge, local crustal motions (and slopes) may change sign during GIA (e.g., Fjeldskaar, 1994).
These mechanisms are extremely difficult to isolate and quantify, and it is therefore not obvious whether (and under what circumstances) each of these acts to accelerate or inhibit the ice flow. As long as the thermomechanical ice sheet model and other companion models (e.g., surface mass balance model and the hydrological model) are dynamically coupled to a comprehensive solid Earth model, however, most of these feedbacks are intrinsically taken into account.

New solid Earth computations
For large timescale (millennia) simulations, most of the ice sheet models (e.g., Pollard and DeConto, 2009;Hughes et al., 2011) capture the solid Earth physics of varying degrees of complexity (cf. Le Meur and Huybrechts, 1996). With the exception of the recent work of Gomez et al. (2013), none of these studies provides the explicit assessment of effects of GIA uplift on several aspects of ice sheet dynamics (e.g., GL migration and gravitational driving stress). Le Meur and Huybrechts (2001) explicitly pointed out the need for the more complete coupling that could be found in the wavelength-dependent relaxation spectra of viscoelastic solid Earth models. In this paper, we quantify two distinct IS/SE feedback mechanisms applied to the AIS on centennial timescales using multiple wavelength-dependent decay spectra. Assuming the equilibrium sea level, we first evaluate how the future bed uplift (governed by the past and future evolution of AIS; cf. Sect. 2) controls the GL migration. Second, we assess the role of bed uplift in the gravitational driving For equilibrium sea level, the GIA uplift due to thinning of the inland ice promotes the GL migration towards the continental shelf. The hydrostatic equilibrium requires that s 1 /b 1 = s 2 /b 2 = −(ρ w −ρ)/ρ, where s i and b i are the surface and bedrock elevations (subscript i = 1, 2 represents the initial and final configurations, respectively), and ρ w and ρ are the water and ice densities. (b) Pinning point, raised by GIA uplift due to thinning of the inland ice, provides basal resistance and buttressing to the ice sheet. Red arrows depict the velocity profiles. Dashed lines in the final configurations (second column) represent the initial geometry. Note that the destabilizing effects of GIA during the periods of inland ice thickening that causes the bedrock to subside can be conceptualized by reverting the direction of mid-arrows.
stress. The overall influence of solid Earth deformation (via changes in GL and driving stress) on the ice surface velocities is also quantified. These assessments, based on reliable models and data, help us to understand whether the GIA effects are significant in controlling the future evolution of AIS on centennial timescales.
The paper layout is as follows. In Sect. 2, we introduce the solid Earth model employed in this research, discuss how the change in ice thickness is translated into the ice load height, provide a detailed account of model tuning procedure, and describe all of the required model input data. In Sect. 3, we present modeling results of the future bedrock topography, and assess the relative importance of the past and future ice load changes. In Sect. 4, we quantify the influence of predicted change in Antarctic bed topography on several aspects of ice sheet dynamics. Finally, in Sect. 5, we summarize key conclusions of broad interest. Ivins et al. (2013) presented a much improved ice loading history for the AIS since the Last Glacial Maximum (LGM). As for the future, the recently concluded SeaRISE (Sea-level Response to Ice Sheet Evolution) project (Bindschadler et al., 2013; provided quantitative projections of the evolution of AIS under ongoing climate warming. We employ the new GIA capability (Ivins and James, 1999) of ISSM (Larour et al., 2012b), hereinafter referred to as the ISSM/GIA model, to combine these data of past and future ice loadings and calculate the change in bedrock topography over the same timescale projections as in the SeaRISE studies. Using appropriate analytical and numerical models, we then evaluate the stabilizing or destabilizing influence of predicted changes in bed topography on the ice sheet dynamics.

Model and data
While the SeaRISE experiments employed state-of-the-art numerical treatments of ice flow, it should be noted that the majority of these models were not coupled to a comprehensive solid Earth model. Furthermore, the SeaRISE experiments do not capture the paleo-evolution of the AIS since the LGM, thus limiting the possibility for the participating ice sheet models (with GIA capability) to perform similar analysis presented in this study.

The ISSM/GIA model
ISSM is a continental-scale, high-resolution, multi-model simulation code developed for understanding the dynamic behavior of large ice sheets (Larour et al., 2012b). This open source finite-element software is capable of simulating ice-flow mechanics of varying degrees of complexity , performing sensitivity analysis (e.g., Larour et al., 2012a), inverting unknown field parameters (e.g., Morlighem et al., 2010), and assessing mechanics of rift propagation and eventual collapse of ice shelves (e.g., Borstad et al., 2012). The semi-analytical GIA solution of 572 S. Adhikari et al.: Antarctic ice dynamics and isostatically evolving bed Ivins and James (1999) is one of several new features being actively implemented in ISSM. Here we briefly summarize the key elements of the model.
We assume that the ice sheet rests on top of the solid nonrigid Earth, which is considered to be a simple two-layered incompressible continuum with an upper elastic lithosphere floating on the viscoelastic (Maxwell material) mantle halfspace. The theory governing the deformation of pre-stressed solid Earth subject to a normal surface traction of ice sheet relies upon the fundamental equations of motion and is discussed elsewhere (e.g., McConnell, 1965;Wolf, 1985;Ivins and James, 1999). For axisymmetric loading problems, it is possible to obtain the semi-analytical solution of vertical displacement at the lithosphere surface (i.e., both the ice-bed and ocean-bed interfaces). This is the essence of the solution for GIA assessment. For the equilibrium sea level hypothesis, however, the GIA solutions perturb the ice sheet only within the area of grounded ice. For the AIS that is surrounded on most of its periphery by floating ice, the extent of the grounded ice may evolve, as we assume in this study, according to the hydrostatic balance between the ice shelf and ocean water.
Given the appropriate ice loading history and choice of the model/material parameters (cf. Sect. 2.3), the GIA solution depends on the size of the ice load itself and the radial distance of the evaluation point from the load center. Assumed axisymmetry implies that the shape of the ice load be essentially cylindrical (e.g., Ivins and James, 1999). In the Cartesian framework of the ISSM/GIA model, we treat the size of the ice load as the property of a mesh element and compute the GIA solution at each node of the element. Individual 2-D (xy plane) mesh elements are defined as the equivalence of a footprint (i.e., projection onto the xy plane) of cylindrical disc loads, ensuring that the corresponding element and disc both share the same origin and plan-form area (cf. inset of Supplement Fig. S1a). The height of the ice load is then assigned to each element such that the average normal tractional force on the corresponding area of ice-bedrock contact is conserved. At each node within the domain, the final GIA solutions are computed by integrating the solutions due to individual disc loads, defined as the property of mesh elements.
The ISSM/GIA model is tested against the benchmark experiments (Wolf, 1985;Ivins and James, 1999) and found to be sensitive to the mesh resolution. For reasonably fine resolution (element size typically on the same order of magnitude as ice thickness, or two orders of magnitude smaller than the characteristic size of an ice sheet), however, the model performs well within the acceptable accuracy. Sample results are provided in Supplement Fig. S1.

Differential ice height
Prior to describing the ISSM/GIA model applied to the AIS, we briefly discuss how the change in Antarctic ice thickness is translated into the height of ice load. In this study, we as-sume that the sea level remains in its present-day state. We define the height of ice load at any time, t, as the differential ice height (DIH), h(x, y, t), with reference to the presentday configuration of the AIS (James and Ivins, 1998).
In the regions where ice is grounded both presently and at the given time t, DIH is simply the change in ice thickness over the corresponding time period. Similarly, in the regions where ice is floating both at present and at time t, DIH is zero as we assume that ice shelves are freely floating over the ocean, respecting the hydrostatic equilibrium. Defining DIH is a bit complicated in the areas over which the GL migrates over the course of time. If ice thickness at time t is smaller than the present-day value (i.e., ice floats at time t, but is now grounded), the DIH is defined as follows: where ρ w is the ocean water density, ρ is the ice density, b(x, y, t) < 0 is the isostatically corrected bedrock elevation (with respect to the present-day sea level) of the marine portions at time t (James and Ivins, 1998), and h(x, y, 0) is the present-day ice thickness. DIH in such a case would be negative. Similarly, in the areas where the ice sheet holds thicker ice at time t than at present (i.e., ice now floating is grounded at time t), DIH is defined as follows: where h(x, y, t) is the ice thickness at time t and b(x, y, 0) < 0 is the present-day bedrock elevation of the marine portions of the ice sheet. DIH in such a case would be positive. As a general convention in this study, we use t < 0 to denote the past, i.e., before present, and t > 0 for the future unless stated otherwise.

Model tuning
We apply the ISSM/GIA model to the AIS. We mesh the footprint of present-day AIS (Bamber et al., 2009) into triangular elements. In order to capture the potentially interesting features, the domain is discretized to consist of a high-resolution mesh around the coast (typical element size of 10 km) than in the interior of the ice sheet (element size of 25 km). This unstructured mesh captures the model inputs (e.g., past or future ice loads) in sufficient detail and provides a reasonable compromise between solution accuracy and computational efficiency. Doubling the mesh density, for example, improves the GIA solution (under present-day ice loading) only slightly (< 0.1 %), but it increases the high computational cost by one order of magnitude.
A comprehensive list of modern global positioning system (GPS) measurements is presented by . However, we tune our model by testing it against 18 highprecision data. Following Ivins et al. (2013), we first eliminate records from the Antarctic Peninsula north of 72 • S due to the associated difficulty of dealing with large elastic and transitional viscoelastic signals present there. We then average the values from stations located within 100 km of one another and eliminate some stations with reported errors greater than the signal amplitude.
In order to make reasonable predictions of the present-day GIA uplift, the slow response of highly viscous solid Earth demands that the evolution of AIS during the past several thousand years be considered in the ISSM/GIA model. There are a few reliable GIA ice loading histories for Antarctica (e.g., Peltier, 2004;Ivins and James, 2005;Whitehouse et al., 2012). These generally describe the timing and magnitude of deglaciation since the LGM based on geological and ice core data. In this study, we employ a much improved loading history discussed in Ivins et al. (2013). By upgrading the loading history of Ivins and James (2005) with recently available geochronological constraint data, the later model was able to provide a more accurate GIA correction to GRACE (Gravity Recovery and Climate Experiment) data measured over the period AD 2003-2012.
Based on Ivins et al. (2013), we have isostatically corrected DIHs available for 11 time stamps in the past (at −1, −2.2, −3.2, −6.8, −7.6, −11.5, −15, −17, −19, −21 and −102 kyr; see Supplement Fig. S2). Note that t = −21 kyr roughly corresponds to the LGM of the AIS, while −15 kyr marks the onset of deglaciation of the WAIS (e.g., Clark et al., 2009). For t = −1 kyr, we consider zero DIH that could be constrained using the recently available surface mass balance data (e.g., Verfaillie et al., 2012;Favier et al., 2013). However, this process is not straightforward because the magnitude and spatial distribution of ice flux during the periods of inferred mass balance are vastly unknown. Furthermore, we consider t = −102 kyr to mark the initial configuration for AIS as being identical to the present-day configuration. This implicitly assumes that the DIHs before the LGM have a minimal impact on the current and future response of the solid Earth. (We demonstrate in Sect. 3.2 that this is indeed a valid assumption.) Note also that the ice loading on the ISSM/GIA model is assumed to vary in a piece-wise linear fashion between the adjacent time stamps.
The model and material parameters considered in this study approximate the preliminary reference Earth model (Dziewonski and Anderson, 1981) Fig. S3) indicate that the majority of changes occur in the WAIS, we consider a 65 km thick lithosphere in our model. We cannot pick the corresponding mantle viscosity from Ivins et al. (2013), as our model does not have two mantle layers. We therefore consider the mantle viscosity as a tuning parameter, such that the difference in the mean between the measured modern GPS data  and modeled current GIA uplift at 18 data points is minimized (Fig. 2b). The optimized solutions for current uplift rate are shown in Fig. 2a. Key characteristic features of our predictions include greater uplift rate around the Mt. Ellsworth territory and a mild rate of bed subsidence in the interior of EAIS. Such spatial patterns of uplift rate essentially reflect signatures of the employed ice loading history (cf. Supplement Fig. S2). Note that the optimal predictions of uplift rate (Fig. 2a) correspond to a mantle viscosity of 7 × 10 20 Pa s. As expected, this magnitude falls in between the upper (2 × 10 20 Pa s) and lower mantle viscosity (1.5 × 10 21 Pa s) recommended for the chosen lithosphere thickness (Ivins et al., 2013). While the architecture of the ISSM/GIA model can capture high-resolution spatial variability in solid Earth material parameters, we do not experiment with lateral inhomogeneities in this study.

Future ice loading
The AIS mass change may become more dynamic in the future due to ice shelf melting (e.g., Pritchard et al., 2012;Depoorter et al., 2013;Rignot et al., 2013). The SeaRISE participating ice sheet models, primarily driven by melt-dominated forcing, quantified the future evolution of AIS under the so-called "R8 scenario", which is the proxy of representative concentration pathway emission scenario 8.5 (RCP 8.5) (Bindschadler et al., 2013;. The RCP 8.5 scenario represents an ongoing rise in emissions throughout the century, reaching 8.5 W m −2 at AD 2100 (e.g., van Vuuren et al., 2011). The radiative forcing associated with the RCP 8.5 scenario is loosely related to the R8 scenario via all three components of the SeaRISE model forcing, namely the surface climate, basal sliding and ice shelf melting. As these forcings are the ones that govern the future evolution of AIS, it is relevant in the present context to provide a brief account of them.
Firstly, the SeaRISE surface climate forcing follows a 1.5 × A1B scenario (IPCC-AR4, 2007) until AD 2200. (The A1B scenario generally describes a future world of rapid growth in economy and population that peaks in mid-century, and technologies that rely equally on both fossil and nonfossil sources of energy.) A mild increase in surface temperature, a total of 0.5 • C, is assumed during the period AD 2200-2500. Secondly, no sliding amplification is considered until AD 2100 assuming that the Antarctic surface temperature will remain below zero, thus ignoring the potential for surface melt-induced basal sliding prior to this time. Thereafter, sliding increases linearly at a rate of 20 % of its original value per century, but only in coastal regions. Inland, the sliding amplification factor decreases linearly as a function of surface elevation such that no sliding enhancement is applied above 1200 m a.s.l. Thirdly, ice shelf melting is assumed to  (Ivins et al., 2013). Black circles locate the position where model results are within the 1-σ uncertainty range of GPS measurements (Fig. 2b). Red circles represent cases where the model overestimates the measurements, whereas blue circles indicate locations where the model underperforms. Big circles are to denote the absolute misfits that are > 0.75 mm yr −1 . (b) Validation of the model against 18 high-precision GPS uplift data . Error bars depict 1-σ uncertainties associated with the GPS measurement.
increase linearly from its present-day value to 60 m yr −1 at AD 2200. Additional 10 m yr −1 melt extends linearly over the next 300 yr. Changes in basal melting conditions are only applied to the Amundsen Sea Sector (90-120 • W) and the Amery Ice Shelf (60-75 • E), not to the Weddell and Ross seas. Such a restriction of ice shelf melting considered in the SeaRISE experiment seems reasonable with reference to current observations (Depoorter et al., 2013;Rignot et al., 2013) that reveal that ice shelves around the Amundsen Sea are the most susceptible to melting. Furthermore, spatial distribution of ocean temperature anomalies (e.g., Pritchard et al., 2012) also supports this hypothesis.
Under the R8 scenario, a total of four ice sheet models simulates the future evolution of the AIS through AD 2500. These models are the Anisotropic Ice Flow Model (AIF) (Wang et al., 2012), the Penn State Ice Sheet Model (PennState) , the Potsdam Parallel Ice Sheet Model (PISM-PIK) (Winkelmann et al., 2011), and the Simulation Code for Polythermal Ice Sheets (SiCoPolIS) (Greve, 1997). These ice sheet models employ different assumptions and methods for solving the full physics involved in simulating ice flow (e.g., shallow-ice vs. first-order flow mechanics; different treatments for basal sliding, subglacial hydrology, and GL migration), they employ different numerics (e.g., different spatial and temporal resolutions), they have unique techniques for dealing with the prescribed model forcing, and so forth (cf. Bindschadler et al., 2013;. Consequently, each ice sheet model produces a unique evolution of the AIS. We extract the ice thickness from each model prediction for five time stamps (at t = 100, 200, 300, 400 and 500 yr) and compute DIHs using Eqs. (1) and (2) as appropriate. For simplicity, we assume b(x, y, t) ≈ b(x, y, 0) in Eq. (1) (i.e., future bedrock elevation is not isostatically corrected while com-puting future DIHs). We anticipate that the associated errors translated into the GIA solutions would be minimal. Examples of future DIHs are provided in Supplement Fig. S3. Again, ice loads in the ISSM/GIA model are assumed to vary linearly between the adjacent time stamps.

Future bed topography
In order to predict the future bed topography for AIS, the calibrated ISSM/GIA model (cf. Sect. 2.3) is forced by an appropriate sequence of ice load changes into the future. As we have four independent sets of future ice loading (cf. Sect. 2.4), we may compute four unique GIA solutions at any evaluation time in the future. Based on these solutions, here we present the estimates of future bed uplift, isolate the role of past and future ice loading, and also evaluate how predicted change in bed topography alters the bedrock slope. Note that in this section we deviate from the general time convention and present our results (predictions) in "AD".

Vertical bed displacement
The GIA solutions for individual future ice loading, combined with the consideration of a lone spin-up loading history, are computed at AD 2100 and AD 2500 and shown in Supplement Figs. S4 and S5, respectively. Although these solutions differ from each other in both magnitude and their spatial distribution, some common features are noteworthy: (i) all models predict minor subsidence in a few places, particularly along the Wilkes Land coast; (ii) the topography of the interior of the EAIS is likely to remain unchanged; and (iii) a pervasive and large uplift is predicted in the WAIS (except for the AIF simulations) and around the Amery Ice Shelf (except for the PennState and PISM-PIK simulations).
We can, of course, offer no one GIA solution as being more reliable than any other. We therefore calculate the average of model solutions (hereafter termed "model-average solutions") ( Fig. 3a and b) and consider these as reasonable estimates of the Antarctic bed uplift. It is possible that our ensemble approach for these predictions is insufficient. Nonetheless, we assert that they provide the correct orderof-magnitude estimates and the likely spatial patterns of the future bed uplift, which are sufficient to evaluate some of the IS/SE interactions outlined in Sect. 1.3.
By AD 2100, the Amundsen Sea Sector and the Amery Ice Shelf may rise by about four and three meters, respectively (Fig. 3a). The rest of the WAIS is likely to rise by up to two meters. The interior of the EAIS is predicted to remain unchanged. The Adelie and Wilkes lands, where all ice sheet models predict large snow accumulation , are likely to subside by about less than one meter. It should be noted that, for the chosen climate change scenario, all ice sheet models but SiCoPolIS predict moderate snow accumulation in the Queen Maud, Ross and Weddell basins, and minimal accumulation in the Amundsen and Amery basins (cf. . Roughly similar spatial patterns of GIA uplift are predicted for AD 2500 as well (Fig. 3b). In this case, the bed may uplift by about 25 m in the Amundsen Sector, and by about 10-15 m in the rest of WAIS and the Amery Ice Shelf. The interior of EAIS may remain mostly unperturbed. Bed subsidence of about four meters is likely to occur along the Adelie and Wilkes Land coasts, as well as along the coast north of the Amery Ice Shelf.
We also find similar spatial patterns for the model-average bed uplift rates (cf. Supplement Fig. S6). At AD 2100, the Amundsen Sea Sector and Wilkes Land are predicted to rise and subside at the highest rates of about 45 and −5 mm yr −1 , respectively. The interior of Marie Byrd Land is predicted to rise at the highest rate of 75 mm yr −1 at AD 2500; The Amundsen Sea Sector also rises at a large rate of about 40 mm yr −1 . The greatest rate of subsidence (of about −15 mm yr −1 ) is predicted along the east coast except around the Amery Ice Shelf.
At both evaluation times, as noted earlier, we obtain different GIA solutions associated with individual future ice loading. Here we briefly outline how well these predictions for the Antarctic bed uplift match one another. The standard deviation shown in Fig. 3c illustrates that the model predictions are generally in good agreement with each other at AD 2100, except around the Amundsen Sector and the Amery Ice Shelf. As depicted in Fig. 3d, similar agreement is found amongst the model predictions at AD 2500. Large deviations are predicted once again around the Amundsen Sea Sector. Moderate deviations can be seen along the east coast including the Amery Ice Shelf. Such deviations amongst model predictions for the Antarctic bed topography predicted both at AD 2100 and AD 2500 can generally be attributed to the limiting values of GIA solutions predicted by the PennState and AIF models (cf. Supplement Figs. S4 and S5). In both years, PennState predicts large uplift in the Amundsen Sector and the least uplift (in fact, subsidence) around the Amery Ice Shelf; the opposite is true in the case of AIF predictions.
Our predictions might slightly overestimate the GIA solutions in the EAIS, as the modeled lithosphere thickness (i.e., 65 km) is much thinner than the more common value (115 km, as reported in Sect. 2.3). Nevertheless, the following general findings should remain unaltered. With reference to the WAIS GIA solutions, (i) Amery Ice Shelf may rise moderately; (ii) the interior of the EAIS may remain mostly unperturbed; and (iii) minor subsidence may occur along the coastal EAIS. It also needs to be mentioned here that we do not model the subglacial erosion in this study, although quite rapid evacuation of soft sediments is now occurring at the bed of the Pine Island Glacier. Erosion has been roughly estimated to cause the topography to lower at a rate of 0.6 ± 3 m yr −1 as ice flow exceeding 1 km yr −1 erodes material in deep longitudinal fjord valleys (Smith et al., 2012). This erosive action typically takes place in 20 km wide valleys of approximately 200 km length. While an important consideration in ice sheet modeling (Kessler et al., 2008), GIA topographic responses occur over much broader length scales exceeding the areal dimensions of fast erosion by nearly two orders of magnitude, thus having an impact on the evolution of the ice sheet on the scale of the drainage basin itself.

Role of past and future loading
Our predictions of bed uplift reflect the combined effects of long-term viscous creep of solid Earth driven by the ice loading history and its short-term (centennial timescale) viscoelastic response to the future ice load change. It might be useful to isolate the contribution of past and future ice loading to the evolution of future bed topography. First, we let the calibrated ISSM/GIA model (driven by the past loading alone) run for the next 500 yr into the future under the idealized condition that the current distribution of ice thickness prevails as is, thus imposing h(x, y, t) = 0 m for all t ∈ [0, 500] yr. We find similar spatial patterns of bed uplift, as shown in Fig. 2a for the current uplift rate, in the future (cf. Supplement Fig. S7b); the total amount of GIA uplift at years AD 2100 and AD 2500 is in the respective ranges of about [−0.1, 0.5] (Fig. 4a) and [−0.5, 2.5] m (Fig. 4b). In such an idealized scenario of the unchanged future AIS, the notable features are that the peninsula, the whole of the WAIS, and coastal regions of the EAIS are likely to uplift, with the highest uplift occurring around the Mt. Ellsworth territory, and that the interior of the EAIS may remain unaffected or subside minimally (cf. Supplement Fig. S7b). These are consistent with characteristic features of the employed ice loading history (cf. Supplement Fig. S2).
Next, we compute model-average GIA solutions at years AD 2100 (Fig. 4c) and AD 2500 (Fig. 4d) by forcing the ISSM/GIA model by the future ice load changes alone. We impose h(x, y, t) = 0 m for all t ∈ [−120, 0] kyr to ensure that the model is properly spun up. Alternative solutions are found by subtracting the GIA solutions associated with the past loading alone ( Fig. 4a and b) from the corresponding final predictions ( Fig. 3a and b). The resulting solutions are essentially identical to those shown in Fig. 4c and d, implying that the principle of linear superposition holds (cf. Supplement Fig. S7).
Comparing the estimates of GIA uplift associated with the past (Fig. 4a and b) and future ice loading alone ( Fig. 4c and d) with those depicted in Fig. 3a and b, we find that the future ice loading dominates and that the contribution of longterm viscous creep associated with the past is only about one tenth the magnitude of the predicted GIA uplift. This suggests that the errors associated with the ice loading history may have minor consequences for the future predictions of the Antarctic bed topography, provided the differing scenarios properly sample the possible amplitude of the ice sheet loading. Significant changes in magnitude and timing of the loading history may, however, require that the different mantle viscosity be employed in the ISSM/GIA model to predict the current uplift rates correctly (cf. Sect. 2.3). This, in turn, may yield the different contribution (not necessarily higher) of long-term viscous creep of the solid Earth to predictions of the future bed uplift. Nonetheless, it is important here to highlight the need for better constraining the DIHs, particularly in the recent past (during the past 1 kyr).

Change in bed slope
Strong spatial variability in the future GIA uplift (cf. Sect. 3.1) implies that the Antarctic bed slope will be modulated in the future (e.g., Gomez et al., 2013;Konrad et al., 2013). We compute current and future bedrock slopes (associated with the model-average GIA solutions) following α b (x, y, t) = α 2 bx (x, y, t) + α 2 by (x, y, t), where α bx (x, y, t) and α by (x, y, t) are the x and y components of the bed slope, respectively. Figure 5b, for example, depicts the present-day bedrock slope of Antarctica. While this plot reveals the degree of bed steepness, it does not provide the information regarding the aspect of slope. It is important to identify whether the bedrock has forward or reverse slope, particularly while evaluating the role of GIA uplift in the marine ice sheet instability (to be discussed later). We therefore plot the current bathymetry of the AIS in Fig. 5a; in order to facilitate the discussion, we only consider the areas with bedrock below sea level. Notice in the figure, for example, the blue color around the interior of WAIS that illustrates the existence of reverse bed slope in those regions. We obtain the future bedrock topography by adding the current bedrock topography and the future GIA solution (Sect. 3.1). From the corresponding bedrock topography, we calculate the bed slope at AD 2100 and AD 2500. The changes in bedrock slope are then computed by subtracting the present-day slope (Fig. 5b) from the future bed slope. Sample results are shown in Fig. 5c for AD 2500. In the figure, we generally notice the reduced bed slope (apart from a few localized regions with enhanced slope) around the Amundsen Sea Sector and the Amery Ice Shelf where large uplift is predicted. The reverse bed in these regions will thus have generally less steep slope in the future. Similar spatial patterns (but small magnitudes) of change in bed slope are obtained for AD 2100 as well (results not shown).
Although the magnitudes of change are larger in the regions where the bed has experienced large uplift, such as in the Amundsen Sea Sector and in the Amery Ice Shelf, percent changes in bed slope are also significant around the Ronne and Ross ice shelves (Fig. 5d). The bedrock slopes beneath these ice shelves, for example, reduce by more than one percent at AD 2500. Under the circumstances when GL is to advance (this actually is the general case in the present context; cf. Sect. 4.2), change in bed slopes around the GL beneath the ice shelves may impact the magnitude of GL migration (cf. Sect. 1.1 and Fig. 1a) and thus ice dynamics.

Implications for ice sheet dynamics
Our predictions for the future evolution of Antarctic bed topography may influence the future dynamics of AIS. In this section, we specifically quantify the potential effect of the predicted change in bed topography on the gravitational driving stress (Sect. 4.1), GL migration (Sect. 4.2), and ice surface velocities (Sect. 4.3). Because our model is not yet capable of computing IS/SE interactions with full dynamic feedbacks, it should be noted that some of the results presented below are obtained by bootstrapping the relevant future bedrock topography. The general procedure includes the following. First, we consider the present-day settings (geometric setting and boundary conditions) of the AIS for our calculations. Next, we upgrade the geometry and relevant boundary conditions (e.g., basal friction while calculating ice surface velocities in Sect. 4.3) to account for the future GIA uplift and perform recalculations. Comparing corresponding results, we finally isolate the influence of GIA uplift.  Fig. 5a. (c-d) Model-average change in bed slope, α b (x, y), at AD 2500. Negative magnitudes of α b (x, y) imply that bedrock will have less steep slopes in the future.
In order to minimize the potential compounding effects on the ice sheet dynamics, we keep the present-day settings (particularly, ice thickness and thermo-mechanical boundary conditions) fixed for the further calculations. We therefore advise caution in overinterpreting any individual results, obtained from the present-day settings of AIS perturbed by the predicted GIA uplift, for AD 2100 and 2500. For clarity, we revert to the general time convention (cf. Sect. 2.2) in order to stress that the following results illustrate the potential effect of the predicted bed uplift, not the predictions of the total GIA effect, on the future dynamics of AIS.

Gravitational driving stress
Here, we discuss the GIA effects on the gravitational driving stress. We compute driving stress following τ d (x, y, t) = τ 2 dx (x, y, t) + τ 2 dy (x, y, t), where τ dx (x, y, t) and τ dy (x, y, t) are the x and y components of the driving stress, respectively. For i = x, y, we define τ di (x, y, t) = −ρ g h(x, y, t) α si (x, y, t) as the ith component of driving stress (e.g., Cuffey and Paterson, 2010), where g is the gravitational acceleration, α si (= ∂ i h+α bi ) is the ith component of surface slope, and ∂ i h is the thickness gradient in the ith direction. In order to isolate the GIA effect, as noted earlier, we use the present-day ice thickness, i.e., h(x, y, t) = h(x, y, 0), and its gradients in all calculations. Hence, the change in bed slope due to the GIA uplift (Sect. 3.3) is responsible for modulating the driving stress. Figure 6 shows the change in driving stress associated with the bed uplift predicted at years AD 2100 and AD 2500. In both cases ( Fig. 6a and b), we notice small but similar trends of change in driving stress. Relatively large changes are predicted at positions of larger bed uplift. Reduction in driving stress is evident around the Amery Ice Shelf, implying that the local surface slopes are likely to flatten in this region. Minor increments in driving stress are predicted in the area around Dome C. Complex patterns are predicted around the Amundsen Sea Sector; an extensive zone of reduced driving stress is surrounded by zones with enhanced driving stress (see particularly Fig. 6b).
Theoretically, for the given distribution of ice thickness, changes in bed slope and surface slope (and, hence, the driving stress) should generally be in phase for the topography with forward slope, but out of phase for cases with reverse bed slope. Due to the complex nature of the Antarctic bathymetry, it is an arduously difficult task to find a robust and consistent relationship between changes in bed slope and driving stress, particularly around the Amundsen Sector (compare, for example, Fig. 5c vs. Fig. 6b). Nonetheless, we generally find the reduction in local driving stress associated with the GIA uplift (compare, for example, Fig. 3b vs. Fig. 6b). Given the small order-of-magnitude predictions for change in driving stress (i.e., several hundreds of Pascal only), it is important here to note that, as will also be shown in Sect. 4.3, minor changes in driving stress may affect the ice sheet dynamics substantially, as ice velocities are directly proportional to the driving stress by about a power of three (e.g., Cuffey and Paterson, 2010).

GL migration
In this section, we evaluate the effects of GIA uplift on the GL. We employ the simple hydrostatic equilibrium criterion to identify the transition points by seeking the floating ice thickness, h f (x, y, t), such that Regions with ice thickness h(x, y, t) > h f (x, y, t) are assumed to be grounded and the rest to be floating. Because we use the present-day ice thickness, i.e., h(x, y, t) = h(x, y, 0), in all calculations, the extent of the current and future grounded ice (i.e., GL) is determined by the corresponding bathymetry, thus highlighting the influence of predicted GIA solutions. Using Eq.
(3), we compute the extent of the grounded ice for t = 0, 100, and 500 yr. Changes in GL are then identified by subtracting, in turn, the first solution from the latter two. Figure 7a, for example, shows the GL migration associated with the GIA uplift predicted for AD 2500. The mask shown in the figure primarily represents GL advance, implying that more ice will be grounded in the future due to the GIA effects. However, we also predict the minor GL retreat in a few scattered locations, particularly along the Wilkes Land coast, where the bedrock is generally predicted to subside partly due to the large snow accumulation simulated under the chosen climate change scenario . Figure 7b, for example, depicts the mask of GL retreat around the Shackleton Ice Shelf. Note that we also obtain similar but less extensive migration in GL associated with the GIA uplift predicted for AD 2100 (results not shown); there is however no evidence of GL retreat in this case.
Based on the measured ice surface velocities (Rignot et al., 2011), we locate more than 2700 ice flowlines (cf. Supplement Fig. S8) to quantify the magnitude of GL migration (mostly advance) associated with the predictions of GIA uplift at AD 2500. Results are shown for two important regions, namely the Ronne (Fig. 7c) and Ross (Fig. 7d) ice shelves. Figures reveal that the GL may advance by more than 100 km in these two ice shelves. Significant GL advance (by tens of km) is also predicted in the Amery Ice Shelf, Amundsen Sector, Larsen Ice Shelf, Brunt Ice Shelf, and in other regions. In a few locations, e.g., the Shackleton Ice Shelf (Fig. 7b), as noted earlier, we predict a minor retreat (by ≤ 10 km) of the GL.
Although the primary control of GL migration in our calculations is bedrock elevation (Eq. 3), it should be noted that the bedrock slope plays an equally important role (e.g., Lingle and Clark, 1985;Gomez et al., 2010), as summarized in Sect. 1.1. Extensive advancement in GL associated with the GIA uplift is generally consistent with what is expected when slope reduces in the reverse bedslope topography (compare Fig. 5a vs. c, for example, around the Ronne Ice Shelf). For bedrock with forward slopes, however, the advancement in GL can be explained by the GIA-induced increment in bed slopes (compare Fig. 5a vs. c, for example, around the Getz Ice Shelf). The minor GL retreat predicted, for example, in the Shackleton Ice Shelf is associated with the flattening (Fig. 5c) of the forward bed slope (Fig. 5a). Although we are able to show a systematic relationship between the bathymetry, change in bed slope, and the direction of the GL migration in a few cases, it is complicated to provide such one-to-one relationships consistently over the entire ice sheet. Nonetheless, the results indicate that the GIA effects generally support the extension of grounded ice (i.e., GL advance) in the future, thereby promoting the stability in the marine portions of the ice sheet that rest on a reverse bed slope.

Ice surface velocities
Finally, we analyze the influence of GIA uplift on the ice surface velocities. By solving the quasi-static thermomechanical problem of ice flow, we calculate the englacial velocity field of the AIS with and without GIA effects. We assume that higher-order physics based on the equations governing mass and momentum conservation (Pattyn, 2003) together with the constitutive relations for isotropic ice (Glen, 1955) describe the internal creep deformation of ice, and that a viscous law of friction (e.g., MacAyeal, 1993) governs basal sliding. We rely on a steady-state thermal problem identical for all simulations, based on present-day conditions. For simplicity, we do not consider the possibility of till deformation underneath the ice sheet. A description of ice rheology (Glen, 1955) and other common assumptions related to ice flow modeling can be found elsewhere (e.g., Cuffey and Paterson, 2010).
For the present-day setting of the AIS, we solve this problem of ice dynamics through diagnostic simulation of the 3-D ice flow capability of ISSM, satisfying a number of boundary conditions (e.g., Larour et al., 2012b). We impose (i) a stressfree condition at the ice-atmosphere interface, (ii) water pressure directing towards the ice sheet at the ice-ocean (peripheral) interface, and (iii) a sliding condition governed by the basal friction at the ice-bed interface. Zero friction is applied at the base of the ice shelf (i.e., free-floating condition), while basal friction under the grounded ice is inferred from InSAR (Interferometric Synthetic Aperture Radar) based surface velocities (Rignot et al., 2011) using a data assimilation technique (e.g., Morlighem et al., 2010). The basal friction pattern is similar to the one described in Morlighem et al. (2013).
We re-run the simulation for two additional cases, associated with the predictions of GIA uplift at years AD 2100 and AD 2500. In each case, we upgrade the bedrock and surface elevations of the grounded ice; the extent of the grounded ice (GLs) is also updated (cf. Sect. 4.2). For floating ice as well, we upgrade both bed and surface elevations so that all floating ice is in hydrostatic equilibrium, which ensures continuity of the driving stress at the GL. Applying the same boundary conditions discussed above (except in the newly grounded or floating areas), we compute the englacial velocity for both cases associated with the future GIA uplift. In areas presently floating that become grounded at t = 100 and 500 yr (cf. Fig. 7a), we update basal friction assuming that it is roughly equal to the gravitational driving stress . Similarly, we impose a freefloating condition in areas presently grounded that float in the future (cf. Fig. 7b). For a given vertical profile of an ice sheet, the maximum velocity is always observed at the surface. We therefore place our emphasis upon the GIA influence on the ice surface velocities. Using the simulation results discussed above, we compute the GIA-induced change in surface velocity associated with the predictions of GIA uplift at years AD 2100 (Fig. 8a) and AD 2500 (Fig. 8b). In both cases, we find similar patterns of change in ice surface velocity. Although the predicted changes are small, about two to three orders of magnitude smaller than the surface velocities themselves (Rignot et al., 2011), a systematic reduction in velocity is evident around the sheet-shelf margins. This suggests that the GIA effects generally contribute to decelerating the flow speed across the GL, and hence promote stability in the marine portions of the ice sheet. Note that we mask out the ice shelves in our figures, for we have no intention of making predictions in the ice shelves.
The predicted changes in surface velocity for the grounded ice can be interpreted as the combined effects of changes in driving stress (cf. Sect. 4.1 and Fig. 6) and the GL (cf. Sect. 4.2 and Fig. 7) associated with the GIA uplift. Around the Ross and Ronne ice shelves, the GIA-induced reduction in surface velocity is consistent with the GL advance. In other regions, e.g., the Amundsen Sea Sector and the Amery Ice Shelf, predicted reduction in ice velocity can be attributed partly to the GL advance and partly to the reduced driving stress. All in all, the effects of GIA on several aspects of ice dynamics (e.g., driving stress, GL, and ice surface velocities) are consistent in that the GIA promotes systematic stability in marine portions of the AIS in the future.

Conclusions
This study has examined the interplay between the ice sheet evolution and the solid Earth responses for the AIS. First, we compute the future uplift of the Antarctic bedrock using the calibrated ISSM/GIA model driven by the inferred and predicted evolution of the ice sheet. Next, we evaluate how such a response of the solid Earth impacts AIS dynamics.
Our calculations are based on several approximations of model physics and numerics; it is important to highlight some of these here. The GIA model describes a simple twolayer representation of the solid Earth; the model and material viscoelastic parameters are kept constant spatially and the viscosity is Newtonian. Our ice sheet model solves the quasi-static thermomechanical flow problem for higher-order mechanics; the GLs are determined by the hydrostatic equilibrium criterion. A more comprehensive exploration of the positive or negative IS/SE feedbacks is warranted in the future. There is much to be learned from GIA models that employ additional GPS data, possibly right in the heart of the Amundsen Sea Sector, where viscoelastic uplift rates may approach 40 mm yr −1 (Groh et al., 2012). Our computations of ice loading after the present-day rely on several ice sheet models driven by the melt-dominated forcing under the R8 scenario. The computed model-average GIA response provides our assessment of its impact on the ice sheet dynamics. While there are limitations to the data and methods employed, our research brings us to the following two important conclusions of broader interest.
First, the short-term viscoelastic response of the solid Earth to the future ice load change, rather than its long-term viscous response to the past loading, governs the future evolution of the Antarctic bed topography. The magnitude and spatial variability in the future bed uplift are therefore determined by the nature of the future evolution of the AIS. A larger uplift is expected, for example, where the ice sheet loses more mass, while its far-field consequences seem to involve a relatively small amplitude subsidence. Our calculations suggest that the Antarctic bed may rise by a few meters and a few tens of meters around the WAIS, particularly the Amundsen Sea Sector and the Amery Ice Shelf, at years AD 2100 and AD 2500, respectively. Minor subsidences of about one meter and a few meters are predicted along the Wilkes Land coast at the respective times, partly caused by the net accumulation in the climate scenario runs . The interior of the EAIS is likely to remain unchanged.
Second, a pervasive and large uplift predicted in the interior of the WAIS, a substantially marine-based ice sheet, has particular significance because it generally corresponds to the flattening of the reverse bed slope. This drives the GL forward and consequently promotes the stability of the ice sheet. Our calculations, based on the present-day setting of the AIS perturbed by the predicted GIA uplift, reveal that the GL may advance by more than 100 km in the Ross and Ronne ice shelves due to the predicted GIA uplift for AD 2500. This may reduce the future ice surface velocities across the GLs by several tens of meter per annum.
The conclusions summarized above indicate a negative feedback between the ice sheet evolution and the solid Earth response for the marine ice sheet. For areas with reverse bed slope, for example, loss in ice mass flattens the bed and drives the GL forward and hence decelerates the rate of mass loss. Although our model is capable of illustrating this mechanism systematically, accurate quantification of its significance requires a dynamically coupled IS/SE model including the intricate details of ice sheet stability to breakup. This negative feedback is consistent with the ice sheet and sea level simulations computed by Gomez et al. (2010Gomez et al. ( , 2013 wherein loss in ice mass reduces the local sea level due to self-gravitation and hence decelerates the rate of mass loss. For accurate simulations of the AIS on centennial timescales under the reasonable climate change scenarios, both the solid Earth and sea level changes proximal to the GL may need to be properly accounted for.
The Supplement related to this article is available online at doi:10.5194/se-5-569-2014-supplement.