Identification of ice-shelf damage

Satellite images (Figs. 1 and 2) are ASAR WSM imagery from ENVISAT, which operated from March 2002–April 2012. Owing to its all-weather, day–night imaging capability, ENVISAT ASAR WSM imagery is available in late autumn and winter, which covers the Wilkins and Voyeykov calving events. A basemap of Antarctica, its islands and the Voyeykov and Wilkins ice shelves was imported from the Making Earth System Data Records for Use in Research Environments (MEaSUREs) Antarctic Boundaries for International Polar Years 2007–2009 from Satellite Radar (Version 2)49. We edited the Wilkins shelf fronts in the MEaSUREs Antarctic Boundaries dataset to match the ENVISAT ASAR WSM imagery. As the Voyeykov is unconfined and terminates into a widespread mélange layer, we edited its shelf-front locations in the MEaSUREs Antarctic Boundaries dataset to match those given by ref. 9. Fracture locations were identified using a method previously applied to the Voyeykov9 and the George VI50. The method uses both ENVISAT ASAR WSM and Landsat-7 Enhanced Thematic Mapper Plus images, with fractures manually mapped in Quantum Geographic Information System (QGIS). For guidance, we used similar work undertaken to map fracture locations for the Wilkins from 1990–2008, using imagery from ENVISAT ASAR, the Advanced Land Orbit satellite and Terra Advanced Spaceborne Thermal Emission and Reflection Radiometer satellites7,31.

The reported percentage changes in ice-shelf area due to calving were calculated by comparing the area of the ice shelf pre- and post-calving. This was achieved by comparing ice-shelf extents over different years using ENVISAT ASAR WSM imagery and manually drawing polygons in QGIS to represent the areas lost to calving. Extended Data Table 1 lists all ice-shelf images used in the study.

Sea-ice lengths

Reported sea-ice lengths (Fig. 3) were calculated from a combination of (1) sea-ice concentration data from the AMSR-E ASI sea-ice dataset51 and (2) fast-ice extent datasets from ref. 25. Sea-ice concentration data are given as a percentage of areal sea-ice cover per grid cell on the ocean surface at a spatial resolution of 3.125 km2 and a daily frequency from June 2002–September 2011. The fast-ice dataset is a binary measure per grid cell at a spatial resolution of 1 km2 and a 15-day frequency from March 2000–March 2018. The fast-ice dataset includes an Antarctic land mask, for which ice-shelf extents are updated every March25. For consistency, we changed the Voyeykov shelf front in the fast-ice dataset to that given by ref. 9, and we updated the Wilkins land mask at each calving event in 2008. To pair the sea-ice and fast-ice datasets, we converted the fast-ice dataset to a 3.125 km2 resolution (to match the AMSR-E grid). AMSR-E sea-ice outputs were imported into QGIS, which produced the sea-ice maps in Fig. 2, where the pack-ice boundaries are corrected for fast ice from the converted fast-ice dataset25.

An algorithm to calculate the presence and age of fast ice and the location of the shelf fronts was applied to the sea-ice and fast-ice datasets on the 3.125 km2 grid. For each shelf front:

(i) On the first day of the fast-ice record (1 September 2000), a counter was initiated to track the number of continuous days fast ice is present in each grid cell. (ii) Starting from 1 September 2002, a daily output containing the coordinates of ocean, land (that is, islands and the ice sheet), ice shelf, shelf front, annual fast ice and perennial fast ice was generated. The ice-shelf edge coordinates were calculated by checking for ice-shelf cells adjacent to ocean or fast-ice cells. (iii) Repeat (ii) for each day until 31 August 2009.

The following algorithm was then applied to calculate the sea-ice lengths. For each shelf front and each day in the record:

