Materials and Methods

GRACE data
We make use of the COST-G RL01 Level-3 GRACE and GRACE-FO V0002 time series obtained from http://gravis.gfz-potsdam.de/antarctica(44, 45 ) at 50 km gridded resolution (noting that GRACE estimates of mass change are correlated over 200-300 km). At the time of download the series spanned Mar 2002 to Mar 2021 inclusive. As described by ref (45 ), post-processing steps include replacement of coefficients C20, C30 (only for the months starting from November 2016), C21 and S21 and its formal standard deviations by values estimated from a combination of GRACE/GRACE-FO and Satellite Laser Ranging (SLR); insertion of geocentre coefficients (C10, C11, S11); and estimation and removal of a periodic 161 day-period signal due to mismodeled and aliased signal of the S2 tide. These series come corrected for glacial isostatic adjustment using the ICE6G_D(46 ) model but we restored this model and removed that of Caron et al.(47 ) with values for both given in Table S1. The choice of GIA model will affect the linear trend in estimated mass change but not variations. Qualitative comparison of time series from an alternative GRACE solution (GFZ RL06) showed modest differences at inter-annual timescales, and generally within uncertainties, but with COST-G time series temporally smoother.
We also considered the supplied time series for 25 individual basins according to the definitions of ref (48 ), taking EAIS to be the combination of basins 302 to 317, WAIS 318 to 323 and 301, and APIS 324 to 325, all inclusive. These basin series benefit from a GRACE data inversion which includes forward modelling of sub-basin-scale mass distribution(44 ). Further details on the GRACE time series are provided at ftp://isdcftp.gfz-potsdam.de/grace/GravIS/COST-G/Level-3/ICE/GravIS_ICE_Technical_Note.pdf
SMB model outputs
We use modelled SMB from the RACMO2.3p2 27km model(30 ) covering Jan 2002 to Feb 2021 in units of kg/m2/month. We computed cell-by-cell SMB anomalies relative to the cell-mean SMB over the full data period. These were then cumulatively summed and converted to units of giga-tonnes (Gt). For differencing with GRACE series, we interpolated the SMB grids to the GRACE grid spacing and then interpolated them to the GRACE time steps.
Climate indices
For SAM, we adopt the station-based index of Marshall et al.(18 ) obtained from http://www.nerc-bas.ac.uk/icd/gjma/sam.html . For ENSO, we use the Niño3.4 index obtained from https://psl.noaa.gov/gcos_wgsp/Timeseries/ based on the HadlSST record over 5°N-5°S, 170°W-120°W. For each index, we truncated the portion prior to the GRACE period, then normalized, cumulatively summed, and renormalized the series to produce the SAMΣ and Niño3.4Σ series shown in Fig. 1 and Fig. S1. In the presence of long-period index changes, the cumulative sum might depend on the chosen reference period which is used to calculate climatological mean and derive anomalies from it. For both SAM and Niño3.4 indices we adopted a uniform reference period of 1971-1999 inclusive which provides the advantage of being a well observed period outside the GRACE data window. The Niño3.4Σ series is largely insensitive to the choice of reference period as the historical Niño3.4 index has little long-term tendency. In contrast, using 1971-1999 as the reference period very likely underestimates the trend in SAMΣ, as the SAM index has exhibited pronounced upward trend since the mid-1970s, largely due to anthropogenic climate change(34 ). Adopting a reference period later than 1971-1999 would result in removal of the shift to positive SAM, which would be unrealistic.
The well document shift toward the positive SAM phase since the 1940s requires that SAMΣ has a long-term positive trend. We discuss below that working with SAM or SAMΣ is mathematically identical in the absence of GRACE data noise but working with SAMΣ has distinct advantages when working with real GRACE data.
To explore the sensitivity to the choice of SAM or ENSO indices, we tested with alternative indices in the multivariate regression that is described below. Instead of the station-based Marshall SAM index we tested an index computed using zonal mean sea level pressure difference between 40°S and 65°S from the JRA55 reanalysis model (referred to as “SAM JRA55” in Fig. S1). The zonal means at 40°S and 65°S were each calculated. Each monthly value was then normalized over 1981-2010 and then the 65°S values were subtracted from the 40°S values and the series renormalized over 1971-2000. Instead of the Niño3.4Σ we tested the Southern Oscillation Index (referred to as SOI in Fig. S1) obtained from https://www.ncdc.noaa.gov/teleconnections/enso/soi, computed as the normalized difference in standardized sea level pressure between Tahiti and Darwin, renormalized over 1971-2000. For consistency with Niño3.4 we use -SOI. The indices are shown in Fig. S1 along with their cumulative sum and detrended cumulative sum. The SAM JRA55 index differs to the station-based SAM index of Marshall et al. (18 ) over inter-annual timescales, but its cumulative sum is in close agreement. -SOI has a small negative mean anomaly over the GRACE period resulting in a small negative trend in -SOIΣ along with some inter-annual differences from Niño3.4Σ. There is no evidence of a long-term trend in ENSO and this is likely a result of a small positive bias in SOI during the reference period. Other reanalysis products are either not available up to the end of the GRACE period or are only available after 1979 and hence do not span our full reference period commencing in 1971.
EOF analysis
We perform a standard EOF analysis using GRACE data on an Antarctic Polar Stereographic grid with 50 km resolution (Mar 2002 to Mar 2021). Before this we removed a linear trend at each grid point (computed using ordinary least squares and considering uneven data sampling) to focus on variability. Periodic terms were not removed before the EOF analysis. Unlike modes 1 and 2 shown in Fig. 1, the remaining modes are dominated by higher-frequency signals and often exhibit less-organized or noise-like spatial patterns, with their explained variances each less than 6%. To reveal processes underlying the two leading modes, we regressed cumulatively summed and linearly detrended SMB from RACMO2.3p2_ANT27, as well as sea level pressure, and each of the two components of 10m wind from ERA5 reanalysis, onto the PCs of the first two modes.
Given the temporal correlations evident in the dominant modes (Fig. 1), we also tested Extended EOF analysis(49 ) (EEOF), which considers temporal correlations in the data. Using a lag of 12 months we computed EEOFs after gap filling and interpolating to a constant monthly timestep and found the resulting two leading EOFs to be almost temporally constant and with PCs that are lag-smoothed versions of the PCs from the conventional EOF analysis. Given this agreement, and for simplicity, we retained the standard EOF analysis in the main text.
Multi-variate analysis
Using ordinary least squares, we solved the coefficients (a, b, c, d, and e) of the functional model describing time-evolving mass (M )
\(M\left(t_{i}\right)=a+b\left(t_{i}-t_{0}\right)+\sum_{k=1}^{3}{\left(c_{k}^{s}\sin{\left(2\pi f_{k}t_{i}\right)+}c_{k}^{c}\cos\left(2\pi f_{k}t_{i}\right)\right)+d\text{SAM}_{\Sigma}+e{Nino3.4}_{\Sigma}}\)(1)
Where fk = [1, 2, 365.25/161] cycles per year, with the third frequency describing the S2 tidal aliasing period. Separate periodic coefficients were estimated for each of GRACE and GRACE-FO. \(\ \)While the S2 term was apparently removed in the GRACE pre-processing(45 ), we found evidence of it in the GRACE-FO period and re-estimated them to provide realistic uncertainties. We adopted \(t_{0}\) as the mid point of the GRACE series.
We explored the impact of the presence of a trend in SAMΣ upon its estimated regression coefficient. To do this we repeated the regression after removing a linear trend from SAMΣ. We found that the SAMΣ regression coefficients were unchanged to one decimal place (compare Tables S1 and S3). This indicates that the variability of the cumulative indices (Fig. S1c) is dominating the regression and that the SAMΣcoefficient (d ) is largely independent from the linear regression coefficient (b ).
The inclusion of SAMΣ in the regression model introduces two types of uncertainties. First, it is very likely that using 1971-2000 as reference period underestimates the trend in SAMΣ. And, so, the signal attributed to SAM might also be underestimated as a result. For example, adopting a reference period of 1961-1990 attributes a larger portion (64% of the total change) of total AIS mass loss to the SAMΣ term in the regression. However, using an earlier reference period would introduce another layer of uncertainty in the quality of earlier data (for example, missing station values, less observations to assimilate in reanalysis products) used to define SAM. These sensitivities would remain if the regression was working with unsummed SAM and time-differenced GRACE. Second, SAMΣ has a quasi-quadratic shape over 2002-2016 (Fig. S1c). This has broken down since 2016 and this allows separation of this term from a pure acceleration such as could be due to monotonically increasing ice discharge. Adding a quadratic term to Equation 1 gives a mathematical correlation of 0.53 with the ENSO term and 0.63 with the SAM term, indicating moderate correlation. As such, there might be some aliasing of any quadratic-like mass change that are not related to SAM and ENSO into the estimated SAM and ENSO terms. We note, however, that recent ice stream speedup has been in rapid response to variable climate forcing rather than being slow and monotonic, especially in West Antarctica where most of the modest velocity change has occurred over the GRACE period(50-52 ). In addition, we note the presence of periodic decadal climate variability in Antarctica(53, 54 ) (see also the SMB decadal variability in Fig. 1b, e, Fig S10-11). Our analysis in the main and supplementary text suggests that the dominant GRACE modes are linked to SAM and ENSO and much of the estimated SAM and ENSO signal is associated with SMB.
To compare with the results of linear regression without the climate indices, we also performed a univariate regression removing the last two terms of Equation 1:
\(M\left(t_{i}\right)=a+b\left(t_{i}-t_{0}\right)+\sum_{k=1}^{3}\left(c_{k}^{s}\sin{\left(2\pi f_{k}t_{i}\right)+}c_{k}^{c}\cos\left(2\pi f_{k}t_{i}\right)\right)\)(2)
The results of evaluating Equation 2 include those shown in Fig. 2d and Fig. S8 (orange).
To examine the potential contribution of SMB to the GRACE-derived regression coefficients, we repeated the above regressions after subtracting from GRACE time series the cumulative, detrended SMB output from the RACMO2.3p2_ANT27 model. We also performed the regression just on the cumulative, detrended SMB from the RACMO2.3p2_ANT27 model outputs.