Grounding line discharge and ice velocity
We estimate grounding line discharge, D, across each flux gate pixel as
$$D={VHw}\rho,$$ (1)
Where V is the gate-normal ice velocity, H is the ice-equivalent thickness, w is the pixel width and ρ is ice density (917 kg m−3). We use “flux gate 1” (or FG1) from ref. 65 for our flux gate, which follows areas of direct bed observations and low bed elevation error, rather than the grounding line, and we resample the flux gate pixels to regular 200 m spacing. All plots in the main text are produced using this flux gate. We tested the sensitivity of our results to flux gate location by migrating each gate pixel along a flowline upstream based on the time-average velocity in 0.5-year increments for 10 years, creating 20 individual flux gates. Ice discharge increases by up to 6.3%, and by 1% on average, compared to that using the FG1 gate, as the gate is migrated upstream (Supplementary Fig. 12), reflecting the difficulty of conserving mass with imperfect thickness and velocity data, rather than algorithmic errors.
The gate-normal ice velocity is given by
$$\sin \left(\theta \right){V}_{x}-\cos \left(\theta \right){V}_{y},$$ (2)
Where \({V}_{x}\) and \({V}_{y}\) are the easting and northing ice velocities, as defined by the South Polar Stereographic grid (EPSG3031), respectively, and \(\theta\) is the angle of the flux gate relative to the same grid. Ice velocity data are compiled from multiple published and freely available sources during 1996–2018 and intensity tracking of Sentinel-1a and 1b images since 2015. For 1996 velocities, we use 450 × 450 m MEaSUREs InSAR-based estimates derived from 1-day repeat ERS-1 imagery4,64, which covers the region spanning Cosgrove to Kohler Glacier, and 200×200 m velocities from ERS offset tracking over the Getz basin, provided by ENVEO (https://cryoportal.enveo.at/data/), which has been filled using an optimised ice-sheet model – BISICLES5. For 2000 and 2005 to 2016, we use 1×1 km MEaSUREs annual velocity mosaics62,63. We also use a MEaSUREs velocity mosaic incorporating velocity estimates between 1995 and 2001, at 450×450 m resolution91. For 1997–2018, we use 120×120 m ITSLIVE annual mosaics65; for 2009 to 2020, we use quarterly ice velocity mosaics from TerraSAR-X (2009–2016) and Sentinel-1 (2015–2020) just over Pine Island Glacier7. From 2015 onwards, we use high temporal resolution (6- to 12-day) estimates of ice velocity derived from intensity tracking of all 6- and 12-day Sentinel-1 image pairs in GAMMA, using standard methods5. We combine a range of cross-correlation window (ranging from ~1.2 × 1.2 km to ~500 × 500 m, each with 75% overlap), which are merged using a signal-to-noise ratio weighted mean, to accommodate the large range in ice velocities in the ASE and to increase coverage. Individual image pair velocities are filtered using a ‘dusting’ approach to remove spatially-isolated pixels and a 2D gaussian filter5, and posted at 100 × 100 m. Image pair velocities are mosaicked to 200 × 200 m monthly-mean velocities after further outlier removal and smoothing using a Discrete Cosine Transform Penalised Least Squares approach, optimised using the generalised cross-validation for fast, unsupervised processing92,93. Each of these velocity products spans a time period; following ref. 94, we ignore this and treat each product an instantaneous measurement with the timestamp given by the central date in the estimate. We extract easting and northing velocities at each gate pixel using nearest neighbour interpolation. Treating each gate pixel as a timeseries, we remove outliers, which we define as data points more than three standard deviations from the detrended timeseries. We fill temporal gaps using a linear interpolation except at the beginning and end of each timeseries, which are back- and forward-filled with the temporally nearest value for that pixel. Each filled timeseries (with no data gaps) is smoothed using a Savitzky-Golay filter with a 6-month window and a second-degree polynomial. Finally, as in previous studies4,94, we assume the depth-averaged velocity is the same as the surface velocity.
We estimate ice thickness, and therefore discharge, at the timestamp of each velocity measurement using a fixed bed elevation and a time-varying ice surface. The assumption of a fixed bed means that we neglect any changes in bed elevation due to erosion or uplift, which we expect to be at least an order of magnitude smaller than the observed ice surface elevation changes. We use BedMachine version 2 bed topography95,96 and the 200 m REMA Digital Elevation Model97 (DEM) to define our baseline thickness estimate. We use observed ice-equivalent thinning rates from 2003–2019 derived from ICESat-1 and ICESat-2 observations66 to adjust the REMA DEM surface, assuming that the REMA DEM is timestamped to 9th May 201597. For velocity measurements prior to 2003, we use the average estimated elevation from 2003 to 2006, rather than extrapolate the thinning rate. Similarly, for velocity measurements after 2019, we use 2017–2019 average ice surface elevation. Our results are insensitive to these thickness adjustments (Supplementary Fig. 13). Ice-equivalent thickness is estimated at each velocity epoch by removing a time-varying firn air content, provided by the IMAU Firn Densification Model98,99, forced by RACMO2.3p254 (Supplementary Fig. 14). Our flux gates are upstream of the grounding line; we remove the effect of dynamic thinning or thickening and SMB changes that occur between the flux gate and grounding line by integrating the observed ice-equivalent elevation change rates and modelled SMB between the flux gate and grounding line100,101. Relative to the final discharge estimate, these corrections are modest but non-negligible: 2.6% for dynamic thinning and 2.5% for SMB, across the full ASE flux gate. For all discharge estimates, we use the MEaSUREs grounding line100,101, and so we ignore changes in grounding line position through time.
Grounding line discharge from each basin is estimated by integrating the pixel-based discharge estimate for all flux gate pixels that fall within the respective basin. In the main text, we use basins from ref. 100. In this paper, we define the ASE as MEaSUREs basins G-H and F-G (i.e. we include the Getz basin in our ASE totals). Basin-scale discharge using alternative basin definitions102 are shown in Supplementary Table 1.
Discharge error
Following ref. 94, we define our discharge error using the minimum and maximum possible discharge based on the errors in the velocity and thickness data at each gate pixel. The pixel-based maximum discharge is defined as
$${D}_{\max }=\left(V+\sigma V\right)\left(H+\sigma H\right)w\rho,$$ (3)
And the pixel-based minimum discharge is defined as
$${D}_{\min }=(V-\sigma V)(H-\sigma H)w\rho .$$ (4)
\(\sigma V\) and \(\sigma H\) are the timestamped, pixel-based errors in the velocity and thickness estimates, respectively. Where available, we use the easting and northing velocity errors provided in each product. Where we have interpolated the velocity, we define the error as 10% of the estimated easting and northing velocity components in each timestamped pixel. We combine the easting and northing velocity errors normal to the flux gate through quadrature. Similarly, the thickness errors are defined as the sum through quadrature of the BedMachine bed elevation error where available and the error in the timestamped surface elevation estimate. We assume a 1 m error in the baseline REMA 200 m DEM, which is equivalent to the 90th percentile of the errors in the mosaic97, and a 0.1 m yr−1 error in the linear surface elevation change66. Following ref. 94, we define invalid thickness estimates as those that are <20 m where the ice is flowing faster than 100 m yr−1, and we estimate the thickness in those pixels using a weak correlation (R2 = 0.55) fit between the time-average velocity and the remaining valid baseline thickness pixels and we assume a 50% thickness error in these locations. Across the ‘FG1’ flux gate, this correction replaces 9.5% of all gate pixels and increases our average FG1 gate thickness by 26 m, resulting in an 1.2% increase in discharge.
In all cases throughout the manuscript, the presented discharge errors are defined as
$$\sigma D=\varSigma \left(\left({D}_{\max }-D\right)+\left(D-{D}_{{{\min }}}\right)/2\right),$$ (5)
where the pixel-based discharge errors are summed across all gate pixels falling within the defined basin. Therefore, we do not perform any mathematical reduction of the errors presented in the manuscript. Across the ASE as a whole, we obtain an average error of 14.9%. If instead we summed our pixel-based errors through quadrature, the average error would be 0.4%.
Surface mass balance
We use SMB output from three regional climate models. We draw predominantly on the Regional Atmospheric Climate Model version 2.3p2 (RACMO2.3p2), which is forced at its lateral boundaries by ERA5 and run at 27 × 27 km resolution54. It has been especially adapted for use in the polar regions by adopting a multilayer snow model that includes routines for melt, percolation, refreezing, runoff, snow albedo, and drifting snow. RACMO2.3p2 has been extensively evaluated and compares favourably with observations of near-surface climate and SMB; in particular, RACMO2.3p2 SMB estimates are similar to and strongly correlated with radar-derived snow accumulation in the Getz and Thwaites Glacier basins54.
These data are supplemented by SMB estimates from the Modèle Atmosphéric Régional (MAR)56,57 and the Danish regional climate model HIRHAM555. MAR is a hydrostatic model specifically designed for polar areas. It is forced at its boundaries by ERA5, runs at 35 × 35 km resolution, and includes routines for meltwater refreezing, snow metamorphism, snow surface albedo, and drifting snow57. MAR monthly-mean SMB output is available for 1980–2021. HIRHAM5 is a hydrostatic model forced at its boundaries by ERA-Interim. It is run at 0.11° resolution and includes a new albedo scheme55, sublimation, surface melt, retention, refreezing, and runoff103 but does not include wind redistribution of snow. HIRHAM5 daily-mean SMB output is available for 1980–2018.
Cumulative SMB anomalies are calculated as follows. In each spatial cell, we cumulatively sum timeseries of anomalies relative to the 1980–2008 climatological mean of that cell. We use 1980–2008 rather than the 1979–2008 mean so that we can use the same climatological period for all SMB datasets. These pixel-based cumulative SMB anomalies are integrated across each glacier drainage basin to obtain regional cumulative SMB anomalies. SMB anomaly errors are taken as the root sum square of the annual errors5,70, which we take as 14.8% of the annual SMB for all products2.
Mass balance
We linearly interpolate both the discharge and SMB timeseries to monthly intervals from January 1996 to December 2021. The mass balance is then the simple difference between these variables. The mass balance error is assumed to be the root sum square of the discharge and SMB errors. Cumulative mass change errors are defined as the root sum square of the annual mass balance errors, excluding errors in initial ice thickness.
Impact of snowfall drought on historical mass balance
Annual grounding line discharge from the ASE is estimated by linearly interpolating the values presented in ref. 4 to annual intervals. Ref. 4 did not present grounding line discharge estimates for Getz. For the purpose of this analysis we estimate grounding line discharge from the Getz basin during 1974–1996 by scaling the annually interpolated values from ref. 4 according to the average relative contribution of the Getz basin to the total ASE discharge presented here since 1996. This approach assumes a synchronous response from all basins which is contrary to observations9. These discharge estimates are therefore used only to illustrate the potential impact of historical SMB changes on ASE-wide mass balance and are not a true estimate of historical ASE grounding line discharge.
Atmospheric rivers
We use an existing ice-sheet wide, 3-hourly catalogue of ARs, derived using a Polar-specific AR detection algorithm45, which is based on contiguous areas of very high (98th percentile) meridional component of vertically-integrated water vapour transport (vIVT). It is calculated using vertically-integrated specific humidity and meridional wind speeds from MERRA-2. From this catalogue, we include ARs only if they intersect the grounded ASE (basins G-H and F-G) and amalgamate AR detections if they are separated by less than 12 hours (accounting for any changes in AR footprint during that time). The precipitation provided by ARs is calculated by integrating daily RACMO2.3p2 precipitation within the landfalling-portion of the AR footprint from the start of the AR detection period and until 24 hours following the end of the event.
Climate indices and historical snowfall drought
We use four climate indices to provide insight into the mechanistic drivers of potential historical periods of snowfall drought. The first of these indices is the Oceanic Niño Index (ONI) is a key indicator of the El Niño Southern Oscillation. It is a rolling 3-month average sea surface temperature anomaly (relative to the mean of the surrounding 30 years) in the east-central tropical Pacific between 120°–170°W (the Niño-3.4 region). An El Niño is considered to be present when the ONI is +0.5 or greater, whereas a La Niña is considered active when the ONI is −0.5 or lower. In Fig. 6, we shade sections of the time series only when the portion of the time-series between each zero crossing contains at least one value that exceeds these thresholds. The second climate index is the Tripole Index for the Interdecadal Pacific Oscillation104 (IPOTPI), of which we use the HADISST1.1 filtered version. The IPOTPI is based on the difference between the sea surface temperature anomaly averaged over the central equatorial Pacific (10°S–10°N, 170°E–90°W) and the average of the sea surface temperature anomaly in the Northwest (25°N–45°N, 140°E–145°W) and Southwest Pacific (50°S–15°S, 150°E–160°W), with anomalies calculated relative to January 1971 to January 2000 period. We also use the Southern Annular Mode (SAM)105, which is a station-based index defined as the zonal pressure difference between 40°S and 65°S. Finally, we also use an ASL index, defined as the annual average of the 500-hPa geopotential height from 75°S to 60°S and 170°E to 70°W, normalised relative to the 1941–1990 period, derived from proxy-based climate reconstructions80.
We define periods of snowfall drought (red and pink shading in Fig. 6) as periods of annual snow accumulation from ice cores or annual SMB from RACMO2.3p2 that are >1.5 median absolute deviations below the 1979–2009 median for at least two years.
Ice-equivalent thickness change due to snowfall drought
The total mass change due to the 2009–2013 snowfall drought was ~260 Gt. We estimate these-equivalent thickness change due to this mass change as follows. When distributed evenly over 1 m2, 1 mm of water of density 1000 kg m−3 weights 1 kg. Therefore, the ice surface height change in water equivalent, H w.e due to the 260 Gt mass change, M, is
$${H}_{w.e}=M/A,$$ (6)