(a) Incorporate the coordinates generated by step (i) of the previous algorithm into the daily AMSR-E ASI sea-ice dataset. Define pack ice as any cell with a sea-ice concentration of ≥15%. In case of a conflict in which two different types of ice or land are present in a cell, priority is given to land, followed by perennial fast ice, annual fast ice and then pack ice. All remaining unclassified cells are then classified as ocean. (b) A pair of cells is chosen, where one cell is along the shelf front and the other cell is 10° N of the mean shelf-front latitude and between 12° W and 5° E of the mean shelf-front longitude (Wilkins RC and Voyeykov) or 12° W of the mean shelf front longitude and between the coastline and 10° N of the mean shelf-front latitude (Wilkins CL) (the 10° N boundary was chosen to avoid sea ice that would contaminate the swell data). A straight line ‘transect’ was calculated between these two cells (white and orange lines in Extended Data Fig. 2a,c). This leads to a higher density of sea-ice transects towards the centre of the shelf fronts (Extended Data Fig. 2a,c). (c) Sea-ice data are extracted for each cell along the transect. If land is present, the transect is discarded. (d) The effective pack-ice, annual fast-ice and perennial fast-ice lengths along the transect, \({L}_{{\rm{pck}}}^{(\mathrm{t})}\), \({L}_{{\rm{ann}}}^{(\mathrm{t})}\) and \({L}_{{\rm{per}}}^{(\mathrm{t})}\), respectively, are calculated. The effective pack-ice length is the sum of the cell lengths containing pack ice multiplied by the mean sea-ice concentration in those cells, whereas the effective fast-ice lengths are simply the sums of the cell lengths containing the relevant fast-ice type. (e) Repeat (b)–(d) for each possible pair of cells.

On a given day, the total effective sea-ice length for a particular transect is calculated as

$${L}_{{\rm{tot}}}^{(\mathrm{t})}\equiv {L}_{{\rm{pck}}}^{(\mathrm{t})}+{L}_{{\rm{ann}}}^{(\mathrm{t})}+{L}_{{\rm{per}}}^{(\mathrm{t})}.$$ (1)

For a chosen shelf front on a given day, the transects with the lowest third of \({L}_{{\rm{tot}}}^{(t)}\) were considered as the weakest part of the sea-ice barrier. The daily values reported for L pck , L ann and L per are the mean values corresponding to these transects (Fig. 3d–i). The daily L tot was obtained from the sum of L pck , L ann and L per . A total sea-ice length in which each of the transects used was first normalized by its total transect length was also used (\({\hat{L}}_{{\rm{tot}}}\); Fig. 4d–f). The number of low sea-ice days in each year for which L tot < 100 km, alongside days with L pck ≤ 50 km, are shown in Fig. 3a–c.

Incoming ocean swell

Peak periods (T p ) and significant wave heights (H s ) were extracted from the Collaboration for Australian Weather and Climate Research Wave (CAWCR) hindcast – aggregated collection52. The CAWCR hindcast was generated using the WaveWatch III v4.08 wave model, which is forced by hourly wind and uses daily sea-ice concentration data from the research data archive operated by the National Center for Atmospheric Research (NCAR) for pre-2011 data. The hindcast is from 1 January 1979 to the present, with hourly outputs on a 0.4° spatial grid of the ocean surface, which were converted to daily maximum values for H s and T p . The operational range extends to 78° S. We used swell data from the hindcast for sea-ice concentrations of ≤25%, as it is unaffected by the presence of sea-ice53.

For each day over the 7-year record and each shelf front, incoming H s and T p were calculated from the hindcast as follows:

(I) Obtain H s and T p for the hindcast cells contained within the bounds defined in (b) from the above algorithm (Extended Data Fig. 2a,c). (II) Retain the median H s and T p values.

This provided daily H s and T p time series for each shelf front (Extended Data Fig. 3).

Model of swell-induced shelf-front flexural stress

The model used to predict swell-induced flexural stress on the Wilkins and Voyeykov shelf fronts (Fig. 4) combines models of swell attenuation over distance due to sea-ice barriers (pack, annual fast and perennial fast)28 with a model of ice-shelf flexure in response to incoming swell30 (Extended Data Fig. 4). The models are two-dimensional (a horizontal ‘wave’ dimension and a depth dimension), time-harmonic and involve interactions between regular incident ocean waves of a prescribed amplitude (A inc ), wave period (T) and a sequence of the different ice covers. They are based on linear potential-flow theory for the water motions combined with thin-elastic-plate theory for the floating ice cover, where the Young’s modulus is set to be 6 GPa for the sea-ice barriers54 and 8 GPa for the ice shelves55.

A first model was used to predict the exponential attenuation over distance of the incident wave through the pack-ice region, such that the amplitude reduces to

