Research ArticleGEOPHYSICS

Rapid postglacial rebound amplifies global sea level rise following West Antarctic Ice Sheet collapse

See allHide authors and affiliations

Science Advances  30 Apr 2021:
Vol. 7, no. 18, eabf7787
DOI: 10.1126/sciadv.abf7787


Geodetic, seismic, and geological evidence indicates that West Antarctica is underlain by low-viscosity shallow mantle. Thus, as marine-based sectors of the West Antarctic Ice Sheet (WAIS) retreated during past interglacials, or will retreat in the future, exposed bedrock will rebound rapidly and flux meltwater out into the open ocean. Previous studies have suggested that this contribution to global mean sea level (GMSL) rise is small and occurs slowly. We challenge this notion using sea level predictions that incorporate both the outflux mechanism and complex three-dimensional viscoelastic mantle structure. In the case of the last interglacial, where the GMSL contribution from WAIS collapse is often cited as ~3 to 4 meters, the outflux mechanism contributes ~1 meter of additional GMSL change within ~1 thousand years of the collapse. Using a projection of future WAIS collapse, we also demonstrate that the outflux can substantially amplify GMSL rise estimates over the next century.


GPS observations of rapid crustal uplift rates near zones of recent mass flux in the Amundsen Sea Embayment of West Antarctica indicate a viscosity in the shallow upper mantle beneath the region (~4 × 1018 Pa s) that is about two orders of magnitude less than typical upper mantle values, and a lithosphere of thickness 50 to 60 km (1). This inference is supported by seismic tomographic imaging of structure beneath West Antarctica, which reveals a thin lithosphere and a shallow mantle characterized by anomalously slow velocities bordering the cold, thick lithosphere of the East Antarctic craton (24). Furthermore, geological evidence such as widespread mafic volcanism and geophysical prospecting of the West Antarctic rift system (57), as well as mantle flow modeling of extension and dynamic uplift of the Admiralty Mountains (8), strongly support a thermal (i.e., high temperature) interpretation of these seismic velocity anomalies.

Following the collapse of grounded ice in marine-based sectors of the West Antarctic Ice Sheet (WAIS), viscoelastic crustal rebound will steadily reduce the volume of the accommodation space for meltwater from the retreating ice sheet and thus increase the total global mean sea level (GMSL) rise within the open ocean (914). Previous geophysical modeling of this effect in marine sectors of West Antarctica used one-dimensional (1D) (i.e., depth varying) viscosity models derived from studies of glacial isostatic adjustment (GIA) in cratonic settings (1012), which are characterized by a thick lithosphere and upper mantle viscosities of ~5 × 1020 to 10 × 1020 Pa s (15, 16). One simulation using these solid-Earth characteristics treated the case of a collapse of WAIS as equivalent to a GMSL rise of 3.26 m (10). They found that the amplitude of the crustal rebound would be small, with an instantaneous elastic contribution adding ~0.05 m to GMSL and an additional viscoelastic contribution that increases slowly from ~0.05 m in 2 thousand years (ka) to ~0.3 to 0.4 m after 10 ka (Fig. 1A, shaded region). Another simulation that only extended 500 years beyond a complete deglaciation of WAIS found that GMSL would differ only marginally from the elastic contribution during this interval (12).

Fig. 1 Time series of GMSL after WAIS collapse.

(A) Evolution of GMSL for a simulation of sea level change driven by WAIS collapse in the scenario B2009 (see Introduction), as described in (10). Solid line: Prediction based on the 3D viscoelastic Earth model V3DSD summarized in Fig. 2 (A and B) and an instantaneous deglaciation. Shaded region: The range of the viscoelastic contribution to GMSL reproduced from (10) for the case of a collapse of WAIS over a 500-year time scale. (B) Evolution of GMSL for two simulations of sea level change driven by the PSU3D1 scenario of WAIS collapse (17), as shown in Fig. 2F, with changes in EAIS masked out. Solid line: Prediction based on the 3D viscoelastic Earth model V3DSD summarized in Fig. 2 (A and B) and an instantaneous deglaciation. Dashed-dotted line: Same as solid line, with the exception that a 2000-year deglaciation was adopted. The y axis in both panels begins at the GMSL computed without the contribution from uplift of exposed marine-based sectors for the associated collapse scenario [3.20 m for B2009 in (A) and 3.76 m for PSU3D1 in (B)].

