Damage accelerates ice shelf instability and mass loss in Amundsen Sea Embayment.

Pine Island Glacier and Thwaites Glacier in the Amundsen Sea Embayment are among the fastest changing outlet glaciers in West Antarctica with large consequences for global sea level. Yet, assessing how much and how fast both glaciers will weaken if these changes continue remains a major uncertainty as many of the processes that control their ice shelf weakening and grounding line retreat are not well understood. Here, we combine multisource satellite imagery with modeling to uncover the rapid development of damage areas in the shear zones of Pine Island and Thwaites ice shelves. These damage areas consist of highly crevassed areas and open fractures and are first signs that the shear zones of both ice shelves have structurally weakened over the past decade. Idealized model results reveal moreover that the damage initiates a feedback process where initial ice shelf weakening triggers the development of damage in their shear zones, which results in further speedup, shearing, and weakening, hence promoting additional damage development. This damage feedback potentially preconditions these ice shelves for disintegration and enhances grounding line retreat. The results of this study suggest that damage feedback processes are key to future ice shelf stability, grounding line retreat, and sea level contributions from Antarctica. Moreover, they underline the need for incorporating these feedback processes, which are currently not accounted for in most ice sheet models, to improve sea level rise projections.

Pine Island Glacier (PIG) and Thwaites Glacier (TG) in the Amundsen Sea Embayment are responsible for the largest contribution of Antarctica to global sea level rise (i.e., ∼5% of global sea level rise) (1, 2). Both glaciers show distinct changes in recent decades driven by changing atmospheric and oceanic conditions that cause enhanced ocean-induced melting of their floating ice shelves (3–5). Due to this enhanced melting, PIG and TG calving fronts retreated (6, 7) and their ice shelves thinned (8), decreasing the buttressing effect they exert on the upstream glaciers. As a result, both glaciers have accelerated and thinned and their grounding lines have retreated (9, 10). Under these conditions and in combination with a retrograde bed, PIG and TG are considered prone to marine ice sheet instability with the potential loss of their ice shelves and with large consequences for sea level rise (11, 12). Yet, quantifying the future timing and magnitude of these instabilities remains difficult as many of the key processes and their boundary conditions are poorly known or not accounted for in ice sheet models (13, 14).

Damage Evolution

In this study, we use time series of satellite imagery to show the rapid development of damage areas on the PIG and TG ice shelves (Fig. 1) that appear crucial to the future of both glaciers, but are not included in current models of their retreat. These damage areas consist of highly crevassed areas near the grounding line and of open fractures that contain dense ice mélange within the shear zones of both ice shelves. Satellite observations over the past two decades show the evolution from lack of crevasses in 1997 to rapidly growing crevasse-damaged areas near the grounding line and in shear zones on both ice shelves in 2019 (Fig. 1).

Fig. 1.

Damage evolution in Amundsen Sea Embayment. (A) Sentinel-2 satellite overview of the Amundsen Sea Embayment area showing PIG and TG glaciers and maximum strain rate (in purple-green) since 2015 derived from Sentinel-1 velocity time series. Zoom boxes, transects [P1, T1], and grounding line evolution from ref. 40 are illustrated in black and spectral colors (1992, purple; 1994, blue; 2000, green; 2011, orange), respectively. (B–D) Evolution of damage areas in Landsat (1997 to 2018) satellite image time series.

For PIG, this damage evolution started near the grounding line in 1999 as has been previously documented (7), but satellite imagery in our study shows how the initial damage has rapidly evolved since 2016 into tearing apart of the southern shear zone of the PIG ice shelf (Movies S1 and S2), whereas the northern shear zone remained largely intact after the unprecedented retreat and disconnection from the northern PIG ice shelf in 2015 (6). For TG, the damage started with the gradual disintegration of the shear zone between its glacier tongue and the eastern ice shelf and the subsequent removal of a large part of the TG glacier tongue as described by ref. 7. Since 2016, however, this TG damage moved farther upstream in the remaining shear zone between the glacier tongue and eastern ice shelf and evolved toward the rapid development of open fractures near the grounding line (Movies S1 and S3).

