Observed rapid bedrock uplift in Amundsen Sea Embayment promotes ice-sheet stability

See allHide authors and affiliations

Science  22 Jun 2018:
Vol. 360, Issue 6395, pp. 1335-1339
DOI: 10.1126/science.aao1447

A quick rebound for Antarctic crust

Earth's crust deforms under the load of glaciers and ice sheets. When these masses are removed, the crust rebounds at a time scale determined by the viscosity of the upper mantle. Using GPS, Barletta et al. found that the viscosity of the mantle under the West Antarctic Ice Sheet is much lower than expected. This means that as ice is lost, the crust rebounds much faster than previously expected. Although estimates of total ice loss have to be revised upward, the surprising finding indicates that the ice sheet may stabilize against catastrophic collapse.

Science, this issue p. 1335


The marine portion of the West Antarctic Ice Sheet (WAIS) in the Amundsen Sea Embayment (ASE) accounts for one-fourth of the cryospheric contribution to global sea-level rise and is vulnerable to catastrophic collapse. The bedrock response to ice mass loss, glacial isostatic adjustment (GIA), was thought to occur on a time scale of 10,000 years. We used new GPS measurements, which show a rapid (41 millimeters per year) uplift of the ASE, to estimate the viscosity of the mantle underneath. We found a much lower viscosity (4 × 1018 pascal-second) than global average, and this shortens the GIA response time scale to decades up to a century. Our finding requires an upward revision of ice mass loss from gravity data of 10% and increases the potential stability of the WAIS against catastrophic collapse.

The viscosity of Earth’s upper mantle has a key role in Earth’s response to surface mass variations, especially visible in glaciated areas, at both local and global spatial scales. Those glaciated areas, including Greenland and Antarctica, contributed over the period 2002 to 2010 to global mean sea-level rise at an estimated rate of +1.22 ± 0.17 mm/year (1). Over the same time interval, the ASE, with an area of less than 4% of the Antarctic Ice Sheet (AIS), contributed to global mean sea-level rise at a rate of +0.3 ± 0.02 mm/year (1, 2), a fourth of the entire cryosphere contribution. Ice flow rates in the ASE region are increasing (3, 4), and the grounding line, controlled by the marine ice-sheet instability (5, 6) and subglacial geology (7), is retreating rapidly (8). Bedrock deforms due to a combination of Earth’s instantaneous elastic response to contemporary ice changes and the delayed viscoelastic response to past deglaciation (fig. S1). Regions underlain by high-viscosity mantle are mostly sensitive to the ancient ice history, whereas in the presence of low viscosity, only the recent ice history is relevant. The response time of ice‐related Earth deformations is determined by the mantle viscosity: Viscosities ~1018 to 1019 Pa s correspond to decadal to centennial time scales (912), and viscosities ~1020 to 1021 Pa s to millennial time scales (13, 14). GIA models for the whole of Antarctica adopt upper mantle (UM) (down to 670 km) viscosities larger than 1020 Pa s and up to 1021 Pa s driven by deglaciation starting 21,000 years ago at the Last Glacial Maximum. One of the most used models, ICE6G (13), is characterized by a UM viscosity with a value of 5 × 1020 Pa s, an intermediate value among the range used by most GIA models, and it predicts the largest present-day uplift rate in ASE among the models.

In glaciated regions undergoing large ice mass losses, such as the ASE, the GIA contribution has a considerable impact on gravity-derived ice mass variation estimates, because the Earth response, driven by the upper mantle viscosity, provides a contribution to the gravity field that, if not corrected for, translates into an apparent mass increase (fig. S2). The upper mantle viscosity has also been found to be crucial in predictions of ice-sheet stability, where bedrock uplift and consequent sea-level fall near the grounding line (fig. S3) have been shown to be an important stabilizing factor in ice-sheet behavior (15), with lower viscosity and consequent more rapid uplift resulting in greater negative feedback (16, 17). The impact of mantle viscosity on ice-sheet stability has been investigated (16, 17) for a wide range of viscosities (1019 to 1021 Pa s) that could not be narrowed further due to of lack of direct constraints on absolute mantle viscosity under the AIS. This issue is especially pronounced in ASE where ongoing ice mass loss is the largest. In fact, relative mantle viscosity variations under the AIS have been inferred based on seismic tomography (1820) but with poor constraints on absolute viscosity. Absolute viscosity can instead be constrained when it is dynamically derived from bedrock uplift rate or relative sea-level history. Finding low viscosity under sensitive areas of the marine ice sheet, such as ASE, would therefore have very important consequences on studies for estimating future ice-sheet stability and substantially reduce uncertainties regarding future ice mass loss projections. The possibility that the mantle viscosity in the ASE deviates significantly from the global average is supported by indications of geologically recent volcanism (21) and Cenozoic rifting in much of West Antarctica and limited evidence pointing to rift-related displacements of Oligocene age between the Thurston Island–Eights Coast and Marie Byrd Land crustal blocks flanking the Amundsen sector (22, 23). In other glaciated regions (Patagonia, Alaska, Iceland, and North Antarctic Peninsula), large uplift rates measured by dense local GPS networks have been used to constrain low viscosity (912). For this reason, we installed GPS on bedrock around the ASE to better constrain the viscosity of the shallow mantle.