$${A}_{{\rm{pck}}}\approx {A}_{{\rm{inc}}}\exp (-{\alpha }_{{\rm{pck}}}{L}_{{\rm{pck}}}/{C}_{{\rm{pck}}}),$$ (2)

where α pck (T) is the attenuation coefficient. The ice cover was modelled as an array of floes with a thickness of 0.52 m and lengths randomized about a mean value of 500 m (ref. 56) and at an average concentration C pck . Attenuation is produced by wave scattering at the floe edges and wave energy dissipation due to drag pressure57, calculated using the Robinson–Palmer model, which gives a frequency dependence similar to observations27 with a damping parameter of 13.5 Pa s m−1 (refs. 11,54,58). The evolution of wave attenuation models, which started in the 1970s, is documented in multiple review articles (for example, refs. 59,60,61). Sensitivity studies have also been conducted, such as that in ref. 62. Three-dimensional versions of the model are available but at a far greater computational cost and without evidence of major differences in the predicted attenuation coefficients63,64.

Annual and perennial fast-ice models were used to predict the reduction in the incident waves due to reflection at the fast-ice edges and the subsequent exponential attenuation of the waves over distance through the fast-ice regions65, such that the wave amplitude reaching the shelf front is

$${A}_{{\rm{per}}}={a}_{{\rm{per}}}{A}_{{\rm{ann}}}\exp (-{\alpha }_{{\rm{per}}}{L}_{{\rm{per}}})\quad \,\text{where}\,\quad {A}_{{\rm{ann}}}={a}_{{\rm{ann}}}{A}_{{\rm{pck}}}\exp (-{\alpha }_{{\rm{ann}}}{L}_{{\rm{ann}}}),$$ (3)

a • denotes the transmitted (non-reflected) components and α • are the attenuation coefficients. The models are based on the one outlined in the seminal study by Fox and Squire66, but include the Robinson–Palmer dissipation model for attenuation, as in ref. 67. The annual fast-ice thickness was set to 2 m (ref. 68) and the perennial fast-ice thickness was 11 m (refs. 18,68).

An ice-shelf model was used to predict the vertical displacement of the ice shelf, η(x), due to flexural-gravity waves, forced by the proportion of the incident wave that reaches the shelf front, such that

$$\eta (x)={A}_{{\rm{per}}}\hat{\eta }(x),$$ (4)

where \(\hat{\eta }(x:T\,)\) denotes the displacement due to unit-amplitude forcing. The form of the model used was pioneered in refs. 42,69, and has since been extended to include an Archimedean ice-shelf draught (90% of the ice thickness)70 and varying spatial profiles of ice-shelf thickness and seabed bathymetry29. The version of model used in this study has been validated in terms of the ratio of flexural-gravity-wave to incident-ocean-wave amplitudes observed on the Ross Ice Shelf30. The model also predicts that swell is far more sensitive to changes in ice-shelf thickness than seabed variations55.

The maximum flexural stress per wave period is calculated from the displacement, such that

$$\sigma =\left\vert \frac{EH}{2(1-{

u }^{2})}{\eta }^{{\prime\prime} }\right\vert ,$$ (5)

where E is Young’s modulus and ν = 0.3 is Poisson’s ratio. Relative to the forcing amplitude, flexural stress tends to increase monotonically with wave period in the swell regime, and is negatively correlated with the shelf-front thickness55. Owing to the underlying assumptions of plane stress and strain71, equation (5) is the only non-zero component of the stress tensor. Three-dimensional models that involve more non-zero components of the stress tensor are currently being developed72,73. However, they require numerical solution methods, which would be prohibitively expensive for the present study, and it is not yet clear in which regimes 2D model predictions of stress require the additional spatial dimension.

For each shelf front and each day over the 7-year record, the model was applied to L pck , C pck , L ann and L per , which were derived as described above. The incident wave amplitude and period are set to be A inc = 0.5H s and T = T p using the median significant wave height and peak period for the relevant shelf front and day, derived as described in the previous section.