The observed damage for both PIG and TG ice shelves occurs typically in the shear zones where the ice shelf is thin (Fig. 2 ) and where high positive maximum strain rates occur (Fig. 1). These high maximum strain rates promote the development of damage through the opening of crevasses and rifts (Movies S1 and S3). This process can be clearly seen in the southern shear zone of PIG and in the shear zone between the TG’s glacier tongue and eastern ice shelf, where high maximum strain rates are observed (Fig. 1). In the northern shear zone of PIG, on the other hand, the observed damage evolution is absent or limited due negative maximum strain rates (Fig. 1) that result in closing of crevasses and rifts.

Fig. 2.

Elevation and thinning rates over PIG area. (A–E) ASTER satellite digital elevation model from January and February 2018 showing the elevation of PIG (A, C, and D) and TG (B and E) ice shelves and the damage areas in the thinner areas. Gray values are masked no data values. Zoom boxes, transects [P1, T1], and grounding line evolution since 1992 (40) are illustrated in black and spectral colors, respectively (see Fig. 1). (F–J) Elevation changes (2010 to 2017) derived from Cryosat-2 satellite altimetry data showing thinning over the PIG and TG glaciers and ice shelves in combination with the local advection of patches of thicker ice. (H–J) Zoom areas corresponding to boxes in F and G. (K) Velocity transects along transect P1 derived from multisource satellite imagery since 1992, where the colors represent the different observation years and the line types show different datasets. (L) Same as K but along transect T1 (Fig. 1).

Simultaneously with the damage development satellite altimetry (Fig. 2 ) shows an average elevation lowering of 0.3 m/y between 2010 and 2017 for PIG and TG glaciers and ice shelves, which locally can reach up to 13 m as described by ref. 15. Additionally, the satellite altimetry data reveal local patches of positive elevation changes, which are the result of the advection of patches of thicker ice (Fig. 2) in an Eulerian reference frame (16). Also the observed velocity gradients across the PIG and TG glaciers have increased in these shear zones with increases up to 30% since 1992 as well as with widening of the shear zones (Fig. 2 ). The largest velocity increase occurred between 2000 and 2010 corresponding to a period of warmer ocean waters that initiated an ice-dynamical response (10, 16–18).

Damage Preconditioning for Ice Shelf Disintegration

We hypothesize that the combination of localized ice shelf thinning, enhanced velocity gradients, and the rapid development of the damage areas in the shear zones of PIG’s and TG’s ice shelves is a sign that these ice shelves are already preconditioned for further disintegration. On the one hand this preconditioning is the result of the fact that the integrity of both ice shelves is compromised as shown in the satellite imagery, similarly to the preconditioning of the shear zones of Larsen B (19–21) prior to its collapse as a result of hydrofracturing (22). An example of this preconditioning can be seen in the interaction of these damaged shear zones with the existing fractures (23–26) perpendicular to flow in the center of PIG’s ice tongue. When these rifts connect with the damaged shear zone, the ice front is no longer stabilized due to the structural weakening, resulting in large calving events. A similar condition happened in September and October 2018 and February 2020 when a large rift from the damage zone developed across the PIG ice shelf (Movie S4), resulting in an unprecedented retreat of the ice shelf front (Movie S5). Due to this retreat, the ice shelf is no longer stabilized by the southwestern tributary ice inflow (6), resulting in a further destabilization of PIG’s ice shelf.

On the other hand, the preconditioning is the result of a damage feedback process, where damage enhances speedup, shearing, and weakening, hence promoting additional damage development, but where loss of buttressing, enhanced shearing, and glacier speedup also enhances damage. Initially, this damage feedback can be triggered by a variety of weakening processes that range from 1) oceanic melting that undermines the shear zones by thinning the ice shelf from below (17, 27); 2) high-velocity gradients that result in surface troughs (28); and 3) ice shelf retreat that results in unpinning, reduced buttressing of existent pinning points (9, 29), grounding line retreat, and glacier speedup (30). Each of these weakening processes could result in the initial development of damage areas in the shear zones. Once the damage is initiated, however, the feedback process kicks in and the weakening of the shear zones results in further speedup and shearing, hence promoting additional damage development.

