American Aerial photographs
We utilize 468 black and white aerial photographs of the Wordie Ice Shelf collected mainly by the US NAVY from 1966 to 1969 (Table 1). The images were captured using a trimetrogon camera setup, wherein left-oblique, on-nadir (vertical), and right-oblique images are captured simultaneously21.
Table 1. Overview of all utilized aerial images
Flightline
|
Date
|
Number of Images
|
1834
|
November 28, 1966
|
150
|
1835
|
November 28, 1966
|
156
|
1843
|
December 27, 1966
|
6
|
1844
|
December 27, 1966
|
51
|
1845
|
December 27, 1966
|
36
|
2167
|
January 12, 1969
|
69
|
To orthorectify the aerial photographs and create combined orthomosaics for each year, we use a standard structure from motion workflow(43) in Agisoft Metashape v. 2.0.2. Photographs from different years undergo separate processing. First, the images are scaled and aligned by detecting the image fiducials using the Detect Fiducial function. All fiducials are manually inspected and misplaced ones are corrected. Subsequently, the images are aligned by extracting up to 50,000 key points per image(44). In certain instances, additional tie points are manually placed to ensure precise image alignment. Next, we generate a dense point cloud at maximum image resolution. The dense point cloud is then gridded into a DEM and used to create a mesh, from which the orthomosaics are generated. As the produced orthomosaics and DEMs are generated in relative coordinates they are georeferenced to satellite images of the region. The orthomosaics and DEMs are generated at similar pixel resolution, which allows us to transfer the ground control points (GCPs) from the ortho-georectification to the DEM. We georeference the orthomosaic using distinct surface features such as nunataks and islands and distribute the control points evenly across the image. To constrain the georeferencing over the shelf area, we use ice rises detectable in the 1974 and 1979 Landsat images. However, we acknowledge that these points may exhibit larger uncertainties compared to the control points over stable bedrock. The georeferencing is performed using an affine transformation, with the number of control points placed varying from 15 to 25 depending on the coverage of the orthomosaic. The georeferenced orthomosaics have a resolution of 4-6 m and a Root Mean Square (RMS) error from the placed control points ranging from 40.1 to 43.2 m. Among the three generated orthomosaics, the one from November 28, 1966, is based on the largest selection of images (Table 1) and provides the most extensive and complete coverage of the entire shelf area, along with the glaciers flowing into it (Fig. S1). For the georeferencing of this orthomosaic a total of 25 control points were placed, achieving an RMS error of 41.8 m.
Following the georeferencing of the orthomosaic, we transfer the GCPs to the corresponding DEM to ensure consistent georeferencing between the two products. Attempts were made to produce historical DEMs with high absolute elevation accuracy for measuring ice surface changes. However, due to the old age of the images, they contain varying degrees of damage, which yielded high vertical inaccuracies. Despite this, the relative elevations of the precisely georeferenced DEMs are still useful for accurately detecting changes in slope.
Grounding line
Historical grounding line positions from 1966 are determined using the break of slope method from Mouginot et al. (2018)44. he grounding line mapping is performed using a slope and a hillshade raster (Fig 3), generated from the November 28, 1966 georeferenced historical DEM at 100x100 m resolution to even out local errors. In addition to the historical grounding lines, we incorporate MEaSUREs 1996 data derived from Differential Synthetic Aperture Radar (DInSAR) data33. However, as no DInSAR grounding line is available for Carlson Glacier, we utilize the data from Antarctic Surface Accumulation and Ice Discharge (ASAID)34 project derived from the break of slope using Landsat 7 imagery and ICESat altimetry data with a temporal coverage from 1999 to 200345.
In an effort to retrieve recent grounding line locations, we used DInSAR measurements from the EU/ESA Sentinel-1 satellites and the conventional method of double-difference interferometry46. We processed 200 image pairs from 2016-2023 with a temporal baseline of 6 days (when available, otherwise 12 days), however, in all cases, interferometric coherence was consistently lost over the full region of interest, meaning that no useful vertical displacement measurements could be extracted. Applying the method of differential range offset tracking47 to the Sentinel-1 measurements equally did not yield meaningful measurements. This limitation may be due to the difficulty of tracking the GL with DInSAR for a relatively small floating section (e.g limited resolution and loss of coherence due to precipitation, surface melting, or fast ice flow and high deformation).
Modern grounding line estimates from 2014 are available for Fleming Glacier30,31, with notable disparities in the findings. Friedl et al. (2018)31 suggest the grounding line to be situated 5-8 km further upstream than Walker and Gardner (2017)30. We have opted to adopt the Friedl et al. (2018)31 estimate, as they calculate the hydrostatic equilibrium of the glaciers center, whereas Walker and Gardner (2017)30 refer to the 1996 grounding line. Moreover, utilizing optical satellite imagery48 we determine approximate 2023 grounding lines for FIP1, FIP2, and Prospect Glacier. This methodology is employed only when the surface shading in the optical satellite imagery provides a clear distinction between the grounded and floating sections of the glacier (Fig S10). As a result, we are unable to determine optical-derived grounding lines for Hariot, Unnamed, Carlson, and Fleming Glacier. By comparing the 1996 DInSAR grounding line to the break of slope from the surface shading in the 1990 optical images (Fig. S10), we observe a close alignment with the largest variation of 1000 m observed at the FIP2. The average variation is ~300m between the grounding line estimates. Thus, we adopt a conservative uncertainty of 500m for our break of slope grounding line estimates.
The grounding zone marks the region where the glacier transitions from being completely grounded to floating. As a result the different mapping techniques map features within the grounding zone that acts as proxies for the true grounding line48. Consequently, some of the changes observed in the grounding line may be attributed to variations in the mapping method, and to the tide magnitude at the time of measurement. However, due to the close alignment between 1996 DInSAR and the 1990 break of slope grounding line (Fig. S10), we anticipate this factor to have minimal influence.
Velocity
We compile a time series of surface ice velocities spanning over 50 years by integrating data from multiple satellites and aerial orthomosaics. Velocities from 1966 to 1968 are obtained from orthomosaics, which are precisely georeferenced prior to feature tracking using stable bedrock areas. For three glaciers (Unnamed, Fleming, and FIP1), we manually track the movement of distinct surface features on the glacier surface between 28-11-1966 and 12-01-1969. The uncertainties of the manual tracking are determined as the combined uncertainty of the georeferencing error, and the manual error attributed with feature positioning which is set to 1 pixel. This approach gives a combined uncertainty of 60 m/yr. For the Carlson glacier, the relatively short time between repeat acquisitions, approximately 1 month (28-11-1966 and 27-12-1966), enables us to employ automatic feature tracking. We first perform a sobel filter in the x and y direction to enhance surface features on images. Finally, we feed the resulting images to the feature tracking algorithm using image correlation windows of 64x64 pixels, search windows of 20x20 pixels with spacing of 10 pixels. From 1986 to 2021, we use time series of surface flow velocities originating from feature and speckle tracking using satellite observations from ESA's ERS-1&2, Sentinel-1a/b and Sentinel-2, RADARSAT-1&2, JAXA's ALOS/PALSAR, NASA's LANDSAT-1-8 (methodology similar to Mouginot et al., 201249, Mouginot et al., 2017 50 Millan et al., 201951). To account for missing years of interest, we augment the time series with data from the ITS_LIVE annual mosaics52. Velocities are extracted along an approximate center flowline of each glacier. The flowlines are positioned as close to the center as possible, with slight adjustments made based on areas with the most complete historical velocity estimates. For each year, we exclude all measurements seaward of the frontal position to avoid comparing velocities over sea ice or melange.
Ice shelf area and pinning points
Frontal positions of the Wordie Ice Shelf and its glaciers spanning 1966 to 2022 are mapped using the generated orthomosaics, optical Landsat images (1974-1991 and 2001-2022), with a preference for images captured during the austral summer. The ice shelf area is determined with respect to the MEaSUREs 1996 grounding line. We also supplemented our analysis with data from ESA's ERS-1 and ERS-2 satellites in the period 1991-2004. Single Look Complex (SLC) images were processed into amplitude and geocoded using the 200 meter resolution DEM from the Radarsat Antarctic Mapping Project (RAMP). Given the lower resolution of ERS backscatter images and distortions due to the different satellite geometries, the frontal positions from 1992-2000 are more difficult to accurately delineate than for the Landsat data.
We map all ice shelf pinning points found in both the 1966 orthomosaic and 1974 Landsat images to exclude ephemeral features, as these may gain contact or unpin due to tidal processes and thus provide limited buttressing 53. Subsequently, we map when the ice shelf unpins from the pinning point either due to thinning or due to retreat (Fig. 1 and Fig. S2). However, during the period from 1993 to 1999, when we rely on ERS radar images, we are not able to determine whether pinning points disappear due to thinning.
Iceberg production
Throughout the entire study period, we generated annual estimates of tabular iceberg production by manually digitizing icebergs with a surface area approximately greater than 0.5 km². We use a combination of both aerial, satellite, optical and radar imagery and map icebergs covering the entire Wordie Bay, spanning from the ice shelf front extending to Mushroom, Puffball and Ramirez Islands. This allows us to study changes in iceberg size, number and timing of release with respect to ice shelf retreat and area changes. This also provides important qualitative information on the status and the changes in ice mélange through time, since it can have an impact on the glacier retreat54. Iceberg sizes can indeed be considered as an indirect indicator on the local thickness of the ice mélange, hence its stabilizing power. Understanding the production of tabular icebergs is also important to serve as a proxy for the dynamics and calving regime of ice shelves, particularly in the context of collapse. Based on the presence of tabular icebergs from year to year, we estimate the residence time inside the bay during the entire study period to be spanning between 1-3 years, depending on the size (and thickness) of icebergs. Thus, we analyze the long-term trend in iceberg production throughout the time-series rather than the year-to-year variations, which might be biased by double counting of elements.
Climate, SMB, ocean, and sea ice concentration of Wordie region
We extract RACMO 2.3p2 mean annual surface mass balance (SMB), snowfall, rainfall, and austral summer air temperatures (December to February) from 1979 to 202255 for the Wordie region (Fig. S5). Additionally, we calculate the mean austral summer temperature 1941-1982 using 0.25° ERA5 data56 and since 1976 using observational records from Rothera Station located approx. 150 from Wordie Ice Shelf extracted from the Scientific Committee on Antarctic Research (SCAR) Met-READER project (https://legacy.bas.ac.uk/met/READER/). Monthly sea ice concentrations in the Wordie Bay area are extracted from the Bootstrap Sea Ice Concentrations dataset57 derived from Nimbus-7 and Defense Meteorological Satellite Program (DMSP) satellites. Subsequently, we calculate the annual mean November to March sea ice concentration in the region.
To track the progression of ocean temperatures in the Marguerite Bay, we used Conductivity Temperature and Depth (CTD) measurements obtained from the Hadley Centre (bodc.ac.uk) spanning the years 1960 to 2019 (Fig. 3). Using the temperature profiles, we calculated the average temperature between 150 and 500 meters. We then used a locally weighted scatterplot smoothing58 to plot the temporal evolution of the warm water layer through time.