Extended Data Fig. 2b,d shows model outputs corresponding to the highlighted Wilkins RC transects in Extended Data Fig. 2a,c. Exponential attenuation of swell through the sea-ice barriers is evident, along with amplitude drops at the fast-ice edges due to partial reflection of the incident waves28. The attenuation rate per kilometre is greatest in the pack ice due to the broken nature of the ice cover, and the attenuation rate is greater in annual fast ice than perennial fast ice. In contrast, the thicker perennial fast ice creates a greater amplitude drop than the annual fast ice. Once the incident wave reaches the shelf front, it imparts flexural stress as it travels through the ice shelf. The stress is responsive to variations in ice-shelf thickness, with a thickness decrease causing a stress increase, although the peak stress is away from the shelf front where the stress values are sampled.

For each shelf front and each day, the flexural stress outputs were reduced to a single value, σ, by taking the mean over 1–2 km from the shelf front for each transect and taking the mean over the transects. The interval of stress locations was chosen to be close to the shelf front but far enough away to avoid the shelf front boundary layer caused by the free-edge conditions. The accumulated stress (Fig. 4) is the weighted sum

$${\sigma }_{N}=\mathop{\sum }\limits_{n=0}^{N-1}\frac{(N-n){[\sigma ]}_{{\rm{day}}-n}}{N},$$ (6)

that is, the stress summed over an N-day window, where the weight of the stress decreases linearly with age until it moves outside the window.

Results in Fig. 4d–f were obtained using 60-day windows, which balance the responsiveness of the accumulated stress to changes in daily stress and capture build-ups over time (Extended Data Fig. 5). The normalized daily accumulated stress (\({\hat{\sigma }}_{60}\); Fig. 4d–f) is the ratio of the accumulated stress to the mean yearly maximum stress over the 7-year record for the shelf front. Corresponding results for climatological sea-ice lengths are also given (\({\hat{\sigma }}_{60}^{{\rm{aice}}}\)), which were calculated using a daily climatological average of L pck , L ann and L per and the same wave conditions as σ 60 .

Topological altimetry

The CryoSat-2 Synthetic Aperture Radar Interferometry Level-2I processing baseline-E mode was used to obtain surface elevation data to calculate ice-shelf thicknesses. The CryoSat-2 satellite was launched on 8 April 2010, and provides coverage to 88° S on a 369-day repeat cycle. For each ice shelf, all available tracks were taken from 2014–2023 within the ice-shelf boundaries from the MEaSUREs Antarctic Boundaries for International Polar Years 2007–2009 from Satellite Radar, Version 249. To account for any ice-shelf growth, a 10 km buffer zone was allowed in front the shelf front. When processing, any elevation with an accompanying ≥30 dB of backscatter was removed, following ref. 74, before outlier elevations (e) were removed if they fell outside the range

$${e}_{25}-1.5\,\times \,{\mathrm{IQR}},\,{e}_{75}\,+1.5\,\times \,{\mathrm{IQR}}.$$ (7)

A tidal-height correction using the Circum-Antarctic Tidal Simulation75, the ocean tide loading from the Finite Element Solutions 2004 tidal atlas76 and geoid from Earth Gravitational Models 1996 were used to correct altimetry data for variations in freeboard elevation. Ice-shelf thicknesses were calculated following refs. 77,78, which used

$$h=\frac{(e-\delta ){\rho }_{{\rm{w}}}}{{\rho }_{{\rm{w}}}-{\rho }_{{\rm{i}}}}+\delta ,$$ (8)

where h is ice-shelf thickness, e as the freeboard elevation, δ is the firn air content79, ρ w = 1,028 kg m−3 is the water density and ρ i = 917 kg m−3 is the ice density.

Individual grids were built for each ice shelf with a resolution of 3 km2 at a yearly temporal resolution. This set-up was chosen to prevent erroneous readings due to nearby ice shelves from influencing the ice-shelf thickness output, and to reduce processing time and counteract a lower density of altimetry tracks near the edge of the ice shelf77. The distance of each grid cell from the shelf front was calculated, with only grid cells within 20 km of the shelf front retained to represent the approximate area lost in each calving. Each grid cell from 2014–2023 was combined, using the IQR method to remove outliers, such as those that picked up the surrounding ice sheet or fast ice. Quartiles were calculated for each ice shelf analysed in Fig. 5b that did not have a defined pre-calving ice-shelf thickness.