Ice sheet modeling studies have also considered the contribution to GMSL from postglacial uplift of marine-based sectors of WAIS [e.g., (17)]. With some recent exceptions (13, 14), these studies commonly refer to volumes of ice mass loss “above floatation,” which typically includes a time-dependent calculation of bedrock uplift beneath the ice and associated reduction in water accommodation space across the simulation. The bedrock uplift adopted in such calculations, however, does not generally include an instantaneous, elastic uplift of the crust, and the time scale of uplift, which is modeled as a simple exponential decay, is too long to appropriately reflect the low viscosity of the region (we return to this issue below).

The presence of anomalously low-viscosity shallow mantle beneath West Antarctica and the recent, direct observations of rapid uplift in response to recent ice mass loss (1) suggest that estimates of the additional GMSL rise driven by postglacial rebound of West Antarctica following WAIS retreat (and any related gravitational effects) require reappraisal. Here, we use these observations to reassess the GMSL response to WAIS retreat, focusing on the implications of the outflux mechanism for GMSL rise during interglacials as well as the impact of the outflux on GMSL for a projection of WAIS retreat over the next few centuries.

Many studies have explored the contribution of the AIS to sea level during interglacial periods (18) or in a warming world (10, 17). For example, the GMSL estimate of 3.26 m for WAIS cited above was based on a simulation that adopted the marine ice sheet instability hypothesis in an analysis of updated maps of Antarctic ice thickness and bedrock topography (10). This estimate is commonly cited as the maximum contribution of WAIS to peak GMSL during the last interglacial (LIG; ~130 to 116 ka) (1923). A total Antarctic contribution to the LIG higher than this would thus imply melt from the East Antarctic Ice Sheet (EAIS) (2426). Furthermore, when combined with independent estimates of GMSL changes associated with the Greenland Ice Sheet (18, 24, 2729), mountain glaciers (30), and thermal expansion (31), the estimate of 3.26 m has led to the view that the lower bound on a widely cited estimate of peak GMSL during the LIG [5.5 to 9 m; (19, 32)] requires no contribution from the EAIS, while it is likely that the upper bound does (33).

A recent and much broader set of possible Antarctic melt scenarios is provided by the Antarctic BUttressing Model Intercomparison Project (ABUMIP), which explored the response of a large number of ice sheet models to the total loss of ice shelves and the associated decrease in buttressing (17). The experiment was designed to explore the maximum possible impact of the marine ice sheet instability on AIS stability over a period of 500 years, and the results yielded a WAIS contribution to GMSL that ranged from 1.91 to 5.08 m. These values reflect the computed GMSL equivalent of the net change in the volume of ice above floatation across the 500-year simulation, and thus, they include a contribution to GMSL from uplift beneath marine-based sectors of WAIS. However, we obtained almost identical GMSL values, assuming no outflux of water from these sectors, using the initial and final states of a subset of the ABUMIP simulations, which indicates that the combination of the 500-year duration and relatively long crustal uplift time scale adopted in the experiments led to negligible uplift.

In the results below, we track the total GMSL changes over time using both the melt geometry constructed in (10), henceforth B2009, and the following subset of five ABUMIP simulations: PSU3D1, SICOPOLIS, f.ETISh, CISM, and GRISLI (17). In our sea level predictions, the latter five scenarios are combined with the bedrock topography models that were used in the associated ice sheet simulations: BedMachine (34) for f.ETISh, and Bedmap2 (Fig. 2C) (35) for the remainder. With the exception of one case discussed below, in which we consider GMSL change following melting from the EAIS, all predictions based on these melt scenarios mask out any ice mass changes within the EAIS. We use this suite of scenarios to explore the possible contribution to GMSL from meltwater outflux driven by uplifting marine-based sectors of West Antarctica, and the time scale across which this contribution is established, for a range of different scenarios for WAIS collapse. Our predictions in these cases will extend 10 ka after the collapse.

Fig. 2 3D viscoelastic Earth models used in calculations.