We analyzed the deformation rates (Fig. 1) using six GPS stations as part of a global geodetic analysis following a recently developed processing strategy (24, 25). We found that the data indicate that ASE is undergoing one of the fastest GIA-related uplift rates ever recorded (up to 41 mm/year). In this work, we made use of solid Earth viscoelastic modeling to invert these data for the first estimate of absolute viscosity under ASE, using the new GPS solutions as constraints. To model uplift, we used a high-resolution approach (12, 26) with a spatial resolution of 1 km, which is much higher than previous studies (911). Our compressible viscoelastic Earth model has a viscosity structure more complex than the single-layer mantle assumed in previous studies (911) and provides insight into the rheology of the whole upper mantle down to 670 km. The paucity of GPS constraints and the scarcity of glaciological information about local ice history forced us to thoroughly explore the impact of those uncertainties on the viscosity model. To accomplish this, we designed a new strategy for incorporating ice history that allows us to combine several independent components (denoted as Hx, with x being an index) with different weights. The most important ice history components are the well-observed, high-resolution, regional, yearly surface distribution of present-day ice changes (27), which have the spatial trend shown in Fig. 1A, and a total mass trend of –130 Gt/year, for the period 2002 to 2014 (denoted as H0), and the component representing ice history for the period 1900 to 2002, constructed by rescaling the spatial pattern in Fig. 1A (denoted as H1) based on increased ice flux since 1970 (4). We performed a grid search over physical parameters by comparing the modeled predictions with the observations (27). The inversion parameters for the Earth model are the thickness of the elastic lithosphere (LT) (ranging from 20 to 90 km) and the viscosities of three layers: shallow upper mantle (SUM) (from the base of the lithosphere to 200 km depth), the deeper upper mantle (DUM) (from 200 to 400 km depth), and the transition zone (TZ) (from 400 to 670 km). We computed 807 different Earth models with about 10,000 variations of those models by linearly combining our results for different scenarios of ice history (27) (Fig. 1B).

Fig. 1 Observations of bedrock motions and ice changes.

(A) The average trend of ice mass changes in m/year of water equivalent for the 13 years 2002 to 2014 derived from satellite altimetry. (B) The cumulative total mass in Gt for the extrapolation of the ice history of the last century derived by rescaling of the pattern in (A). The black line represents the cumulative mass loss estimated from the yearly grids observed for the period 2002 to 2014, that is, the H0 ice history contribution. The blue solid line represents the cumulative mass loss for ice history contribution H1 for the period 1900 to 2002. We used the sum of the contribution for H0 and H1 as our initial estimate of mass change history. The black dashed line represents the ice history contribution obtained by extending the present-day rate to the past. The light blue and the red lines represent half and twice the rate of the H1 ice history contribution respectively. (C and D) show vertical (squares with error bars) and horizontal (arrows) observed GPS bedrock motion (black), respectively (reported in table S2). Horizontal motions are expressed in a frame that fixes, or nearly fixes, East Antarctica (27). Blue arrows show predicted elastic velocities (as in column 5 of table S2) using the observed ice mass changes shown in (A). Green arrows show the residual motions (observed uplift minus elastic predictions). Purple arrows—only for the vertical (C)—show the uplift prediction of the ICE6G (13) model. The ellipses in (D) represent uncertainties (±1 SD) for the associated observed horizontal velocity vectors. The background of (D) is the ice mass change in m/year of (A).

