The influence of realistic 3D mantle viscosity on Antarctica’s contribution to future global sea levels

The response of the Antarctic Ice Sheet (AIS) to climate change is the largest uncertainty in projecting future sea level. The impact of three-dimensional (3D) Earth structure on the AIS and future global sea levels is assessed here by coupling a global glacial isostatic adjustment model incorporating 3D Earth structure to a dynamic ice-sheet model. We show that including 3D viscous effects produces rapid uplift in marine sectors and reduces projected ice loss for low greenhouse gas emission scenarios, lowering Antarctica’s contribution to global sea level in the coming centuries by up to ~40%. Under high-emission scenarios, ice retreat outpaces uplift, and sea-level rise is amplified by water expulsion from Antarctic marine areas.


INTRODUCTION
With nearly 700 million people living in coastal areas and estimates of the cost of sea-level rise reaching ~14 trillion dollars by the end of the century, understanding the future contribution of the polar ice sheets to global sea levels in coming years is of critical societal importance (1,2).Ice-sheet changes take place on a range of spatiotemporal scales and depend on interactions with the surrounding atmosphere, ocean, and solid Earth (3).Large sectors of the Antarctic Ice Sheet (AIS) feature glaciers terminating in the ocean with their base lying below the local sea surface height (4).These marine-based sectors have long been considered susceptible to rapid, runaway retreat (5)(6)(7)(8)(9)(10).Satellite and in situ observations indicate that marine sectors of the AIS are currently changing with possible runaway retreat already occurring (11).Modern instrumental records do not resolve the types of major AIS changes we may expect in the future, but geological evidence suggests that substantial marine ice-sheet retreat has caused multi-meter global sea-level rise at times in Earth history when the climate was warmer than present (3,12,13).Projections of future ice loss in these regions differ substantially among studies, however, and dominate the uncertainty in future sea-level projections (1,(14)(15)(16)(17).
Past and long-term future marine ice-sheet evolution in Antarctica depends on the evolving elevation of the solid Earth beneath the ice and on sea levels at the grounding line (18)(19)(20)(21)(22)(23).Coupled ice sheetglacial isostatic adjustment (GIA) models (also termed coupled ice sheet-sea level models in the literature) show that the future evolution of the ice sheet is sensitive to the underlying Earth structure, with models characterized by lower viscosity mantle and thinner lithosphere producing more and earlier uplift and relative sea-level fall beneath ice-sheet grounding lines across Antarctica (20,21,23).These outcomes are more effective at slowing or stopping future ice-sheet retreat on millennial timescales because ice flux is strongly sensitive to the thickness of ice at the grounding line, which is in turn proportional to water depth, and a fall in sea level therefore reduces the ice discharge across the grounding line-an effect termed "the sealevel feedback" (Fig. 1) (18).Model results have also highlighted the importance of including low viscosity beneath individual glaciers on centennial timescales (24,25).These models adopted a spherically symmetric Earth structure that is unlikely to represent the strength of ice-Earth-sea level feedbacks across the Antarctic continent, or even across West Antarctica.Global Navigation Satellite System (GNSS) inferences and seismic imaging indicate large lateral [three-dimensional (3D)] changes in Earth structure across the Antarctic continent (26)(27)(28)(29)(30)(31)(32), including estimated mantle viscosities as low at 10 18 Pa s and a thinner lithosphere in parts of West Antarctica leading to viscous deformation on decadal timescales, rather than the millennial scales commonly adopted in Antarctic-wide modeling studies [e.g., (33,34)].In contrast, seismic tomography models and tectonic reconstructions (31,(35)(36)(37) suggest that much of the East Antarctic Ice Sheet (EAIS) sits atop a craton with cooler, more viscous mantle and thicker than average lithosphere (Fig. 2).
A recent study with a sea-level model with deformation of a purely elastic Earth coupled to a dynamic ice-sheet model showed the importance of elastic effects for the strength of the sea-level feedback on ice-sheet evolution (22).Not included in the study, however, is viscous deformation associated with modern and future ice loss, which has been shown to contribute up to 60% of the total predicted GIA by the end of the 21st century in low-viscosity regions such as the Amundsen Sea Embayment (38).Literature applying coupled ice sheet-GIA models with radially varying Earth structure and simplified Earth deformation treatments suggests that viscous deformation has the potential to affect future grounding line dynamics in West Antarctica (20,21,(23)(24)(25)39).However, Earth structure varies markedly in Antarctica across small spatial scales, and up to now, no studies have robustly assessed the impact of 3D Earth structure on the AIS and global, spatially variable sea level incorporating both the elastic and viscous responses.
Here, we bring together the latest advances in coupled ice sheet-GIA modeling and solid Earth geophysics, informed by geophysical and geological observables, to assess the influence of GIA and 3D Earth structure on AIS evolution and its contribution to global sealevel changes in the coming centuries.To capture lateral variations in Earth structure, we combine the most recent seismic and geodetic records in Antarctica (26,28,29,31) with a global seismic tomography model (40) to build a global 3D model of Earth's mantle viscosity and lithospheric thickness with 3-km spatial resolution beneath key areas of ice loss in Antarctica (Fig. 2; see Materials and Methods).This 3D viscoelastic Earth model serves as input to a 3D GIA model (41,42) that includes gravitational effects and rotation and deformation of an elastically compressible, viscoelastically deforming Earth (see Materials and Methods).The model includes migrating shorelines and state-of-the-art treatment of water expulsion (Fig. 1) in uplifting marine sectors of Antarctica that includes its impact on calculations of Antarctica's contribution to global sea levels (43,44).To capture the sea-level feedback on grounding line migration (Fig. 1), we couple this 3D GIA model to the dynamic PSUICE3D ice sheet-shelf model (15) with a new coupling algorithm based on (42) but with the addition of nesting in the ice model to reach 5-km resolution across West Antarctica, matched with grid refinement in the GIA model to 3 km (resulting in a global computational grid with ~28 million grid nodes) to accurately represent the GIA response (Fig. 2) (38).The ice-sheet model setup (15) is informed by modern satellite constraints on mass balance as well as paleo-sealevel records from past warm periods, and adopts improved climate forcing from regional climate modeling (see Materials and Methods) with more finely resolved meteorology relative to previous coupled ice sheet-global GIA simulations.Results exploring resolution and Earth structure dependence of GIA predictions (38) are used to inform our GIA model setup (see Materials and Methods).Our coupled modeling framework represents a substantive advance in computational state-of-the-art, capturing future ice-Earth-global sea-level interactions at unprecedented detail and complexity.
We perform a suite of simulations considering a range of Earth rheologies (rigid, elastic, 1D viscoelastic, and 3D viscoelastic; see Materials and Methods), climate forcings (Representative Concentration Pathway scenarios RCP2.6,RCP4.5, and RCP8.5), and ice physics assumptions.We repeat all scenarios with consideration of ice shelf hydrofracture and the marine ice cliff instability mechanism [henceforth MICI; (15)] and considering only marine ice-sheet instability without MICI (henceforth MISI).By combining state-ofthe-art treatment of sea level and ice physics, our findings quantify the importance of Earth rheology and structure for future projections of grounding line dynamics, solid Earth deformation, ice mass changes, sea levels, and coastlines around the globe.

Projections of Antarctica's future contribution to global sea-level changes
The contribution of the West Antarctic Ice Sheet (WAIS) and associated local bedrock elevation changes to global mean sea-level change (ΔGMSL) for the full range of simulations are shown in Fig. 3 and (left column) illustration of the sea-level feedback mechanism that dominates under low climate forcing whereby uplifting bedrock and lowering sea surface beneath marine sectors underlain by low viscosity reduce ice loss across the grounding line, leading to reduced far-field sea-level rise (bottom panel) compared to a scenario where these effects are excluded (top panel).(Right column) illustration of the water expulsion mechanism whereby exposed ocean areas in Antarctica continue to rebound (bottom panel), expelling water out of Antarctica and increasing far-field sea level more relative to the scenario with a rigid bed (top panel).note that both mechanisms are active in general, but under low climate forcing, the sea-level feedback dominates and GiA reduces WAiS contribution to GMSl rise.On the other hand, under strong climate forcing, the ice-sheet retreat is too fast to be strongly sensitive to the uplifting bedrock and the sea-level feedback is weaker so the ice sheet retreats in a similar way regardless of whether a rigid or 3d viscoelastically deforming earth model is adopted, but substantial ice loss and bedrock uplift occurs in exposed marine areas with the 3d deforming earth model.hence, water expulsion is the dominant effect, enhancing Antarctica's contribution to global mean sea-level rise.Figure by e. Goblet.table S1.With 3D viscoelastic Earth structure, under the RCP2.6-MISIscenario, the WAIS is projected to contribute 0.21 m by 2150 and 0.70 m by the end of the simulation at 2500.The contribution increases with higher emission scenarios, reaching 1.03 m for RCP4.5-MISI and 3.05 m for RCP8.5-MISI by 2500.When hydrofracture and ice cliffs are accounted for (i.e., with MICI), as in (15), the AIS contribution to ΔGMSL increases substantially, especially after the 21st century (Fig. 3).By 2500, the ΔGMSL contribution from WAIS reaches 1.68 m under RCP2.6-MICI,2.15 m under RCP4.5-MICI, and 5.75 m under RCP8.5-MICI.The ice loss in all scenarios largely comes from the WAIS, whereas the EAIS experiences a small ice gain due to increased accumulation (fig.S2), except for RCP8.5-MICI in which EAIS ice loss is activated and the contribution from the whole ice sheet exceeds 15 m by 2500.
The projected Antarctic ice loss in all scenarios produces a geographically and temporally variable pattern of global sea-level change due to gravitational, Earth rotational, and deformational (GRD) effects that amplify sea level rise along low to mid latitude coastlines (Fig. 4).For example, the ΔGMSL contribution for the whole AIS (including WAIS mass loss and small EAIS mass gain) under RCP2.6-MICI is 0.27 m by 2150 and 1.46 m by 2500, while peaks of up to 0.33 and 1.8 m occur in the Pacific and Indian Oceans at these times when accounting for GRD effects.A similar pattern of amplification occurs across the rest of the simulations, except for RCP8.5-MICI in which both EAIS and WAIS on either side of the rotational pole undergo ice loss, reducing Earth rotational effects (45) and bringing the regions of peak sea-level rise closer to the equator (Fig. 4F).Note that all simulations produce higher than average sea-level rise due to Antarctic ice loss throughout the next 500 years for all Small Island Developing States (also called Large Ocean States), in agreement with previous studies (46,47).This geographic vulnerability will be further amplified due to higherthan-average sea-level rise associated with Greenland ice loss [e.g., (47)].

Sensitivity of projected Antarctic sea-level contribution to Earth deformation and climate scenario
Comparing simulations adopting a range of different Earth structure models in Fig. 3, we find that GIA begins to influence the contribution of the WAIS to ΔGMSL on multi-century timescales for all climate and ice physics scenarios considered, starting as early as the end of the current century for RCP2.6.Uplift of the solid Earth and gravitational drawdown of the sea surface have the strongest effect when 3D viscoelastic Earth structure is adopted, where zones of thinned lithosphere and low mantle viscosity (Fig. 2) show faster, more localized bedrock uplift beneath areas of active WAIS marine ice retreat (fig.S6).With the 3D viscoelastic Earth structure, GIA reduces the WAIS contribution to ΔGMSL in all scenarios except RCP8.5-MICI.Differences between simulations with GIA and 3D Earth structure and those with fixed bedrock elevation (i.e., labeled "rigid"; from stand-alone icesheet model simulations; see Materials and Methods) reach over 50% during the simulations and are up to 40% for RCP2.6 and up to 25% for RCP4.5 cumulatively by 2500 (Fig. 3 and table S1A).Assuming a purely elastic Earth structure or a reference 1D viscoelastic model representative of continental average structure (see Materials and Methods and Fig. 2) neglects/underestimates the growing viscous uplift in low-viscosity zones, thereby weakening the sea-level feedback on ice dynamics and overestimating the WAIS contribution to ΔGMSL.ΔGMSL contribution is up to 34% and 17% higher with elastic and 1D viscoelastic Earth structures, respectively, relative to 3D Earth structure simulations.
The impact of GIA and 3D Earth structure on the ice sheet is larger and plays a role on ice dynamics sooner for lower emission scenarios.This is especially evident with MICI included, whereby the inclusion of 3D viscoelastic relative to rigid Earth structure reduces Antarctica's GMSL contribution by >50% with RCP2.6-MICI, and up to 25% with RCP4.5-MICI (table S1).Under RCP8.5, grounding line retreat outpaces bedrock uplift, with retreat across WAIS occurring within the first two centuries regardless of the adopted Earth model.However, after ice has completely retreated from WAIS basins, solid Earth uplift continues (fig.S4), driving water expulsion and increasing the WAIS contribution to ΔGMSL by 15% (0.8 m) by 2500 relative to the simulation with a rigid bed (Fig. 3F).In summary, we find GIA effects reduce the WAIS contribution to GMSL rise under lower emissions through the sea-level feedback on marine ice-sheet grounding line dynamics (18), but amplify it under high-emission scenarios through water expulsion (Figs. 1 and 3)

Influence of Earth deformation on line dynamics
The impact of the sea-level feedback on the grounding for a given ice physics and climate forcing scenario varies spatially between peripheral glaciers (e.g., compare the retreat of Thwaites and Pine Island glaciers in Fig. 5 and movies S1 to S8).Ice flow and grounding line migration depend on regional bedrock and ice etry, including local bedrock elevation changes beneath retreating grounding lines.In all simulations, especially those with no MICI, grounding line retreat rates are variable rather than monotonic, accelerating when backing steep reversed slopes, and slowing or pausing when reaching bumps pinning points (illustrated by grounding line positions in Fig. 5 and movies to S8).With the deformable bedrock, the grounding line remains on bedrock ger, slowing the retreat compared to simulations a rigid Earth.
example, with RCP2.6-MISI, the Thwaites grounding line retreats 134 km less 2500 than with the rigid bed.
Consistent with previous models (15), invoking MICI results in and more extensive Thwaites grounding line retreat for warming scenario compared to only scenarios.When GIA included in lower emission scenarios with MICI, accelerated rates of ice loss at the grounding line begin later and remain lower compared to rigid-bed simulations (fig.S3).Uplift of interior bedrock stabilizes the grounding line sooner than when the remains at its initial evation.RCP2.6-MICI with rigid bed, the grounding line retreats 40 km further inland by 2100, and 100 km further inland by 2200 compared to the simulation with 3D deformation.By 2300, the grounding line has retreated across most of the region in both simulations, but the retreat stops sooner on high located ~425 km inland from the initial grounding line with deformation included, while the grounding line continues to retreat beyond this high bedrock uplift is neglected (Fig. 5, A and D, and fig.S3).While extensive ice loss under RCP8.5-MICIleads to greater bedrock uplift than in other scenarios (figs.S4 to S6), the ice cliff retreat outpaces uplift, and most of the additional viscous uplift occurs after the grounding line has passed through the region, causing water expulsion and enhancing ΔGMSL (Fig. 4I and fig.S6D).

Impact of viscous deformation on bedrock elevation and sea-level change at the line and globally
Viscous effects contribute substantially to bedrock uplift and relative sea-level fall in Antarctica on decadal to centennial timescales (fig.S6).In simulations adopting a purely elastic Earth structure, peak uplift in the Amundsen Sea Embayment reaches 5 to 7 m by 2150 without MICI and 24 to 74 m with MICI.Peak deformation increases by 20 to 26 m without MICI and by 51 to 268 m with MICI at the same time points in simulations with 3D viscoelastic Earth structure (fig.S6).In scenarios with less/slower ice loss (e.g., fig.S6, A and C), substantially more uplift occurs with the 3D viscoelastic model compared to the elastic model during the period when the grounding line retreats across the basin.Conversely, with MICI and stronger climate forcing (fig.S6D), much of the grounding line retreat across the basin occurs before differences in the amount of uplift emerge between the different Earth models.Under RCP2.6 and with 3D viscous effects included, sea-level rise is up to 0.20 m lower in 2150 and 0.80 m lower in 2500 along global coastlines away from Antarctica (Fig. 4, C and D) relative to the Elastic Earth model or reference 1D viscoelastic Earth model (fig.S7).
Our results highlight the importance of considering lateral variability in Earth structure (Fig. 2A) in Antarctic-wide simulations (Fig. 6).With RCP8.5-MICI, substantial ice loss occurs in both East and West Antarctica.In West Antarctica, and especially the Amundsen Sea Embayment, low viscosities lead to greater uplift (Fig. 6A), reducing ice loss (Fig. 6B) compared to the simulation adopting an average radially varying Earth structure.In East Antarctica, higher than average mantle viscosities lead to less uplift beneath marine sectors, allowing for more ice loss compared to the average model.In particular, the bed remaining at lower elevations in the Wilkes and Aurora Subglacial Basins with the 3D Earth model allow for increased ice loss in these regions.The effects of this East to West contrast have also been explored in (23).Thus, purely radially varying Earth structure models do not accurately capture the deformation of the Earth and impact on the ice sheet in either East or West Antarctica.

Importance of the method of computing ice-sheet contribution to ΔGMSL
In agreement with recent studies (48,49), the method adopted to compute ice-sheet ΔGMSL contribution affects the degree of influence of GIA on sea-level projections.Here, we use a method based on total water conservation (see Materials and Methods) that accounts for the effects of ice loss and water expulsion.This method produces the same results as the "modified volume-above-floatation (VAF)" method defined in recent studies (48,50,51), correcting an important calculation flaw in the basic VAF method widely used in earlier literature.
Differences in ΔGMSL due to the method chosen are discussed in the Supplementary Materials.Figure S1 shows ΔGMSL projections using the basic VAF method, which differ substantially from those in Fig. 3 (except in the rigid simulation), especially early in the simulations.Other differences later in the simulations are due to the water expulsion effect (43,44), which is included in our method but not in basic VAF or (49) (see Supplementary Materials).Caution should thus be taken when comparing calculations of ice-sheet ΔGMSL contributions in the literature, and it is important to explicitly state the method adopted in a given study.

DISCUSSION
Our results indicate that GIA drives two competing effects that alter the AIS contribution to global sea levels: (i) the sea-level feedback on grounding line dynamics (18,19), in which solid Earth uplift and sealevel fall reduces grounding line retreat in marine basins in Antarctica, and (ii) the water expulsion effect (43,44,48), in which uplifting marine sectors freed of ice expel water out of Antarctica, increasing sealevel rise (Fig. 1).Here, we examine these effects together under a range of climate forcing and ice physics scenarios.When 3D variations in viscoelastic Earth structure across Antarctica are included in coupled ice sheet-GIA model simulations for the future, the first effect dominates in low-emission ice loss scenarios.Under RCP2.6, the sea-level feedback can reduce grounding line retreat by Ο(100 km) and reduce the contribution of the AIS to global sea-level rise by up to ~40% starting around 2100.These findings are consistent with regional studies and simplified treatments of GIA (23)(24)(25)(26), finding that rapid viscous deformation has the potential to reduce retreat of individual glaciers in the Amundsen Sea Embayment on centennial timescales.For high-end projections under RCP8.5, the rapid ice-sheet retreat across marine sectors outpaces GIA and the ice-sheet dynamics are relatively insensitive to the choice of Earth model.In this case, water expulsion is the dominant effect, increasing global sea-level rise on multi-century timescales.Our findings suggest that sea-level physics should be accounted for in future generations of ice-sheet model projections informing Intergovernmental Panel on Climate Change Assessment Report 7 (IPCC-AR7) (1,14).
Uncertainty in the AIS's contribution to future global sea-level projections remains substantial, with recent estimates ranging from minimal change to tens of meters of change on multi-century timescales.Uncertainties mainly stem from the choice of climate forcing (52,53), ice-sheet model differences ( 14), and Earth system feedbacks such as the solid Earth-ice interactions that are the focus here.The goal of our study is not to update future sea-level projections but rather to understand the influence of GIA feedbacks across a range of icesheet retreat timescales.Our work shows that GIA effects and the 3D rheological structure of the solid Earth play important roles in governing the response of the AIS to a warming climate with implications for impacts on global coastlines in the coming centuries.
Viscous deformation in our simulations takes place on decadal timescales in WAIS marine sectors (fig.S6), exceeding the contribution from elastic deformation to sea-level fall (bedrock elevation increase) near grounding lines and reaching tens to >100 m of uplift by 2150.We find that neglecting lateral heterogeneity in Earth structure results in substantive bias in the amplitude and spatial patterns of the projected ice loss, deformation, and projected global sea level in continent-wide simulations (Fig. 6 and fig.S7).Future work is needed to clarify to what degree radially varying Earth structure models can be applied to represent deformation and ice-Earth interactions accurately across regional and basin scales, and ambiguity remains on how best to define an optimally applicable regional 1D Earth model.Alternative inferences of 1D Earth structure to approximate GIA across regions of laterally heterogeneous Earth structure, such as West Antarctica, yield substantial differences depending on the methodology adopted and the observations considered (26,54) structure from the literature [extended data figure 8 from (15)] can produce different levels of grounding line stabilization.Recent seismic surface wave tomography suggests that there is heterogeneity in Earth structure on <100km spatial scales in the Amundsen Sea Embayment (55) that can produce notable differences in predicted GIA at basin scale modern and future grounding lines relative to continental 3D and 1D Earth models (56).Details of the local bedrock setting also matter, with spatiotemporally variable retreat and ice cliff calving rates depending on the detailed bed geometry (Fig. 5 and fig.S4).Comparison of ice-Earth interactions at different glaciers from our results (movies S1 to S8) suggests that the degree of applicability of 1D Earth models may be basin specific.
Our findings suggest that further constraints on Earth structure from improved seismic imaging and bedrock elevation changes from ongoing geodetic (GNSS) observations will be needed to provide further insight on this matter.Improved length and spatial coverage of observations will be important for both reducing uncertainty in ice-sheet projections and improving the interpretation of geophysical and geological observations in Antarctica.Achieving better resolution of Earth structure in coastal West Antarctica where seismic observations are sparse, and particularly in areas of active ice loss in the Amundsen Sea Embayment region, is critical to understanding of grounding line retreat patterns.
Sea-level changes along coastlines in our simulations are spatially variable, peaking in the Caribbean, along the western coast of Africa and in the Pacific and Indian Ocean basins (Fig. 4).In many of these regions, when accounting for 3D viscoelastic Earth structure globally, sea-level rise due to Antarctic ice loss reaches up to 1.7 m for low-end projections under RCP2.6 and up to 19.5 m by 2500 for high-end projections.Our results further support recent findings (46,47) that lowlatitude islands and coastal sites already being affected by sea-level rise will experience higher than average sea-level rise associated with Antarctic ice loss, regardless of the ice loss scenario.This finding highlights the climate injustice toward nations whose emissions are low, while their exposure and vulnerability to sea-level rise is high (46).We have found that the effect of GIA is to reduce the Antarctic contribution to global sea-level rise for low emissions and amplify it for high emissions.Therefore, reducing greenhouse gas emissions will allow the rebound of the solid Earth to play a greater role in preserving more of the AIS and avoiding the worst and most inequitable impacts of future climate change on global coastlines.

Experimental design
We use an ice sheet-sea-level coupling approach detailed and applied to the last deglaciation in (42) and adapted to future simulations here to model AIS thickness change and global GIA (i.e., including bedrock elevation changes relative to the geoid beneath the AIS, and gravitationally consistent global sea-level change) for a period of 550 years, from 1950 to 2500.In this approach, the PSUICE3D icesheet/shelf model (15,57) is coupled to the SEAKON 3D, global GIA model (41,42) that incorporates deformation of a range of 3D and 1D Earth structure models.We describe the two models and the coupling procedure below, along with the adopted Earth structure models below.

Ice-sheet model
We simulate AIS dynamics using the PSUICE3D model (58,59).The model uses a heuristic combination of shallow-ice and shallow-shelf approximations according to model equations and specifications described in detail in (58,59), with climate forcing and setup for future simulations described in (15).Key elements of the model setup, climate forcing, and details of the coupling procedure are described below.The simulations are run to 2500 as described in (15), following the extended RCP scenarios provided by (60).
The ice simulations are performed in two different modes.First, an initial continental-wide run with 10-km resolution is performed followed by a nested 5-km resolution simulation over West Antarctica (Fig. 1, boundary of the nested domain is marked with a black box).The latter resolution is the finest feasible over such a large area, and nested tests over the Amundsen Sea Embayment in (15) showed minimal resolution dependency with grid spacing ranging from 1 to 10 km.Both continental and nested modes start from a modern state in 1950 and extend until 2500, with 1-year time stepping.The continental simulations provide the necessary lateral boundary conditions for ice thickness and velocity in the nested model runs.Basal sliding coefficients are derived from surface velocities using the inverse methodology outlined in (61).Bathymetry and ice surface elevation at the start of the simulations are taken from Bedmap2 (62) [note that ice model tests not shown here show only minor differences between simulations adopting Bedmap2 and BedMachine (4) bedrock topography in these types of large-scale runs].Atmospheric and oceanic forcing for the ice-sheet model are provided as described in (15).Annual surface mass balance is calculated from monthly mean precipitation and surface air temperatures in a regional climate model (63), adapted to Antarctica.Meteorological climate forcing and subsurface temperatures determining sub-ice-shelf melt rates in the model follows three future Representative Concentration Pathway (RCP) scenarios: RCP2.6,RCP4.5, and RCP8.5 (60).For consistency with previous work, ice-sheet initial conditions (ice thickness, bed elevation, velocity, basal sliding coefficients, and internal ice and bed temperatures) follow the same procedure described (15), with a 100,000-year spin-up using observed climate forcing.Ice-sheet initial conditions are identical in all simulations.
Marine ice-sheet instability (MISI) dynamics are captured in all simulations.Half of the simulations incorporate the effects of hydrofracturing and mechanical failure of marine-terminating ice cliffs (MICI) as described in (15).These later simulations rely on two key, adjustable cliff-calving parameters (one controlling the ice thickness penetrated by hydrofracturing and the other controlling the maximum horizontal ice wastage rate of the cliff face), which we select based on the average calibrated values that are best able to reproduce Last Interglacial, Pliocene, and modern observations of ice loss as in (15).

3D global GIA model
To compute global sea-level changes and evolution of bedrock elevation beneath the AIS, we adopt a 3D, finite volume GIA model (41) that solves a generalized form of the sea-level equation that is described in detail (45), accounting for time-varying migration of shorelines (64,65), load-induced changes in Earth rotation (66), and the viscoelastic deformation of a self-gravitating, elastically compressible Maxwell (viscoelastic) Earth with 3D rheological structure.As input to the 3D GIA model, spatiotemporal evolution of grounded ice is provided by the ice-sheet model (see Materials and Methods) and a range of different rheological structures of the Earth are adopted, as discussed below.
Computations are performed on a global tetrahedral grid with triangulated spherical surfaces, and the model includes grid-refinement capabilities to achieve higher resolution across regions of West Antarctica characterized by active ice loss and low viscosity [see the supplementary materials from (42)].Results of sensitivity experiments in our previous work from (38) indicate that GIA predictions associated with modern and future ice loss converge for surface resolutions finer than 4 km, with small differences for surface resolutions less than 15 km, with the limiting factor the ability of the grid to capture the ice loading changes.On the basis of these findings, our computational grid reaches ~3km surface resolution over the refined area indicated by the green contour in Fig. 2, ~7km surface resolution over the rest of the Antarctic continent, and 12 to 15 km over the rest of the globe.Lateral resolution decreases with depth to ~50 km at the core-mantle boundary.Radially, there are 67 spherical layers spaced at similar distances to the lateral resolution of the corresponding triangulated surface to ensure relatively regular geometry for the tetrahedrons.The resulting grid is composed of ~28 million grid nodes and ~160 million elements.

Earth structure
We consider purely elastic (green lines in Fig. 3, labeled ELAS), radially varying (1D) viscoelastic and 3D viscoelastic models of Earth structure in our simulations.We also perform "rigid" simulations (blue lines in Fig. 3, labeled RIGID) with the ice-sheet model alone, where bedrock and sea surface elevations remain fixed in time.
As in (38), elastic and density structure in all GIA model simulations is radially varying and based on the seismic model STW105 (40).The reference 1D viscoelastic Earth model (adopted in simulations represented by orange lines in Fig. 3, labeled 1D) is characterized by uniform viscosities of 10 20 Pa s in the upper mantle (from the bottom of the lithosphere to the depth of 670 km) and 5 × 10 21 Pa s in the lower mantle (from 670-km depth to the core-mantle boundary) and a lithospheric thickness of 96 km, representing the average thickness across Antarctica in the 3D viscoelastic models described below.The mantle viscosities in the 1D model serve as the reference profile for perturbations in viscosities in the 3D viscosity model described below.
The 3D viscoelastic model we adopt (shown in Fig. 2, and adopted in simulations shown by red lines in Fig. 3, labeled 3D) is first derived and described in detail in (38) (see description of model EM3D_L).Viscosity variations relative to the 1D reference model described above are based on the ANT-20 seismic shear-wave speed model of the upper mantle described in (31) beneath Antarctica and the S362ANI seismic tomography model (40) across the rest of the globe.The S362ANI model has a resolution of (1000 km).The ANT-20 model resolves structures with a spatial scale of ~100 km down to 410-km depth and higher at greater depths.Variations in lithospheric thickness are based on the model of ( 67) globally and the model of (35) over Antarctica.As described in (68), the lithospheric thickness model is scaled to an average of 96 km over Antarctica, reaching a minimum thickness of 40 km in Marie Byrd Land.
Seismic velocities from the global and regional Antarctic tomography models are converted to viscosity variations following the procedure described in (41) (see their equations 27 and 29), by first converting seismic velocity anomalies into temperature anomalies, and adopting a scaling factor, ϵ, to relate temperature to viscosity.A scaling factor of 0.033°C −1 is adopted in Antarctica (ANT-20) and 0.04°C −1 for the global model (S362ANI) based on comparison to GNSS inferences of vertical viscosity structure beneath the Amundsen Sea Embayment (26), Fleming Glacier in central Antarctic Peninsula (29), and northern Antarctic Peninsula (28).More detail on this procedure is provided in (38).

Coupled ice sheet-3D GIA model simulations
We perform a suite of coupled ice sheet-global GIA model simulations, repeating results with the range of different adopted models of Earth structure described above, under extended RCP emission trajectories RCP2.6,RCP4.5, and RCP8.5 (60), both with and without MICI-related processes included.Note that sea-level changes due to Greenland and mountain glacier ice loss and thermosteric effects are not included in our model and will also alter sea levels around Antarctica.However, these contributions will be small at marine grounding lines compared to the local drawdown of the sea surface and uplift of the solid Earth associated with local ice loss (68,69).
The coupling procedure is as follows.First, the ice thickness changes across the simulation time from 1950 to 2500 are calculated from continental and nested simulations with the ice-sheet model alone with 1-year time stepping (see above).Then, the results are combined, with nested simulations defining WAIS ice cover changes and continental runs predicting ice loss over the rest of Antarctica, and passed as ice loading input for the 3D GIA model.The GIA model then computes the global GIA changes based on these ice cover changes at 2-year time increments across the simulation time, and these results are linearly interpolated to 1-year time increments.The 2-year time increments in the GIA model substantially reduce computation time, and initial testing of this procedure showed negligible differences between computing GIA every year, and every 2 years and then linearly interpolating.Changes in sea level (i.e., changes in elevation of the geoid relative to the solid Earth surface) are calculated globally and passed to the ice-sheet model over the ice-sheet model domain, where they are then used to provide changes in bed elevation every year in a new ice-sheet simulation.We repeat the coupling process three to four times to allow for convergence of predicted ice volume changes, following findings on convergence in (42).Figure S1 from (42) demonstrates convergence with our method for paleosimulations and shows that results with the iterative method are comparable to those produced with the standard interactive coupling procedure from earlier work with a radially varying GIA model (70).Note that while these tests suggest that our coupling procedure is appropriate for capturing the continental-scale ice-GIA interactions that are the focus of this study, further testing may be necessary for future applications focusing on more localized ice cover changes.

Global mean sea-level calculations
Estimates of the contribution of the AIS to global mean sea-level change (ΔGMSL), shown in Fig. 3, are calculated from the ice-sheet model results using a method based on the conservation of global water mass.It is appropriate in simulations with no major changes in ice sheets or shorelines in the rest of the world.Neglecting relatively small changes in other reservoirs such as groundwater (71), changes in total ocean, and ice mass within the Antarctic ice-sheet model domain (ANT) are equal in magnitude and of opposite sign to those in the rest of the world (REST).Also neglecting steric changes (71), this implies where A RO is the modern oceanic area outside the Antarctic ice model domain (REST) and ΔGMSL is the increase in sea level, i.e., ocean column thickness or the equivalent where there is floating ice, averaged over REST.The left-hand side with invariant A RO assumes that any changes in shorelines and ice distribution within REST are negligible.ΔGMSL is taken to be the Antarctic "contribution" to global sea-level change, because in our simulations there are no prescribed ice changes in the rest of the world.Σ is the sum over all grid cells in the relevant Antarctic region ANT (a weighting term equal to gridcell area is omitted for clarity in all summations here).The region ANT can be limited to West or East Antarctic sectors (separated by the pink line in Fig. 2A).Δ indicates the change from initial modern to current time.h is ice thickness (grounded or floating), G is the elevation of the geoid (coinciding with the sea surface where no floating ice), and h b is bedrock elevation, both relative to a fixed Earth reference level.The term in square brackets is the column thickness of liquid ocean water, equal to the radial distance from the seafloor to the base of floating ice, or to the sea surface if no ice.ρ i = 910 kg m −3 is ice density and ρ w = 1028 kg m −3 is ocean water density.Equation 1 yields exactly the same results as the "modified VAF" method recommended, for instance, in (48) Equation 2 is equivalent to that used in (48,50,72), except for small "max" corrections that are needed for generality.The expressions in large parentheses on the right-hand sides of Eqs. 1 and 2, i.e., the total water equivalent volume and the modified volume over floatation for individual grid points, are exactly equal for all points (terrestrial or marine, grounded, or floating ice).For grounded ice points [where (ρ i /ρ w ) h > G − h b ], both are equal to (ρ i /ρ w )h.For floating ice points or open ocean, both are equal to G − h b .For land points with no ice, both are zero.Therefore, the summed totals over ANT in each equation are equal.
For reference, the method using "basic" VAF is which is just Eq. 2 without the additional bathymetric term (+ max [G − h b , 0]) on the right-hand side.This calculation is applied to the global ocean area, i.e., A RO + A AO on the left-hand side (where A AO is the oceanic area in the Antarctic domain, using its modern value as in previous studies).The basic VAF method has been used as a standard in many previous publications.However, recent studies (48,49) show that without the bathymetric term, the basic VAF method is flawed and produces ΔGMSL estimates substantially different from the more accurate methods.The main flaw is that if marine ice columns remain grounded, changes in bedrock elevation with constant ice thickness cause basic VAF to change when there should be no effect on sea level.For comparison, fig.S1 shows ΔGMSL's for our simulations calculated using basic VAF.As expected, there are notable differences from Fig. 3. Differences between simulations adopting different Earth models generally emerge sooner and are larger when using the basic VAF method.Sensitivity tests show that this is nearly all due to the flaw in basic VAF described above.Therefore, early in our simulations before large grounding line retreat, changes in bedrock elevation below grounded marine ice spuriously affect ΔGMSL projections in basic VAF (fig.S1) but not in our results (Fig. 3).
Because our method includes the bathymetric term for all points including open ocean and floating ice, there is another difference from basic VAF: The "water expulsion" effect is included.This effect is produced by rising or falling bedrock elevations under open ocean or floating ice within the Antarctic domain, which displaces water into or out of the rest of the world oceans (43,44).The impact of water expulsion is particularly pronounced for RCP8.5-MICI with Earth deformation, causing ΔGMSL projections to be higher than those with rigid beds after ~2200 (Fig. 3C).In contrast, using basic VAF (without water expulsion), the ΔGMSL projections for RCP8.5-MICI with Earth deformation remain lower than those with rigid beds (fig.S1C).
Adhikari et al. (49) also correct the main flaw in basic VAF for ice that remains grounded, effectively in the same way as above (in their equation 11).However, they do not include the bathymetric term for open ocean or floating ice so that water expulsion is excluded from their calculation and their ΔGMSL estimate is an average for the entire world ocean area.This is consistent with the observation that the net effect of water expulsion on global mean sea level is zero, because it only redistributes water regionally and does not change global oceanic (plus floating ice) mass.The effects of water expulsion in Fig. 3 can be isolated by subtracting the ΔGMSL values from those using the method of (49) (not shown).The differences are small compared to the total in most cases [consistent with (44)], except for the later stages of RCP8.5-MICI with deforming Earth when water expulsion is pronounced, as mentioned above.
Whether or not to include the water expulsion effect in ΔGMSL calculations may depend on the focus of a particular study.We include it to assess the impact of changes in Antarctic ice, solid Earth, and sea levels on the rest of the world's coastlines; however, excluding the effect may be more appropriate for other applications, such as icesheet model intercomparison efforts, and to normalize predictions of spatially variable impacts of gravitational, rotational, and deformational (GRD) effects on global ocean depths.Regardless, it is important to explicitly state the method adopted in a given study.
Note that for the simulations with rigid beds (black curves), the ΔGMSL estimates for corresponding simulations in Fig. 3 and fig.S1 are nearly identical, i.e., the basic VAF method yields essentially the same results as the more accurate methods.That is because the bathymetric correction term needed to convert "basic" to "modified" VAF (+ max [G − h b , 0]) is invariant in time for simulations in which bedrock and geoid elevations are held constant; therefore, the ΔGMSL estimates with basic VAF are the same as with modified VAF and our calculations.In other words, the main flaw inherent in basic VAF does not arise if bedrock and sea-surface elevations remain constant [also noted by (49)].Furthermore, the water expulsion mechanism is not operative in this case.[There are still slight differences between the rigid-bed curves in Fig. 3 versus fig.S1, because the resetting of seasurface elevations in the simulations themselves is unphysical and does not conserve global water mass; i.e., the changes in sea level G − h b in Eq. 1 and the last term in Eq. 2 for Antarctic oceanic areas are (unphysically) zero, and so do not reduce the meltwater supplied to the rest of the world as they should, causing a slight overestimate of ΔGMSL in those equations.] Although not shown in Eqs. 1 to 3 above, all ΔGMSL calculations here include a "density correction" that accounts for the difference between fresh meltwater and ocean water (48,49), which amounts to subtracting the total change in solid ice volume (floating or grounded) multiplied by ρ i /ρ f − ρ i /ρ w from the right-hand sides, where ρ f = meltwater density 1000 kg m −3 .This correction is minor, ~5% or less of ΔGMSL and generally uniform across simulations, in the type of future Antarctic simulations performed here [sensitivity tests not shown, see also (48)].

Supplementary Materials
The PDF file includes: Figs.S1 to S7 table S1 legends for movies S1 to S8 Other Supplementary Material for this manuscript includes the following: Movies S1 to S8

Fig. 1 .
Fig. 1.Schematic illustrating the sea-level feedback and water expulsion mechanisms.(leftcolumn) illustration of the sea-level feedback mechanism that dominates under low climate forcing whereby uplifting bedrock and lowering sea surface beneath marine sectors underlain by low viscosity reduce ice loss across the grounding line, leading to reduced far-field sea-level rise (bottom panel) compared to a scenario where these effects are excluded (top panel).(Right column) illustration of the water expulsion mechanism whereby exposed ocean areas in Antarctica continue to rebound (bottom panel), expelling water out of Antarctica and increasing far-field sea level more relative to the scenario with a rigid bed (top panel).note that both mechanisms are active in general, but under low climate forcing, the sea-level feedback dominates and GiA reduces WAiS contribution to GMSl rise.On the other hand, under strong climate forcing, the ice-sheet retreat is too fast to be strongly sensitive to the uplifting bedrock and the sea-level feedback is weaker so the ice sheet retreats in a similar way regardless of whether a rigid or 3d viscoelastically deforming earth model is adopted, but substantial ice loss and bedrock uplift occurs in exposed marine areas with the 3d deforming earth model.hence, water expulsion is the dominant effect, enhancing Antarctica's contribution to global mean sea-level rise.Figure by e. Goblet.

Fig. 2 .
Fig. 2. 3D Earth structure beneath the AIS.(A and B) the constructed 3d earth viscosity structure based on seismic tomography(31) in Antarctica and (40) globally, plotted (A) at 120-km depth and (B) along a vertical cross section beneath the thwaites Glacier (shown by red lines in frame (A), and discussed in Figs.3 and 4), using the same viscosity color scale.the vertical axis above the solid black line in (B) has been amplified to show variations in ice thickness and bedrock elevation.the gray area below the solid black line represents the lithospheric thickness along the cross section in the 3d model, while the dashed black line represents the 96-km lithospheric thickness of the 1d reference model.viscosity variations are relative to a reference 1d radially varying viscoelastic model representing an Antarctica average with a viscosity of 10 20 Pa s in the upper mantle from the base of the lithosphere down to 670-km depth.the black box in (A) shows the higher-resolution nested ice-sheet model simulation region (see Materials and Methods), while the green line outlines the area of grid refinement in the GiA model.the pink line defines the WAiS-eAiS division adopted in sea-level contribution calculations presented in Fig.3and figs.S1 and S2.(C) lithospheric thickness in the 3d earth model in kilometers.the initial grounded ice thickness in meters and location of floating ice shelves at year 1950 in the ice-sheet model are shown over the ice-sheet model's (D) nested and (E) continental domains.Black lines in (d) and (e) indicate grounding line location.

Fig. 3 .Fig. 4 .
Fig. 3. West Antarctic ice loss and contribution to global GMSL for a range of adopted Earth structure models, ice physics, and climate warming scenarios.values given in the time series [plots in (A to F)] are from nested ice-sheet model simulations, just for the region of West Antarctica bounded by the pink line in Fig. 2A.Projected global mean sea level change (ΔGMSl) associated with WAiS evolution for six different ice physics and climate warming scenarios: RcP2.6 (A and d), RcP4.5 (B and e), and RcP8.5 (c and F), including Mici (A to c) and excluding Mici and considering MiSi only (d to F). colored lines represent simulations adopting four different assumptions of the earth structure as indicated in the legend.details of each of these models are described in Materials and Methods.ΔGMSl calculations shown here are described in Materials and Methods, and ΔGMSl calculations with a basic vAF approach are plotted in fig.S1.Beneath each time series plot in (A) to (F) are ice thickness for the 3d earth model simulations (left, with corresponding color bar on the top right of the figure) and difference in ice thickness between 3d and rigid simulations (right, with corresponding color bar on the bottom right of the figure) in 2500.Green and black contours correspond to the grounding line position with the 3d and rigid earth models, respectively.

Fig. 5 .
Fig. 5. Impact of solid Earth rheology on grounding line migration in the Thwaites Glacier basin in nested ice-sheet model simulations.(A to L) cross section along thwaites Glacier (location shown in Fig. 2A) indicating changes in the ice surface elevation and bedrock elevation associated with all modeled climate and ice physics scenarios as labeled.As discussed in the text, Mici indicates that marine ice cliff physics and ice shelf hydrofracturing are included in the ice model, while MiSi simulations exclude these physics and consider only ice-sheet physics.Simulations adopt either (A to c and G to i) the viscoelastic earth structure shown in Fig. 2 or (d to F and J to l) a rigid (i.e., nondeforming) earth.contours are plotted at 10-year time intervals from dark to light blue, with faster ice loss and bedrock uplift indicated by greater spacing between contours.colored vertical lines mark grounding line positions at 2100, 2200, 2300, 2400, and 2500 in the simulations.the gray vertical lines indicate the initial 1950 grounding line position.

.Fig. 6 .
Fig. 6.Impact of 3D Earth structure on predicted bedrock elevation and ice-sheet changes in continent-wide ice-sheet model simulations.(A) the difference in bed elevation change and (B) ice thickness between simulations adopting the 3d viscoelastic earth structure versus reference (average) 1d viscoelastic structure at 2500 for scenario RcP8.5-Mici, which includes full collapse of marine ice sectors in West Antarctica and substantial retreat of ice in east Antarctica.Gray letters indicate the locations of the Amundsen Sea embayment (ASe) and the Wilkes (WSB) and Aurora (ASB) subglacial basins.