(A) Elastic lithospheric thickness and (B) mean viscosity from the base of the lithosphere to 400-km depth across Antarctica and the Southern Ocean for the Earth model V3DSD used in the majority of results described in the main text. Labels EL and MBL indicate the location of Ellsworth Land and Marie Byrd Land, respectively. The symbols indicate the location of GPS sites in inferences of mantle viscosity by (49) in the southern Antarctic Peninsula (squares) and (1) in the Amundsen Sea Embayment (triangles). (C) Antarctic bedrock topography from Bedmap2 (35). (D and E) As in (A) and (B) but with the alternative 3D viscoelastic Earth model V3DRH discussed in Materials and Methods. (F) Final grounded ice configuration in the ice sheet simulation PSU3D1 (17), with the red line showing the initial ice extent.

Last, we also consider the potential impact of the outflux mechanism on a projection of GMSL rise over the next few centuries. Specifically, we compute sea level changes associated with AIS evolution generated using the PSU ice sheet-shelf model (36) forced by the Representative Concentration Pathway (RCP) 8.5 scenario (37), with climate forcing implemented as in (38), and hydrofracture and ice cliff physics turned off such that minimal ice loss occurs in marine sectors of the EAIS. The PSU model adopts a hybrid combination of the shallow ice and shallow shelf approximations with a parameterization of flux across the grounding line from (39) and has been widely applied in recent paleo and future studies of the AIS, including the PSU3D1 scenario in (17). Our simulation produces 21st century WAIS ice above floatation within the range of model projections considered in (40) and a nearly full collapse of marine sectors of the WAIS within ~500 years.

Lateral variations in mantle viscosity exert a strong influence on the GIA process, and we consider a set of 3D Earth models to illustrate this sensitivity. The models are constructed by combining a chosen model of lithospheric thickness variations with a viscosity field derived by scaling high-resolution seismic shear wave tomography models for Antarctica. In the first model, described in detail in (41), lateral variations in elastic lithospheric thickness are taken from (42) within Antarctica and (43) elsewhere. The model lithospheric thickness varies from ~50 to 70 km across much of West Antarctica to greater than 200 km in cratonic regions of East Antarctica (Fig. 2A). We generate lateral viscosity variations using shear wave velocity (VS) anomalies from two high-resolution regional tomography models [(44) in WAIS; (45) in EAIS] and a global background model S40RTS (46). For each seismic model, VS is converted into viscosity using a depth-dependent linear scaling factor (47), and the resulting viscosity fields are merged using a smoothing algorithm applied across a region within ~200 km of the interface between the fields (41, 48). The radially averaged viscosity from the base of the lithosphere to a depth of 400 km beneath West Antarctica varies from 6 × 1018 to 2 × 1020 Pa s (Fig. 2B).

The average viscosities of this model in the vicinity of the Amundsen Sea Embayment (~1019 Pa s) and in the southern sector of the Antarctic Peninsula (~1020 Pa s) are consistent with the inferences based on analyses of GPS observations at sites within these regions in (1) and (49), respectively. Barletta et al. (1) decomposed the shallow mantle viscosity into two layers: the first extending from the base of the lithosphere to a depth of 200 km, and the second, from a depth of 200 to 400 km. Their inferences of viscosity in these depth zones were found to be strongly correlated, and thus, only the average, ~1019 Pa s, was uniquely constrained by data. Note that the low-viscosity areas beneath West Antarctica coincide with the grounded marine-based sectors of WAIS, including the Amundsen Sea Embayment, Ellsworth Land, and Marie Byrd Land (Fig. 2, A to C). We henceforth label this 3D viscosity model as V3DSD.

The results below also include sea level predictions based on three other 3D viscoelastic Earth models. The first two are identical to the Earth model V3DSD described above, but the mapping from seismic velocities to viscosity is altered so that peak lateral variations in mantle viscosity from the background 1D viscosity model (see Materials and Methods) are decreased (V3DSD–) and increased (V3DSD+) by a half an order of magnitude, respectively. A final Earth model (V3DRH) is constructed from an independent set of seismic constraints and a different procedure for mapping seismic velocity anomalies into viscosity and lithospheric thickness, as described in Materials and Methods (Fig. 2, D and E).


WAIS collapse above a low-viscosity upper mantle

We compute sea level changes due to GIA over 10 ka using a theory that accurately accounts for water influx into regions vacated by grounded, marine-based ice, as well as perturbations in the Earth’s gravity field, crustal elevation, and rotational state (see Materials and Methods). To begin, we adopt two melt scenarios characterized by a nearly full collapse of marine-based sectors of WAIS (Fig. 2C). The first is the scenario that we have termed as B2009 [(10), see their fig. S2], and the second is the PSU3D1 scenario in (17). The GMSL rise that we compute under the assumption of no uplift of exposed marine sectors of West Antarctica is 3.20 m for B2009, which approximately agrees with the value of 3.26 m cited in (10). The analogous value for the PSU3D1 scenario is 3.76 m, and the postcollapse ice cover in this second scenario is shown in Fig. 2F.