We first verified that the instantaneous elastic contribution to the uplift rate caused by the present-day ice mass changes in the ASE region are only ~20% of the recorded signal at best (Fig. 1, C and D), in marked contrast to the generally elastically driven high-uplift rates observed in Greenland (24). The spatial pattern of the high uplift rates at the ASE excludes the possibility of an alternative geological origin, such as volcanic activity (27), and is also incompatible with Earth deformation models driven by deglaciation starting at the Last Glacial Maximum 21,000 years ago. Fig. 1C shows the predictions from model ICE6G (13), which are the largest among recent GIA models, in which a stiff mantle viscosity (>1020 Pa s) produces a smooth, uniform deformational response at the spatial scale of the ASE. This is in sharp contrast to the more spatially variable uplift that we measured with GPS (Fig. 1C), and it is evident in the residual signal (simply residual, in the following) obtained subtracting the modeled elastic uplift rates from those measured by the GPS. Horizontal GPS displacements, which are heterogeneous in both magnitude and direction (Fig. 1D), show similar spatial behavior. The horizontal motions depend strongly on the location of ice mass changes, and the residual motion vectors have directions that are close to those of the elastic response to modern-day ice changes (Fig. 1D). This indicates that these large GPS rates (Fig. 1) are a response to ice mass loss with a pattern that is consistent with present-day ice changes (27).

By comparing the residual with predictions from our viscoelastic modeling, we found that a viscosity range for the SUM between 1018 and 6.3 × 1018 Pa s (log1018.0 to log1018.8) is suitable to explain the deformation rates observed in ASE (27). Unexpectedly, we also find a better-fit viscosity model by decreasing the viscosity of deeper mantle layers (note the range of low viscosity in the axis of Fig. 2). The viscosities of SUM and DUM are not independent (blue line in Fig. 2A): that is, if we chose to fix DUM viscosity at excessively high values, a good fit is obtained for excessively low values for SUM, but still worse than our best fit to observations. Thus, the observations collected by our GPS network configuration are indeed sensitive to both shallow and the deeper mantle structure (27). Even more surprisingly, we found that the resolving power of the GPS network extends to the TZ (Fig. 2B). The relationship we found between DUM and TZ (blue line in Fig. 2B) shows that the TZ layer has a minor but detectable influence on observed bedrock motion (27).

Fig. 2 Viscosity of the mantle to 670 km.

(A) (Left) The misfit χ (color scale from 2 to 16) as a function of the shallow upper mantle (SUM) viscosity (x axis) and the deep upper mantle (DUM) viscosity (y axis) for transition zone (TZ) viscosity fixed at 7 × 1019 Pa s, for experiment E2 (27). Estimates of lithospheric thickness, ranging from 30 to 70 km, are shown by circle size (right), with solid circles showing the acceptable estimates within the 95% confidence interval. The blue line represents the best-fitting relation between the DUM and SUM viscosities. The yellow line divides the region where the SUM has a higher viscosity than the DUM, a configuration considered less likely. The best fit for this experiment is for LT 60 km, SUM log10 18.6 Pa s, DUM log10 19.0 Pa s (TZ log10 19.84 Pa s) with a χ = 2.36. (B) (Left) The misfit χ (color scale from 2 to 9) as a function of the DUM viscosity (x axis) and the TZ viscosity (y axis) for three values of SUM viscosity (squares) and for fixed lithospheric thickness of 60 km for experiment E3 (27). (Right) Confidence interval, similar to the right plot in (A), applied to experiment E3 (27). The best fit for this experiment is for LT 60 km, SUM log10 18.6 Pa s, DUM log10 19.2 Pa s, and TZ log10 19.4 Pa s with a χ = 2.32.

We strongly constrained the upper mantle to a narrow range of low viscosity with our grid search (table S4 and Fig. 2); even considering uncertainties due to ice history reconstruction (27), fitting the observed signal with models results in R2 = 0.96 (Fig. 3A). Our preferred model [lowest root mean square deviation (RMSD)], using the initial estimate for ice history (H0 + H1 in Fig. 1B, using a mass loss rate of 130 Gt/year for the period 2002 to 2014 and of 32.5 Gt/year for the period 1900 to 2002), has a 60-km-thick lithosphere, a shallow upper mantle of ~4 × 1018 (log1018.6) Pa s; a deeper upper mantle that is about four times stiffer, ~1.6 × 1019 (log1019.2) Pa s; and a transition zone value of 2.5 × 1019 (log1019.4) Pa s (Fig. 3A, Best 2). These values provide a good fit to the vertical and horizontal nonelastic signal (Fig. 3A) that is within the confidence intervals for the sites with the highest signal-to-noise ratio. In all cases, we found a model preference for an elastic lithosphere that is thicker than 50 km, with a clear, rapid deterioration of the fit for a thinner lithosphere (Fig. 2A). The best fits for elastic lithosphere thicknesses are between 50 and 60 km, which are consistent with seismological studies of the region (18).