Modeling of Damage Feedback

To assess the importance of this damage feedback in the shear zones, a continuum damage model (CDM) was coupled to the BISICLES ice sheet model (31) in an idealized setup to illustrate the impact of damage on ice sheet response. Earlier simulations (31) with this model have already shown that inclusion of damage in the model corresponds to a flow enhancement factor of 10 or ∼20 K difference in ice temperature, whereas earlier work with another damage model (21, 32) has shown that it is possible to use damage models to assess the mechanical weakening of ice shelves. In this study, however, we opted for a different approach with the goal to quantify the potential importance of the damage feedback process relative to the process of ocean-induced melting, which is typically modeled for future ice sheet scenarios. Therefore, the model was applied on an idealized marine ice sheet geometry with retrograde bed slopes and strong lateral stresses, following the Marine Ice Sheet Model Intercomparison Project (MISMIP+) setup, which has strong similarities to PIG conditions (Fig. 3 and ). We deliberately chose such an idealized setup as it gives a better experimental control, and a real-world experiment may suffer from the chosen initial conditions that—in this particular case—are poorly constrained.

Fig. 3.

Ice sheet model with damage feedback result (see Movie S6 for animated version). (A) Initial ice thickness in BISICLES ice sheet model with initial grounding line (black dotted line). Grounding lines after 100 y model simulation for different model forcings are illustrated in colors corresponding to key in B. Dashed lines are used for no initial damage and solid lines for damage via introduction of crevasses at the location of star in C on both sides of the shear margin, while the red dotted-dashed line corresponds to channelized melting scenario. Black dashed-dotted lines illustrate location of profile in B. (B) Profile showing the geometries and grounding lines locations for different model forcing corresponding to different colors (see key) where dashed lines represent no initial damage, solid lines show runs for damage via introduction of crevasses at the location of the star in C, and the red dotted-dashed line corresponds to a channelized melting scenario. (C) Integrated damage for 25 m/y oceanic melt at grounding line with 20 m enhanced damage. Upper half shows simulation without the local damage enhancement at the grounding line, whereas Lower half shows simulation with damage enhancement at the grounding (at location of star) via introduction of crevasses. (D) Integrated damage for 25 m/y oceanic melt at grounding line with 20 m enhanced channelized melting in the lower shear zone. Gray values represent areas where the ice thickness is 0. (E) Maximum (Max) strain rate for simulation without/with (Upper/Lower half) enhanced damage at the grounding line. (F) Max strain rate for simulation with channelized melting. (G–L) Changes in damage fraction (G and H), thickness (I and J), and velocity (K and L) in 25 m/y melt run without/with 20 m vertically integrated crevasse depths implemented in the shear zone at grounding line (G, I, and K) and without/with channelized melting (H, J, and L), respectively, where negative/positive values represent thinning/thickening/increased values or slowdown/speedup/decreased values due to damage introduction and channelized melting, respectively.

In this idealized setup we carried out several time-dependent simulations with different ocean-induced melting and damage parameterizations to assess the importance of damage in the shear zones relative to ocean melting. In this framework, several damage scenarios were implemented by locally enhancing damage in the shear zone at the grounding line via the introduction of crevasses/damage in the model, where damage is expressed as vertically integrated crevasse depth in meters () or by the introduction of channelized melting that may weaken ice shelf margins as in ref. 28.

The model results in Fig. 3 and Movie S6 highlight the importance of the damage feedback as a driver for ice shelf instability as the initial weakening due to localized damage or channelized melting results in enhanced damage that further weakens these ice shelves. Locally enhancing damage at one point in the shear zones (Fig. 3) causes the increase and expansion of damage across the entire shear zone (Fig. 3) together with the increase of maximum strain rates (Fig. 3), resulting in the weakening of the complete shear zone and ice shelf farther downstream. Locally enhancing channelized melting in the shear zone has a similar effect with also the development of weaker shear zones (Fig. 3), enhanced maximum strain rates (Fig. 3), and an increase of overall damage fractions relative to the undisturbed model simulations.