To isolate and assess the time scale over which viscoelastic uplift of marine-based sectors of West Antarctica will continue beyond the melt phase, we initially adopt an instantaneous deglaciation of these sectors. To compute GMSL change for each scenario, we track the total volume of meltwater leaving West Antarctica over the course of the simulation and divide this by the area of the open ocean, where the open ocean is defined as the area outside the initial region covered by grounded marine-based ice. The instantaneous elastic rebound of exposed, marine-based sectors of West Antarctica increases the net GMSL at the end of the melt phase by ~0.2 m for both melting scenarios (Fig. 1, A and B, solid lines). Subsequent viscoelastic uplift of these sectors adds another ~0.8 m. Thus, the net contribution to GMSL due to the outflux of meltwater from uplifting marine sectors of West Antarctica is ~1 m, bringing the total GMSL rise for the B2009 and PSU3D1 scenarios to 4.20 and 4.78 m, respectively.

Following deglaciation, viscoelastic uplift of marine-based sectors of WAIS occurs in localized areas coinciding with the zones of ice melt. Figure 3 tracks this uplift pattern for the PSU3D1 scenario. As a consequence of the low-viscosity setting of the West Antarctic region, the peak magnitude of the predicted sea level fall in West Antarctica in this case rapidly increases to 435 m in 500 years and ultimately plateaus at a value of 610 m (Fig. 3, B to D). This rapid uplift is reflected in the GMSL rise (Fig. 1B), which reaches 81% of the total (1.02 m) contribution from the outflux mechanism in just 1000 years (0.83 m).

Fig. 3 Postglacial sea level change after WAIS collapse.

Snapshots of predicted sea level changes (A to D) 0, 500, 2000, and 4000 years after an instantaneous collapse of marine-based sectors of WAIS, based on the PSU3D1 scenario [(17); see Fig. 2F]. Calculations were performed using the 3D viscoelastic Earth model V3DSD summarized in Fig. 2 (A and B). The thin blue line on each frame shows the areal extent of ice that melts in the PSU3D1 ice sheet simulation [(17); see Fig. 2F].

The consistency between the predicted contribution to GMSL from the outflux mechanism using the B2009 and PSU3D1 melting scenarios indicates that any situation involving a full collapse of marine-based sectors of WAIS will raise GMSL by ~1 m above previous estimates over a time scale of ~1000 years. This is significantly greater than previous estimates of the contribution (e.g., Fig. 1A, shaded region). We conclude that the peak GMSL rise of 3.3 m that is commonly cited in the LIG literature for this collapse (1923) is a substantial underestimate of up to 30%.

The sea level fingerprint of WAIS collapse

To investigate the geometry of WAIS meltwater distribution, we next consider the total sea level change across the global ocean immediately following the instantaneous melt event and 2 ka later, predicted using the PSU3D1 melt scenario (Fig. 4). The meltwater redistribution in the first case reflects a purely elastic response, in which sea level falls in the vicinity of the melting ice sheet due to the combined effects of postglacial elastic uplift of the crust and the reduced gravitational pull of a smaller WAIS on the ocean (Fig. 4A) (1012). Beyond this zone, sea level generally increases by progressively larger amounts with increasing distance from WAIS due to the gravitationally driven migration of excess meltwater away from Antarctica (i.e., meltwater that is not captured within local marine sectors). The far-field sea level rise reaches a maximum of ~5.1 m in the northeast Pacific Ocean and has a secondary peak of ~5.0 m in the Indian Ocean. The location of these peaks largely reflects the feedback on sea level of load-driven changes in the orientation of Earth’s rotation axis, which is characterized by a ~300-m shift of the southern rotation axis toward West Antarctica. Two thousand years later (Fig. 4B), viscoelastic adjustment produces a localized zone of subsidence around the periphery of WAIS characterized by a peak sea level rise of ~31.4 m (note that the color bar is saturated at this magnitude). Outside of this area, viscous effects significantly reduce the gradients in sea level that were predicted in the purely elastic response, and a peak far-field sea level change of ~5 m is maintained in both the northeast Pacific and Indian Oceans.