Fig. 3 Vertical and horizontal bedrock motion predictions.

(A) The vertical (left) and horizontal (right) observed GPS minus elastic modeled residual velocities (green dot in left plot, green arrows in right plot) as in column 5 of table S2. Vertical and horizontal velocity predictions from three viscosity models providing good fit (R2 = 0.96) are shown by the blue (Best 1), purple (Best 2), and dark green (Best 3) symbols. These models are our three best-fitting models (the total χ is reported in the legend of the right plot). In the text ,we indicate Best 2 as our preferred model because it has the lowest RMSD for the uplift rates. The error bars (left) are the weighted errors (error bar = SD/weight). The ellipses in (A) (right) represent uncertainties (±1 SD) for the associated residual horizontal velocity vectors. The weights are 1 except for LPLY and the horizontal for TOMO. Assuming that the current rate of ice thinning will be constant in the ASE throughout the next century, the evolution of uplift rates (B) and uplift (C) at the six GPS sites and points A and A′ of the profile along Pine Island Glacier are shown in (D). Ice mass loss in Pine Island Glacier (D, inset), and the elevation of the profile A-A′ for bedrock and ice [black and blue line from BEDMAP2 (32), respectively]. Bedrock variations along the profile AA′ from 2004 to 2114 with respect to 2014 are labeled.

We found that ASE has a rheological structure very different from those commonly used in global and regional dynamic Earth models (13, 14), with the upper mantle viscosity being about two orders of magnitude lower. Low viscosity provides independent support for a hot mantle, as suggested by seismological results (19) and an elevated geothermal heat flux (28) in the region. Low viscosity is also supported by the characteristics of the subglacial geology that provide important boundary conditions for the evolution of the ice sheet (7, 29). The low viscosity we found in deeper layers could be an indication of a plume-like thermal structure such as the one hypothesized in nearby Marie Byrd Land (19, 30). Yet it can more confidently be attributed to a mantle chemically altered by hydrous fluids and volatile silicate melts generated by protracted Paleozoic-Cretaceous subduction (31). Low viscosity could characterize other areas of West Antarctica where geological and seismic evidence similarly suggests that the presence of a hot mantle (19, 28, 30) is more widespread under areas where surface bedrock exposure is rare or absent. The viscosity we estimated under ASE is dynamically constrained without any a priori geological assumptions and can be used to set a reference value for a regional Earth model based on seismic tomography, which has large absolute uncertainties when interpreted for viscosity.

Observed ASE bedrock uplift rates, due to the low-viscosity mantle and short time scale Earth response to ice changes (Fig. 3A), are higher than expected (13, 14). The uplift rate is also predicted to accelerate with time (Fig. 3, B and C), producing an increasing impact on the gravity field. Both ice mass loss and the flow of mantle material perturb the time-variable gravity field observed by the Gravity Recovery and Climate Experiment (GRACE) satellite pair (fig. S2). The solid Earth effect results in an apparent mass increase that must be removed to estimate ice mass changes from observations of gravity changes. GRACE studies employ a constant rate of apparent mass change due to GIA of between +3.5 and +5.5 Gt/year water equivalent for the combined Pine Island and Thwaites glacier basins (2). Using our best local model, we calculated a much larger present-day apparent mass change of +16.7 Gt/year (27) for GRACE-derived mass trends over 2002 to 2014. This effect is, furthermore, not constant with time, as considered in previous studies (1, 2), because it depends on recent ice mass changes. We estimate that apparent mass change will increase by >2% per year in coming decades even if the ice mass loss rate does not increase. We thus contend that published GRACE-derived ice mass loss estimates for ASE, for example, ~108 Gt/year (1) (drainage basins of Pine Island and Thwaites and Smith glaciers), are systematically underestimated by between 10.0 and 13.9 Gt/year, which is more than 10% of the total ice mass loss estimate for the ASE.