The model results also indicate enhanced thinning at the grounding line and in the shear zones and even open gaps in the ice shelf that decouple the ice stream and could be filled with open water or mélange. This thinning and decoupling reduce ice shelf buttressing, leading to accelerated ice transport to the ocean and consequent thinning of grounded ice upstream. Moreover, the thinning and decoupling in the shear zones can result in the advection of patches of relatively undamaged, thick ice (Fig. 3 ). This modeled pattern of thinning in the shear zones and advection of patches of undamaged thicker ice corresponds to the elevation changes observed by satellite altimetry (Fig. 2) and the visual pattern of crevasses and rifts that resemble the strongly thinned or open patches (Movies S1 and S3). Finally, the model shows a speedup of the glacier tongue as a result of the weakening, which results in an increase in maximum strain rate (Fig. 3 ) and in velocity gradients across the shear zone (Fig. 3 ) similar to the increasing velocity gradients across PIG and TG observed by satellites (Fig. 2 ).

The damage also has an important impact on the modeled grounding line retreat as the enhanced damage scenarios in the model initiate an enhanced grounding line retreat. For both the local initialized damage and the channelized melting this grounding line retreat is equivalent to quadrupling the oceanic melting near the grounding line from 25 to 100 m/y (Fig. 3 ). This indicates that weakening the shear zones and enhancing damage are powerful mechanisms for triggering strong grounding line retreat as all enhanced damage model runs result in a stronger grounding line retreat (up to 200%) than the model runs with only strong oceanic melting and no enhanced crevassing; for example, a 20 m/y oceanic melting combined with crevasses of 20 m deep in the shear zone at the grounding line causes stronger grounding line retreat than 100 m/y oceanic melting without this enhanced damage.

Although the results of idealized model output show similarities with observed damage, thinning, and velocity evolution, it is important to stress that the idealized experiments do not allow us to directly evaluate the observed changes at PIG and TG. First, the idealized model may not include the potential feedbacks that might be important when interpreting the observations. For example, the rapid thinning near the grounding line as a result of reduced buttressing may result in a larger ice flux from upstream that reduces grounding line thinning and slows down grounding line retreat. Second, the observed changes in reality are the result of the interplay between oceanic forcing and bathymetry [e.g., the stagnation of grounding line retreat over PIG (10, 16) due to reduction of oceanic melting (5, 17) or the changes in individual pinning points (6, 29), sea water intrusions (15), or tidal/wave-induced stress (33, 34)]. The observed changes in damage, thinning, and velocity gradients are therefore not expected to be the result of damage only, but also include these other drivers including the abated ocean forcing since 2011 (17). In both the observations and the idealized experiment it is not possible to distinguish between the cause and effect of the damage feedback on the observed speedup and thinning, but our idealized experiments do, however, allow us to assess the importance of the weakened shear margins and damage process. This illustrates that weakening these glaciers at their most vulnerable locations, as can be currently observed in the satellite observations, is a very effective way of introducing grounding line retreat, increased ice flux, and hence mass loss. It is therefore crucial to take these weakening processes into account when modeling the evolution of PIG and TG.

Damage Implications and Conclusion

Our satellite results show that the initial damage as described in refs. 7 and 28 is now rapidly expanding and that the structural weakening of the Pine Island and Thwaites shear margins, which was simulated for Thwaites by ref. 11, is already extensively happening under current conditions. These quickly weakening shear zones precondition the PIG and TG ice shelves for future area losses and possible collapse, similar to the preconditioning of the shear zones of Larsen B (19–21) prior to its collapse as a result of hydrofracturing (22). Although the potential of such a collapse through hydrofracturing in the Amundsen Embayment may be restricted due to the limited projected surface melt in this region (35), the damage makes the future response of PIG and TG ice shelves more sensitive to varying and extreme future atmospheric, oceanic, and sea ice conditions (4, 5, 36, 37) and could trigger a nonlinear response (17). The weakening could lead moreover to changes in calving patterns (6) and to changes in stabilizing pinning points and ice inflow (6), which could result in large reductions of PIG and TG ice shelf area. Nevertheless, even without a collapse, our damage model results imply that these damage areas close to the grounding lines could have important implications as the observed damage in these shear zones makes them vulnerable to enhanced mass loss and grounding line retreat.