Fig. 4 Global sea level changes after WAIS collapse.

Predicted sea level changes (A) 0 and (B) 2 ka after an instantaneous collapse of marine-based sectors of WAIS based on the PSU3D1 scenario (Fig. 2F) (17) and computed using the 3D viscoelastic Earth model V3DSD summarized in Fig. 2 (A and B).

This result has implications for the relative sea level change at specific sites across the interglacial, as illustrated by predictions at a far-field site that do, or do not, include the expulsion of water driven by crustal uplift within West Antarctica (Fig. 5, solid and dashed lines). As noted above, gravitational and deformational effects associated with the modeled collapse would initially lead to a migration of meltwater away from West Antarctica (Fig. 4A), explaining the early sea level rise of ~5 m at the site in both predictions. Subsequent global-scale viscous adjustments and, in particular, the subsidence of the peripheral bulge offshore of West Antarctica would reduce far-field sea level over time (50). This relaxation leads to the rapid, monotonic sea level fall predicted in the case without water expulsion (Fig. 5, dashed line). In contrast, when the water expulsion effect is included, the sea level fall due to the global-scale viscous relaxation is almost entirely compensated by a sea level rise associated with the ongoing water flux out of West Antarctica, and the predicted sea level remains relatively constant across the simulation (Fig. 5, solid line). This raises an important issue. One would expect that far-field sea level rise at a time well after deglaciation (~4.73 m in Fig. 5, solid line) should provide a relatively accurate reflection of the total GMSL. The GMSL rise computed by incorporating the outflux mechanism (4.78 m) does so, but a calculation of GMSL that ignores this contribution (3.76 m) does not.

Fig. 5 Sea level change at Eleuthera, Bahamas, in the far field of WAIS.

Solid line: Predicted sea level change at Eleuthera, Bahamas, in the far field of WAIS for the simulation based on the PSU3D1 melt scenario (Fig. 2F) (17) and the 3D viscoelastic Earth model V3DSD (Fig. 2, A and B). The dashed line is analogous to the solid line, with the exception that only ice above floatation before the deglaciation is melted, and no outflux of water from rebounding marine sectors of WAIS is permitted (in this case, the GMSL rise is 3.76 m across the entire simulation; see The sea level fingerprint of WAIS collapse).

Sensitivity tests

In this section, we investigate the sensitivity of our results to various aspects of the sea level predictions. First, we repeat the prediction based on the 3D Earth model V3DSD and the PSU3D1 scenario for WAIS collapse shown in Fig. 1B (solid line) but adopt a melt duration of 2 ka (dashed-dotted line on the same figure). This time scale is commonly observed in climate simulations of the LIG (18). GMSL once again peaks at 4.78 m in this case, and by the end of the melt phase, the contribution to GMSL associated with the outflux mechanism (0.91 m) has already reached 89% of the final contribution (1.02 m).

Next, we consider the sensitivity of predictions based on the WAIS collapse scenario PSU3D1 to variations in the adopted 3D viscoelastic Earth model (Fig. 6A). As we have noted, the Earth models V3DSD– and V3DSD+ were constructed from V3DSD by decreasing and increasing the magnitude of peak-to-peak lateral viscosity variations superimposed on the background 1D viscosity model by half an order of magnitude, respectively. Thus, in generating the Earth model V3DSD+, viscosity values in V3DSD that are higher than the background 1D model are increased, while those below the background value are decreased. The predictions of GMSL are relatively insensitive to this level of perturbation in the viscosity field (Fig. 6A). The sensitivity is higher when the Earth model V3DRH is considered. A comparison of Fig. 2B and Fig. 2E indicates that the viscosity of the latter is somewhat higher than the former across the region of melt in the PSU3D1 scenario (Fig. 2F; see also fig. S1), and the result is a moderately slower-paced contribution to GMSL from the outflux mechanism.

Fig. 6 Time series of GMSL after WAIS collapse for different Earth models and melt geometries.