The extremely low upper mantle viscosity that we constrain supports the possibility of increased stability of the WAIS with respect to previous studies (16, 17). Lower mantle viscosity leads to faster bedrock uplift in response to ice mass loss. Rapidly rising bedrock shallows the ocean at the grounding line and reduces the buoyancy forces experienced by the edge of the ice sheet while reducing the slope of the bed beneath the ice sheet (fig. S3). A laterally homogeneous viscosity of 1019 Pa s for the upper mantle has been used in the Earth modeling (16, 17) and found to be sufficient to prevent collapse of WAIS under moderate climate warming scenarios and to delay it for the more extreme scenarios. The timing of the Earth response is completely determined by the viscosity structure, and we numerically estimated (27) the response-speed increase when using our preferred model. Even under the assumption that the current rate of ice thinning will be constant in the ASE over the next century, we found that the uplift rate will increase (Fig. 3, B and C). In 100 years, uplift rates at the GPS sites will be between 2.5 and 3.5 times larger than currently observed (Fig. 3B), and the bedrock at the Pine Island Glacier grounding line will have risen by about 8 m compared to the present (Fig. 3, C and D). This is about three times higher than values shown to reduce runaway ice surface velocities within 100 years (15). The time required to build sufficient deformation to trigger the stabilization effect is much shorter (27) than in (16). Under low and medium climate forcing, with the onset of the stabilization feedback about two times faster (27), the condition for ice-sheet stability and its possible readvance can reasonably be expected to occur much earlier than predicted in previous studies (16, 17). Based on our estimates, it might produce a deformation large enough and early enough in the deglaciation phase to prevent the complete collapse of the WAIS even under strong climate forcing. To investigate this last point, a full ice-Earth coupled simulation as in (16, 17) that incorporates revised Earth properties is necessary. Earth structure can be expected to become increasingly relevant to the next generation of ice-Earth coupled models, where lateral variations under WAIS will be incorporated.

The weak Earth structure that we found under ASE and the related rapid stability feedback have a strong direct impact on WAIS evolution at the centennial time scale. The presence of such a low viscosity under ASE has effects on the millennial time scale also, with important implications on understanding of the potential WAIS collapse scenarios, as well as on the development of paleo ice-sheet history reconstructions and the associated long-term Earth deformation models (13, 14), hopefully helping to improve the accuracy of future climate change prediction and sea-level projection.

Supplementary Materials

Supplementary Text

Figs. S1 to S13

Tables S1 to S4

References (3385)

References and Notes

  1. Materials and methods are available as supplementary materials.
Acknowledgments: V.R.B. is indebted to C. Madison for important support in writing the manuscript. The authors thank DTU Computing Center (DCC), where the modeling simulations were performed, for the precious support provided. This research is a contribution to the SCAR SERCE program. Funding: This research was supported by the European Space Agency (ESA) in the framework of the project GOCE+ Antarctica - Dynamic Antarctic Lithosphere, the U.S. National Science Foundation Office of Polar Programs, and the U.S. Antarctic Program. VE-CL0V3RS was developed by V.R.B and A. Bordoni in a self-funded project. Author contributions: V.R.B. and M.B. planned the modeling strategy. V.R.B. computed the solid Earth predictions, performed the comparative analysis with the data, interpreted the results, prepared the figures, and wrote the original manuscript. M.B. provided GPS analysis, reference frame realization, trajectory analysis, interpretation, and writing. B.E.S. processed the altimetry data. T.W. supervised the work and provided input on the analysis, interpretation, and writing. A. Bordoni computed the Earth model Green’s functions and advised on analysis, interpretation, and writing. S.A.K. computed elastic deformation for a double-blind test and contributed to writing. M.R.-N. provided ice history uncertainty assessment. M.W. advised on altimetry data and contributed to writing. A. Brown and E.K. performed GPS data analysis. D.J.C., I.D., E.K., S.K., R.S., M.W., and T.W. performed GPS fieldwork. D.J.C. provided GPS instrumental support. S.K. advised on the comparative analysis. D.A.W., R.C.A., and A.N. provided seismological Earth parameters and contributed to interpretation. T.W. and I.D. worked on tectonic framework issues. Competing interests: None declared. Data and materials availability: The POLENET GPS data are open and immediately available from the Data Archive Interface (DAI, of UNAVCO. The grids of ice mass changes derived with satellite altimetry are available at Other relevant data are available at Code availability: The code VE-HResV2 computes the vertical and horizontal deformation given an input grid load and given the elastic or the viscoelastic Green’s functions for an Earth model. It is available as Fortran code upon request to the author. The code VE-CL0V3RS v.3.5.6 computes the elastic and the viscoelastic Green’s functions given an Earth model. The part of the code for the elastic Green’s function (or elastic Love numbers) is available as executable upon request to the authors. The part of the code for the viscoelastic Green’s function has been ported to high-performance clusters. It is still under development, and it is not user friendly at this stage. It is recommended that you ask the author to have it configured and run for you with your preferred Earth model. Both codes already have been used in other studies (12).

Stay Connected to Science

Navigate This Article