In the future, this mechanical weakening and increased velocity gradients due to loss of frictional gradients at the ice shelf margins are not expected to trigger negative feedbacks that counterbalance the damage as damage healing is expected only for negative maximum strain rates, which are limited for ice shelves (38). As such, it is different from other ice shelf weakening processes such as surface or subshelf melt as these can be counterbalanced by changes in atmospheric or oceanic processes, which are prone to climatic variability (3–5). Therefore, the damage process and mechanical weakening in the shear zones have similar far-reaching consequences for ice shelf stability as localized ice shelf thinning in basal channels (27, 39). This sensitivity suggests that incorporating damage processes in future ice sheet models in combination with accurate knowledge of ocean forcing, bathymetry, bedrock topography, ice velocity, and surface melt is crucial to assess the future sea level contributions from major Antarctic glaciers.

Materials and Methods

Satellite Imagery of Damage Evolution.

Multisource satellite imagery was used to show the development of the damage areas. First, annual optical satellite composites from Landsat and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) (panchromatic/blue band with spatial resolution of 15/30 m for images after/of 1997, respectively; downloaded from the US Geological Survey Earth Explorer https://earthexplorer.usgs.gov) time series were processed to show the development of the damage areas in polar stereographic projection (Fig. 1 and Movies S2 and S3). Second, a time series of Copernicus Sentinel-1 synthetic aperture radar (SAR) backscatter imagery (ground range detected [GRD] images of interferometric wide [IW] with spatial resolution of 10 m) was processed to show the recent development of damage areas at PIG and TG since 2014 including all images during the austral winter (Movies S1 and S4).

Surface Elevation and Elevation Changes.

A surface elevation dataset was retrieved by mosaicking individual ASTER digital elevation models derived from AST-L1A data acquired in January and February 2018 (freely available at https://search.earthdata.nasa.gov/) using the open-source Ames Stereo Pipeline (ASP) (41, 42) (Fig. 2).

The elevation dataset was complemented with a map of elevation change at 500 m horizontal resolution. Rates of elevation change for 2010 to 2017 were derived from interferometric measurements from Cryosat-2. European Space Agency L1b waveforms retrieved by the satellite in its synthetic aperture radar interferometry (SARIn) mode were processed following the swath processing approach of Gray et al. (43). The resulting dense set of time-dependent elevation measurements was then used to derive elevation changes in a Eulerian framework at a 500-m resolution following the method presented in Wouters et al. (44).

Ice Velocity, Strain Rates, and Grounding Lines.

Time series of ice velocity data were retrieved by combining different available velocity datasets: 1) ice velocity data from feature tracking Copernicus Sentinel-1 since 2014 (45), which are available at 200-m resolution via the Enveo Cryoportal website (46) and which are shown as solid lines in Fig. 2 ; 2) ice velocity data from Making Earth System Data Records for Use in Research Environments (MEaSUREs) Annual Antarctic Ice Velocity Maps between 2005 and 2017 (47), which are available at 1-km resolution via National Snow and Ice Data Center (NSIDC)’s data portal (48) and which are shown as dotted lines in Fig. 2 ; 3) ice velocity data from MEaSUREs InSAR-Based Ice Velocity of the Amundsen Sea Embayment (9), which are available for 1996, 2000, 2002, and 2006 to 2012 at 1-km resolution via NSIDC’s data portal (49) and which are shown as dashed-dotted lines in Fig. 2 ; and 4) ice velocity data from European Remote sensing Satellite tandem interferometry/offset tracking, which are available at 500-m resolution via the Enveo Cryoportal website (50–52) and which are shown as dashed lines in Fig. 2 . All velocity data were subsequently averaged per dataset and year of acquisition. Subsequently, maximum strain rates were calculated by deriving the first principal stress component from strain rates derived from the two-dimensional velocity data (Fig. 1), where negative/positive maximum strain rates result in healing/enhancement of damage, respectively. Grounding line evolution (40) was derived from MEaSUREs Antarctic Grounding Line from Differential Satellite Radar Interferometry, Version 2 (53).