(A) Evolution of GMSL for four simulations of sea level change driven by WAIS collapse scenario PSU3D1 (Fig. 2F) (17). Each simulation adopts a different viscoelastic Earth model. Solid, dashed, and dotted blue lines: Earth models V3DSD (Fig. 2, A and B), V3DSD+, and V3DSD–, respectively (see Sensitivity tests). Solid purple line: Earth model V3DRH (Fig. 2, D and E). (B) Evolution of GMSL due to uplift of marine-based sectors for five scenarios of WAIS collapse taken from (17), as labeled, and the 3D viscoelastic Earth model V3DSD (Fig. 2, A and B). The GMSL computed without the outflux mechanism (i.e., without uplift of exposed marine-based sectors) for each simulation is listed in the legend.

Last, we repeat the calculation of GMSL change based on the V3DSD Earth model for four additional WAIS collapse scenarios described in (17). Figure 6B shows predicted time series of the contribution to GMSL from the outflux mechanism for each of these scenarios, together with the result for the geometry PSU3D1. These results demonstrate, once again, that any scenario involving a near-full collapse of marine-based sectors of WAIS (e.g., PSU3D1, SICOPOLIS, and f.ETISh) will be characterized by a rapidly increasing contribution to GMSL from the outflux mechanism that peaks at ~1 m, or ~27 to 28% greater than each GMSL computed without this outflux. Expectedly, melt scenarios that only involve a partial collapse of WAIS (e.g., CISM and GRISLI) experience a smaller increase in GMSL associated with the outflux mechanism, for example, ~0.43 m or ~23% of the GMSL computed without accounting for uplift of any exposed, marine-based sectors of West Antarctica in the case of GRISLI. Analogous numbers for the CISM melt scenario are ~0.34 m or ~19%. Nevertheless, the time scale over which this contribution is established remains rapid and of the order ~1 ka.

Relevance to future melting scenarios

Our results also have relevance for projections of GMSL change in the future, warming world. In previous millennial time scale projections [e.g., (51, 52)], estimates of GMSL rise associated with full collapse of WAIS will be amplified by ~1 m over a time scale comparable to the duration of the collapse event. Moreover, the results in Figs. 1 and 6 indicate that the outflux mechanism may also be important over centennial time scales. To consider this issue, we perform a final calculation of sea level change associated with an ice sheet model (36) projection of AIS retreat under the RCP 8.5 greenhouse gas concentration trajectory (37, 38). The model simulation has a 900-year duration beginning in 1950, with much of WAIS retreating due to marine ice sheet instability by ~2500. Ice thickness at the start and end of the simulation are shown in Fig. 7A and Fig. 7B, respectively.

Fig. 7 Projected GMSL change in response to AIS mass flux simulated with RCP 8.5 climate forcing.

(A and B) Ice geometry at the start and end of an ice sheet simulation forced by a high-emission (RCP 8.5) scenario with the PSU ice sheet model (36). The red line in (B) shows the perimeter of the ice in frame (A). (C) Solid line: Evolution of the total GMSL computed using this ice history and the 3D viscoelastic Earth model V3DSD (Fig. 2, A and B). Dashed line: GMSL rise when the outflux mechanism is not included. The inset focuses on the GMSL results within the period 2000–2100 CE.

The total GMSL rise computed using the 3D viscoelastic Earth model V3DSD is plotted in Fig. 7C, together with the GMSL time series in the case where the outflux mechanism is not included. At the end of the simulation, the total GMSL reaches a value of 3.68 m, ~19% (0.59 m) higher than the case without crustal uplift and water outflux (3.09 m). The inset in Fig. 7C focuses on results for the 21st century. Across this century, the outflux mechanism amplifies the GMSL computed in the absence of this mechanism by 18%. While a more comprehensive consideration of ice physics, climate forcing, and ice and Earth model parameters would be required to fully assess the implications of this effect on future GMSL, our results suggest that even century-scale projections of GMSL change should account for the volume of water driven out of West Antarctica by crustal rebound.


The additional contribution to GMSL from water flux out of West Antarctica has been neglected or significantly underestimated both in discussions of the GMSL rise budget during the LIG (1922) and in the values of GMSL change cited in coupled climate ice sheet studies of WAIS stability during the LIG [e.g., (18, 24)]. This suggests that efforts to bound contributions of individual sources to peak GMSL during the LIG may require reappraisal. We have also demonstrated that the outflux mechanism has a substantial impact on site-specific relative sea level histories (Fig. 5). An accurate treatment of the signal in sea level modeling of the LIG requires not only that the exposure of marine-based sectors be included in the modeling (912, 18) but also that areas of low viscosity in the shallow mantle beneath much of WAIS be accounted for. These conclusions also hold for previous interglacials in which marine-based sectors of the AIS destabilized (53), including those in EAIS.

A comparison of Antarctic bedrock topography (Fig. 2C) with the viscosity fields associated with the 3D Earth models V3DSD and V3DRH (Fig. 2, B and E) indicates that the marine-based sectors of East Antarctica, and in particular the Wilkes and Aurora basins, are relatively close to areas of low viscosity in the shallow mantle. To explore this issue further, we tracked GMSL following the two melt events for EAIS that occur in the PSU3D1 and SICOPOLIS scenarios (in these calculations, any melt from WAIS is masked out of the scenarios; Fig. 8). The GMSL rise computed without including the outflux mechanism is 1.44 m for PSU3D1 and 6.40 m for SICOPOLIS. The additional GMSL that accumulates from the uplift of exposed marine sectors is 0.08 m for the former and 0.57 m for the latter. In contrast to WAIS collapse scenarios, a substantial fraction of this outflux occurs in concert with the instantaneous melting: 0.05 and 0.28 m for the PSU3D1 and SICOPOLIS scenarios, respectively (Fig. 8). In addition, the subsequent rise in GMSL is slower paced than predictions associated with WAIS collapse events (e.g., the solid lines in Fig. 1), reflecting the comparatively less extreme viscosity values in the shallow mantle beneath East Antarctica.

Fig. 8 Time series of GMSL after EAIS collapse.

Evolution of GMSL due to uplift of marine-based sectors for two scenarios of EAIS collapse taken from (17), as labeled, and the 3D viscoelastic Earth model V3DSD (Fig. 2, A and B). The GMSL computed without the outflux mechanism (i.e., without uplift of exposed marine-based sectors) for each simulation is listed in the legend. All mass changes within WAIS for the two melt scenarios are masked out.

We have shown that GMSL rise following any period of collapse of marine-based sectors of AIS would be rapidly amplified by water outflux driven by postglacial rebound. This mechanism is especially important in West Antarctica, where a full collapse of marine-based sectors during the LIG would lead to an additional ~1 m of GMSL rise or approximately 30% higher than values commonly cited in the literature. In our simulations, this contribution is largely established within ~1 ka or less (dashed-dotted line, Fig. 1B) of the end of the melt event. However, we interpret this time scale as an upper bound, since our calculations do not include any continuing uplift of marine sectors in response to the prior glacial cycle. Even over periods in which a full collapse of marine-based sectors of WAIS does not occur, the results for the CISM and GRISLI collapse scenarios (Fig. 6B), as well as the projection forced by the RCP 8.5 scenario, demonstrate that the water outflux mechanism will be important. The magnitude and rate of water expulsion in these cases are dependent on the detailed geometry of marine exposure and the underlying mantle viscosity structure. Whether on century or millennial time scales, this additional contribution to GMSL, as well as the rapid time scale in which it is reached, is a consequence of the unique setting of WAIS: The ice sheet is grounded on bedrock that is largely below local sea level, and it is underlain by an anomalously hot, low-viscosity mantle and thin lithosphere.


Sea level calculations

The sea level calculations are based on a gravitationally self-consistent ice age sea level theory that accurately accounts for all deformational, gravitational, and rotational effects, as well as shoreline migration due to changes in the perimeter of grounded, marine-based ice sectors (5456). Our sea level results are generated using a finite volume numerical scheme (47) that solves for the response of an elastically compressible, Maxwell viscoelastic Earth model with 3D variability in mantle viscosity structure to an arbitrary surface mass loading history.

For this study, we use the same code as in (57) and a similar grid, but with a larger volume of regional refinement. It is constructed first as laterally uniform to honor all first- and second-order Preliminary Reference Earth Model (PREM) interfaces (58). It has 67 radial layers across the mantle and lithosphere and a total of 17 million grid points. The internodal spacing is variable: 12 to 15 km from the surface to the Mohorovičić Discontinuity (MOHO) at 24.4-km depth, ~25 km from the MOHO to 220 km, and ~50 km down to the core-mantle boundary. For a more precise mapping of viscosity and lithospheric thickness variations in the load near field, a single mid-edge nodal refinement is applied under a surface covering all of Antarctica plus a seaward margin bounded between 49°S and 57°S and extending all the way from the surface to the base of the mantle. This refinement doubles the resolution described above, resulting in 133 radial layers under the Antarctic region and a total of 27 million grid nodes globally.