Ice Sheet Modeling.

The BISICLES ice sheet model with a CDM was implemented to assess the importance of the damage feedback [see Sun et al. (31) for a complete description]. Ice flow velocity in the BISICLES-CDM ice sheet model is computed by solving the vertically integrated stress balance equation, according to the shallow shelf approximation (SSA) where the temperature is assumed constant throughout the simulations. The damage continuum model considers the conservation of damage (54) due to downstream ice advection and local sources of damage (55, 56). This includes a time-dependent enhancement factor in Glen’s flow law (54–56) due to damage that affects the ice viscosity. In this way, both the development of damage and the ice flow field are strongly coupled. This allows carrying out idealized numerical experiments examining the interaction between damage and large-scale ice sheet and shelf dynamics (31). The model was run for 100 y using an adaptive mesh refinement framework with four levels from 0.5 to 2 km, which allows to use a nonuniform, evolving mesh during simulations. Damage was introduced by introducing vertically integrated crevasses in the model and by modifying Glen’s flow law based upon the crevasse opening formula of Nye (55). Damage is consequently expressed as vertically integrated crevasse depth in meters, which can be converted to a unitless fraction of the vertically integrated ice thickness () after dividing by the ice thickness. consequently takes values between 0, for fully intact ice, and 1, for ice that is cracked through its full extent similar to the isotropic scalar damage used by refs. 21 and 32.

The BISICLES-CDM model was applied on the MISMIP+ geometry (57, 58), where ice flows along an 800-km-long and 80-km-wide submarine bedrock trough, from an ice divide at one end to an ice shelf and calving front at the other and which show strong similarities with PIG conditions. We deliberately opted for such an idealized model setup to prevent shortcomings due to unknown initial conditions of the ice shelf that could obliterate the mechanisms at work. As such, the setup leads to a greater control over the experiment to delineate the impact of damage adjacent to ice shelf weakening due to subshelf melt. As such, it enables us to comprehend the physical mechanism at work in conjunction with observed features.

To quantify the impact of the damage feedback on ice sheet/shelf thickness, velocity, strain, and grounding line retreat, a set of numerical experiments was implemented where we applied varying ocean forcing rates ranging from 25 to 100 m/y and where we enhanced damage in the shear zone at the grounding line by the local introduction of crevasses or by the introduction of enhanced channelized melting of 20 m/y in the southern shear zone. The choice for the local introduction of damage is motivated by the fact that the observed damage at the grounding line (e.g., due localized melting) (15, 59) is not represented in the idealized setup of Sun et al. (31). In this framework we implemented different crevasse depths for different damage scenarios. These scenarios varied from no enhanced damage to 100-m vertically integrated crevasse depths. The location of initiated damage is illustrated in Fig. 3 by the black star and is applied on both sides of the shear margin and this location evolves with grounding line migration (i.e., if the grounding line retreats, the location of enhanced damage also retreats). These locations compare well to the locations of observed damage origin and correspond to the observations where the damage is constantly initiated locally close to the grounding line while the damaged zones are subsequently advected downstream (Movie S1). As such, the setup also differs from just advecting a single crevasse as we constantly reimplement a crevasse at the grounding line, resulting in a damaged area that expands while being advected downstream (Movie S6). The set of varying oceanic melting rates and enhanced damage values allows us to assess the interplay between ocean forcing and damage feedback and by analyzing the changes in ice sheet/shelf thickness, velocity, strain, and grounding line retreat we quantified their relative importance. Additionally, a simulation was performed where we combined the 25-m/y ocean forcing introduction of enhanced channelized melting of 20 m/y in the southern shear zone. The choice for this channelized melting scenario is motivated by the results of ref. 28 that show that the shear zone might also be weakened by channelized melting.