Construction of the 3D viscoelastic Earth models

To model instantaneous elastic deformation of the solid Earth, we adopt bulk and shear moduli from PREM (58). As realistic lateral variations in the elastic moduli have an insignificant effect on GIA, we use this purely 1D elastic structure in all simulations. In contrast, lateral variations in viscosity exert a strong influence on GIA, and we have constructed two individual 3D Earth models to illustrate this sensitivity. We first define a radially averaged viscosity structure, which comprises an elastic lithosphere that is 96 km thick, an upper mantle viscosity of 1 × 1020 Pa s, and a lower mantle viscosity of 5 × 1021 Pa s that extends from 670 km down to the core-mantle boundary.

The primary 3D Earth model used in the study (V3DSD; Fig. 2B) is described in the main text. To construct a second model (V3DRH; Fig. 2E), we use the SL2013sv seismic tomography model (2) in the upper 400 km, which contains a large number of surface wave constraints that are particularly sensitive to structure in the uppermost mantle. VS is first converted into temperature, T, using an inverse calibration scheme described in (59) that includes the effects of anelasticity observed in laboratory deformation experiments. Lithospheric thickness is taken to be the depth of the 1175οC isotherm (Fig. 2D) (4). Temperatures outside the lithosphere are converted into viscosity variations using the scheme in (60), which assumes that dislocation creep viscosity, η, is exponentially dependent on temperature according toη=η0exp[H(TT0)R T0T]where η0 is the background radial viscosity profile described above, H = 682 kJ/mol is a constant activation enthalpy, R is the molar gas constant, and T0 is the radially averaged mantle temperature outside the lithosphere (which closely follows the 1333οC adiabat). This conversion yields mean upper mantle viscosities beneath the southern Antarctic Peninsula that are consistent with the inference in (49) and are approximately half an order of magnitude more viscous than the results in (1) beneath the Amundsen Sea Embayment.

Beneath 400 km, we adopt the whole mantle SEMUCB-WM1 tomography model (61). VS is first converted into temperature by assuming a pyrolite composition and using the Perple_X Gibbs free energy minimization software (62) and a thermodynamic database (63) to obtain the elastic moduli. We correct for anelastic effects using the Q5 radial attenuation profile (64) and an appropriate lower mantle solidus (65). The resulting temperature variations are converted into viscosity variations using the diffusion creep activation enthalpy profile in (59) between 400 and 660 km, followed by the lower mantle profile in (66).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This study was supported by Fonds de Recherche du Québec–Nature et technologies (to L.P.), Star-Friedman Challenge (to L.P. and K.L.), NSF Graduate Research Fellowship Grants DGE1144152 and DGE1745303 (to E.M.P.), John D. and Catherine T. MacArthur Foundation (to J.X.M.), NASA grant NNX17AE17G (to J.X.M., E.M.P., and M.J.H.), NSF grant OCE-1702684 (to J.X.M.), American Chemical Society Petroleum Research Fund grant 59062-DNI8 (to M.J.H.), Harvard University (to L.P., E.M.P., K.L., J.X.M., and M.J.H.), Oregon State University (to J.R.C.), Natural Sciences and Engineering Research Council (to N.G.), and Canada Research Chairs program (to N.G.). Author contributions: L.P., E.M.P., K.L., and J.X.M. led the study with significant input from all the remaining authors. Revisiting calculations of GMSL emerged from their discussions with J.R.C. and P.U.C. K.L. worked with E.M.P. and N.G. to construct the 3D viscoelastic Earth model described here and with M.J.H. on the model described in Materials and Methods. K.L., in collaboration with L.P. and E.M.P., performed the sea level simulations, and they worked with J.X.M. to perform benchmark calculations on 1D Earth models. L.P., E.M.P., and J.X.M. wrote the original text, which was revised with input from all other coauthors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials and are available at Codes adopted here have been used in numerous previous studies. The Fortran code necessary to reconstruct sea level fingerprints on elastic Earth models (as in Fig. 4A) is available at the same link. Calculations of GMSL using any alternate AIS history or 3D viscoelastic Earth model can be run by and with the assistance of K.L., the developer of the code. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article