We estimate mass changes of the Antarctic Ice Sheet from 2002 to 2020 based on data from the GRACE (2002–2017) and the GRACE-FO (launched 2018, operational) satellite missions that provide nearly continuous monthly measurements of the time-varying gravity field of the Earth. Based on the Stokes potential coefficients (Level-2 data) issued by the missions’ Science Data System (SDS), which are a spherical harmonic representation of the Earth’s gravity field changes, mass redistribution on the Earth surface can be derived (Tapley et al., 2004, 2019). We use Release 06 (RL06) Level-2 data from the three GRACE/GRACE-FO SDS teams; German Research Centre for Geosciences (GFZ) (Dahle et al., 2019), Centre for Space Research University of Texas, Austin, United States (CSR) (Bettadpur, 2018; Save, 2019) and Jet Propulsion Laboratory, California Institute of Technology (JPL) (Yuan, 2019), as well their combination (AV RL06) (Sasgen et al., 2020). A comparison of time series of mass change obtained from the four Level 2 GRACE/GRACE-FO data sets is shown in Supplementary Figure S1.

The data capture 191 out of 218 possible months during the missions’ operation; 163 monthly solutions from GRACE covering April 2002 to June 2017 and 28 solutions from GRACE-FO covering June 2018 to November 2020. To estimate uncertainties in the derived mass change time series associated with processing choices and corrections, we additionally consider Level-3 mass change gridded products for Antarctica provided by the Technical University of Dresden (TUD) published via the GravIS data portal (Sasgen et al., 2019a) and the CSR RL06 Mascons (Save et al., 2016). All GRACE/GRACE-FO Level 2 and 3 data used in this study are publicly available (Data availability).

We invert the monthly GRACE/GRACE-FO gravity fields for mass changes in Antarctica using a spectral inversion method similar to Jacob et al. (2012), modifying the forward modelling approach in the spatial domain previously applied to Antarctica (Sasgen et al., 2013). For this, we define spatial patterns of mass change, from either satellite observations or modelling, within 25 drainage basins of the ice sheet (Figure 1), calculate their spectral gravity field signature, and adjust these spectra by scaling in order to minimize the difference to the observed spectrum represented by the GRACE/GRACE-FO Stokes coefficients. Prior to adjustment, we limit the representation of the GRACE/GRACE-FO coefficients to the region of Antarctica by employing a spectral mask of the ice sheet using a buffer zone of 200 km (Martinec, 1989). Finally, we obtain the total mass change within in each drainage basin by spatial integration of the scaled mass distribution underlying the forward models. The optimal adjustment factors are estimated by weighted least-squares, allowing us to account for the noise structure in the data set of GRACE/GRACE-FO coefficients.

We compose a forward model by assuming that the spatial pattern of mass change corresponds to CryoSat-2 rates of mass change (Sasgen et al., 2019b) for trends and accelerations, and the mean surface mass balance (van Wessem et al., 2018) for the remaining, de-trended temporal components. We find that evaluating GRACE/GRACE-FO data at the basin scale is relatively insensitive to the mass distribution underlying the forward model. We quantify the uncertainty of our results associated with these processing choices by comparing with the independently derived Level-3 mass change products TUD COST-G (Sasgen et al., 2019a) and the CSR RL06 Mascons (Save et al., 2016) (Data availability). A sensitivity analysis using a spectral circular mask extending from the South Pole to 60°S, as well as using no masking produces consistent results.

An advantage of the spectral inversion method is that uncertainties of the GRACE/GRACE-FO coefficients can be used directly as weighting factors in the least-squares adjustment without propagation to the spatial domain. Here, we estimate these uncertainties of the GRACE/GRACE-FO coefficients from the root-mean squared residual of each coefficient’s time series. We note that the spectral inversion estimates are very similar to previous results performed in the spatial domain (Sasgen et al., 2013; Sasgen et al., 2019b), however, show a slightly lower noise level and better agreement with the expected fluctuations of the surface mass balance. The results fall within the range of a multi-source ensemble of GRACE mass balances for Antarctica undertaken as part of the Ice Sheet Mass Balance Inter-Comparison Exercise (IMBIE2) (Shepherd et al., 2018) (Supplementary Figure S2).

As a final step, we correct the trends of mass change for each basin for the signal induced by glacial-isostatic adjustment (GIA) using the arithmetic average of the models IJ05 R2 (Ivins et al., 2013), AGE1 (Sasgen et al., 2013) and ICE-6G_D (Peltier et al., 2015) (Supplementary Figure S3). These GIA models are a subset of the larger ensemble investigated for the IMBIE2 (Shepherd et al., 2018) and selected as they have shown to reproduce GPS-inferred uplift rates in Antarctica (Martín-Español et al., 2016), except for high uplift rates measured in the Amundsen Sea Embayment region (Barletta et al., 2018). At the same time, the models are constructed using independent approaches and largely different input data (Supplementary Table S2 in Shepherd et al., 2018), allowing us to reflect regional variability in the GIA patterns. In addition, it has been shown that applying this GIA correction produces GRACE mass rates consistent with CryoSat-2 mass rates for the time period 2011–2017 (Sasgen et al., 2019b). As the consistency with GPS-inferred uplift rates implies similar response magnitudes and patterns, the uncertainty of GIA derived from the model differences is smaller compared to the full IMBIE2 ensemble.

We acknowledge two limitations related to our GIA correction. First, regional optimized models with improved fit to the GPS uplift rates exist for the Amundsen Sea Embayment region (Barletta et al., 2018), the Antarctic Peninsula (Nield et al., 2014), and the Siple Coast (Nield et al., 2016), yielding GIA-induced apparent mass changes of 17 Gt yr⁻¹, 3 Gt yr⁻¹, and a range of ±6–8 Gt yr⁻¹, respectively. However, these predictions are based on regional Earth structure models and cannot be simply superimposed with continent-wide GIA simulations adopting a global average Earth structure. Moreover, it has been shown that adopting a thin elastic lithosphere may produce large localized uplift rates without a large impact on the gravity field change at the spatial resolution of GRACE/GRACE-FO (Sasgen et al., 2017). However, we note that mass loss rates in the Amundsen Sea Embayment region could be underestimated by 10 Gt yr⁻¹ to 14 Gt yr⁻¹ Barletta et al. (2018).

Another limitation is the assumption that GIA is approximately linear over short time periods and has henceforth no impact on accelerations measured with GRACE/GRACE-FO. Barletta et al. (2018) have shown that for Amundsen Sea Embayment region, underlain by a low-viscous upper mantle, the GIA rate increases over time. The authors state an increase of the GIA rates of about 2% per year, equivalent to a per year increase of the GIA-induced apparent rate of mass change of +0.3 Gt yr⁻¹. However, the acceleration of GIA in this region remains an issue of investigation, as it strongly depends on the poorly known load changes during the past century. In addition, adding the effect of this local response based on a specific viscosity structure to the continent-wide GIA corrections would introduce an inconsistency we can not resolve within this study.

We estimate the uncertainty of the monthly GRACE/GRACE-FO mass changes at the basin scale from two components; the time varying, normalized noise level of the GRACE/GRACE-FO gravity field solutions, and the residual variability of the mass change within each basin. To obtain the noise level of each monthly solution, we calculate the coefficient’s residual time series   C j m R e s ( t i )   after removing bias, trend, acceleration, annual, and semi-annual and temporal variations longer than 4 months (using a moving average filter). Then, we cumulate the degree power of the residual in the noise-dominated spectral range 40–60, ∑ j = 40 60 ∑ m = − j j C j m R e s ( t i ) . Finally, we normalize the cumulative degree power by the respective mean of all solutions, meaning that low-quality solutions show noise levels larger than one, high-quality solution lower than one. Then, we apply this time dependent noise factor to the root-mean-squared residual variability of mass change within each drainage basin. The resulting monthly uncertainty estimates for the SDS data sets (GFZ RL06, CSR RL06 and JPL RL06) and their combination (AV RL06) used for error propagation when estimating trends and accelerations (Supplementary Figure S4).

Apart from the stochastic noise associated with the GRACE/GRACE-FO coefficients, long-term mass changes may contain systematic differences arising from processing choices, background models and corrections, such as the GIA correction (Martín-Español et al., 2016; Shepherd et al., 2018) or the ocean and atmosphere de-aliasing during Level-2 data processing (Dobslaw et al., 2017; Hardy et al., 2017). Here, we empirically estimate the uncertainty of the trends and accelerations from the differences of our mass change time series with the publicly available Level-3 mass change products TUD COST-G (Sasgen et al., 2019a) and the CSR RL06 Mascons (Save et al., 2016) (Data availability). To increase the sample size, we estimate the linear trend and acceleration by randomly selecting 163 out of 191 monthly solutions from each data set and repeating this procedure 1,000 times in a bootstrap approach. In addition, we estimate the uncertainty of the GIA adjustment correction based on the standard deviation of model realizations entered into IMBIE2 (Shepherd et al., 2018) (Supplementary Figure S3). The total uncertainty of the GRACE/GRACE-FO trend and acceleration components is estimated as the quadratic sum of the respective propagated monthly uncertainties (stochastic) and the cross-solution empirical uncertainties, including uncertainties of the correction for GIA in the linear trend component (systematic). The standard deviations of the linear trends and accelerations and their stochastic and systematic components are shown in Supplementary Figure S5.

We approximate the SMB over the grounded part of the Antarctic Ice Sheet by net accumulation, calculated as total snowfall, minus snowmelt and evaporation based on the re-analyses product ERA-5 of the European Centre for Medium-Range Weather Forecasts (ECMWF) for the time period January 1979 to September 2020 (Copernicus Climate Change Service, 2017). The data are provided on a regular 0.25° × 0.25° geographic grid, yielding a spatial resolution of approximately 15 km along the coast of West Antarctica (latitude ∼73°S). In addition, we use 35 km × 35 km output of regional climate model MARv3.6 (Agosta et al., 2019) for quantifying uncertainties related to the surface mass balance component (see following section). In contrast to the Greenland Ice Sheet, which experiences extensive summer melt, snowfall constitutes almost the entire surface-mass balance over grounded ice in Antarctica. Our analysis of ERA-5 data shows that snowmelt amounts to less than 1% of the regionally averaged surface mass balance, with some exception along the Antarctic Peninsula (up to 4%), confirming previous results (Lenaerts et al., 2012; Agosta et al., 2019). We estimate storage changes as measured by GRACE/GRACE-FO by integrating monthly estimates of SMB, hereafter referred to as cumulative SMB, S M B ( t ) , over the time periods 2002–2017 and 2002–2020, respectively. The cumulative SMB is summed up within each basin or region (Figure 1).

We estimate the uncertainty of cumulative SMB based on four model realizations, which involve the ERA-5 reanalysis and three MARv3.6 simulations with different lateral forcing (ECMWF ERA-Interim, MERRA2 and JRA-55) (Agosta et al., 2019). Thus, the uncertainty derived from this ensemble encompasses deviations in the large-scale atmosphere circulation, as well as inherent modelling differences. We separate the uncertainty analysis into three temporal components, reflecting differences between de-trended monthly estimates, the linear trend and the acceleration of the cumulative SMB estimates, respectively. The differences are analyzed for moving 15-years time windows within the time period covered by all four realizations (1979–2015). While deviations of detrended monthly cumulative SMB estimates tend to reflect random fluctuations introduced by differences in the atmosphere forcing and climate modeling, deviations of long-term changes (e.g., trends and accelerations) will arise from potential biases and different sensitivities in the responses of the models to the same forcing. Thus, similar to the GRACE/GRACE-FO analysis, we consider the random monthly uncertainties of ±55 Gt (Supplementary Figure S4) and their propagated uncertainties of the trends and accelerations (hereafter named stochastic; ±10 Gt yr⁻¹ and ±0.5 Gt yr⁻², respectively). As a second component, we estimate systematic uncertainties of the SMB trends and accelerations, amounting to ±9 Gt yr⁻¹ and ±0.4 Gt yr⁻², respectively (Supplementary Figure S5). The total uncertainty is obtained by quadratic summation of the propagated monthly and the systematic uncertainties.

For each basin, we decompose the time series of mass change from GRACE/GRACE-FO and cumulative SMB using least-squares linear regression with a quadratic temporal model consisting of offset, M ¯ , linear trend, M ˙ , and acceleration, M ¨ , M ^ ( t ) = 1 2 M ¨ t 2 + M ˙ t + M ¯ (2)

In the regression, we employ time variable (stochastic) monthly uncertainties for GRACE/GRACE-FO, which reduces the influence of low-quality monthly solutions, for example at the end of the GRACE mission in 2017 (Supplementary Figure S4). We assume constant uncertainties and therefore equal weighting for the monthly cumulative SMB estimates, in the absence of representative time-variable monthly uncertainties. Note that throughout this paper, accelerations are presented as the temporal derivative of the rates (factor 1/2 in Eq. 2).

The offset M ¯ in Eq. 2 reflects the mean mass distribution on and within the Earth in the region of Antarctica, which is not of interest here. The respective linear trend M ˙ represents the average rate of mass change either of the ice sheet (GRACE/GRACE-FO), the rate of cumulative SMB or the rate of cumulative discharge, S M B ˙ ( t ) and D ˙ ( t ) , respectively. The acceleration term M ¨ in the GRACE/GRACE-FO time series represents, depending on the sign, a long-term decrease or increase of the mass balance of the ice sheet. The acceleration of cumulative SMB, S M B ¨ , indicates trends in the rate of mass input, caused by changes in the rate of snowfall accumulation due to atmospheric variability and climate trends. The residual of the GRACE/GRACE-FO mass change and quadratic fit of Eq. 2, M ( t ) − M ^ ( t ) , is dominated by mass changes related to interannual SMB variability in addition to noise. Although dynamic discharge may show year-to-year variability for individual glaciers and ice streams, these fluctuations are below the detection limit of GRACE/GRACE-FO, more so because regional averaging further reduces the possibly uncorrelated variability within a larger entity.

The comparison of GRACE/GRACE-FO mass changes, M ( t ) , and cumulative SMB mass changes, S M B ( t ) , shows very high correlation and agreement in magnitude at interannual time scales (Supplementary Figure S6), providing justification for calculating discharge by differencing both data sets. To further increase the robustness of our results, we group the 25 basins to ten regions as shown in Figure 1, effectively decreasing the spatial resolution and, thus, noise in the GRACE/GRACE-FO data. For this, we identify blocks of co-varying interannual changes of the cumulative SMB, as reflected by inter-basin correlations (Supplementary Figure S7). These covariations are a consequence of a common regional atmospheric forcing causing synchronous fluctuations and regional multi-year trends in snow accumulation (Mémin et al., 2015), discussed later in this paper.

The mass changes M ( t ) of the Antarctic Ice Sheet result from the difference between the cumulative surface mass balance S M B ( t ) minus the cumulative ice discharge over the grounding line D ( t ) (Rignot et al., 1997, 2008; Shepherd et al., 2018; Mankoff et al., 2020). Here, we reorder the mass balance equation (Eq. 1), to estimate cumulative discharge as D ( t ) = M ( t ) − S M B ( t ) , where M ( t ) is obtained from GRACE/GRACE-FO, and S M B ( t ) is obtained from the snow accumulation product of the re-analysis data (ERA-5). In the following, we focus on the long-term temporal components for specific time periods; namely the rate and acceleration of cumulative discharge, D ˙ and D ¨ , respectively, and likewise its rate and acceleration of cumulative SMB, S M B ˙ and S M B ¨ , respectively.

For fast-flowing ice streams and fixed geometries, D ˙ is proportional to the surface-ice velocity (Rignot et al., 1997, 2008; Mankoff et al., 2020), therefore an increase in surface-ice velocity implies a component of accelerated discharge D ¨ . On the other hand, an acceleration of S M B ( t ) may reflect both long term and temporary trends in the snowfall rate, for example, due to increased moisture content in the atmosphere or prevailing circulation, respectively. We estimate the rates and accelerations independently for the time series of M ( t ) and S M B ( t ) and infer accelerations of D ( t ) by subtraction. Uncertainties of the trends and accelerations of M ( t ) and S M B ( t ) are added to obtain uncertainties of the respective temporal component of D ( t ) (Supplementary Figure S5).

Note that S M B ( t ) is defined as positive if adding mass to the ice sheet system; it has positive values for all basins. Therefore, a positive acceleration of S M B ( t ) denotes an increase in mass gain, indicated by a positive regression coefficient for the acceleration. Likewise, D ( t ) is defined as positive if adding mass to the ice sheet system and has negative values for all basins. In the following, an acceleration of D ( t ) refers to an increase in mass loss out of the system, indicated by a negative regression coefficient for the acceleration. The mass change of the ice sheet M ( t ) , can have either sign. Therefore, we specify an increase in mass loss, as acceleration of mass loss (negative regression coefficient), a reduction in mass loss as deceleration (positive regression coefficient). In this definition, an acceleration in S M B ( t ) can compensate an acceleration in D ( t ) in terms of mass loss M ( t ) .

We compare rates and accelerations of D ( t ) obtained with the indirect method with estimates based on surface-ice velocity measurements (Rignot et al., 2019a), referred to as “direct” estimates. As basin definitions slightly differ between both estimates, we map the direct estimates onto our basins by scaling with the respective area, in the few cases of intersecting basins. For consistency, we limit the evaluation period to the common time interval 2002–2017, for which we calculate the unweighted mean of the annual rates of cumulative discharge from the direct method (Rignot et al., 2019a).

In the following, we will extrapolate sea level rise caused by the recovered mass balance trends and ice-dynamic accelerations to the year 2100 (here; 80 years). As a sensitivity experiment, we distinguish between the extrapolation based on the GRACE/GRACE-FO mass balance M ˙ only, and the extrapolation additionally including acceleration of the cumulative discharge, D ( t ) , denoted D ¨ , obtained by the indirect method. We evaluate the temporal model of Eq. 2 and extrapolate the total mass change to the year 2100, M t ( t ) = 1 2   D ¨ t 2 + M t ˙   (3)

The quadratic extrapolation is supported by Ice Sheet Model Intercomparison 6 (ISMIP6) of the Coupled Model Intercomparison Project Phase 6 (CMIP6) projections of Antarctica’s contribution to sea-level rise until the year 2100, which suggest a time evolution well-approximated by a second-order polynomial, across models and model set ups, as well as using various climate forcing (Payne et al., 2021). Note, however, that despite similar temporal characteristics of the projections, the complex interplay of physical processes in response to the climate forcing governs the actual path of each simulation. We chose to extrapolate to the year 2100 as the typical time frame of ice sheet model projections. In addition, (Lowry et al., 2021), show that ice-sheet model projections forced by RCP2.6 and RCP8.5 scenarios become significantly different only after 2100 (at the 68% confidence level).

We would like to emphasize limitations of this extrapolation. First, mass balance and dynamic acceleration are assumed to prevail unchanged until the end of this century, implying no changes in the external forcing. Next, we neglect possible long-term changes in SMB, which have been shown to have an important impact on the sea-level contribution Payne et al. (2021). In addition, we extrapolate present-day mass balance components based on a relatively short time period of 19 years. Wouters et al. (2013) have estimated the uncertainty in mass trend estimates due to stochastic variability of the Antarctic Ice Sheet as function of observation length. For the contribution of SMB and a 19-years time series, Wouters et al. (2013) estimate a stochastic variability of ±3.4 Gt yr⁻² for the acceleration (denoted according to Eq. 2) and ±3.9 Gt yr⁻¹ for the rates (Figure 3 in Wouters et al., 2013). In addition, the authors suggest a contribution from ice dynamic variability an order of magnitude lower. The uncertainty of the rates (±3.9 Gt yr⁻¹) and the dynamic acceleration (±0.3 Gt yr⁻², denoted according to Eq. 2) caused by stochastic variability is added to the propagated uncertainties of the indirect method in the extrapolation of the sea-level contribution. Note that the uncertainty in the extrapolation is based only on the propagation of the uncertainty of the presented-day trends and accelerations; it does not included process-related uncertainties, which could possibly be estimated in a Monte Carlo approach based on the CMIP6 output. In the following, we assume a volume density of water of 997 kg m⁻³ and the area of the ocean of 361·10⁶ km², such that 360 Gt of ice-mass loss causes 1 mm of global mean sea-level rise.

Regional mass changes and their associated linear trends are estimated based on 191 monthly GRACE and GRACE-FO gravity solutions for the Antarctic Ice Sheet (Figure 2). For the entire ice sheet, rates of mass change, M ˙ , amount to -131 ± 23 Gt yr⁻¹, of which -138 ± 8 Gt yr⁻¹ originate from West Antarctica, -17 ± 9 Gt yr⁻¹ from the Antarctica Peninsula and +24 ± 20 Gt yr^-1 from East Antarctica. Within the study period from April 2002 to September 2020, the ice sheet contributed approximately 6 mm to global mean sea level rise.

The three main subdivisions of the Antarctic Ice Sheet (Figure 1) show different mass evolutions. The year 2009 marks a transition to mass gain in East Antarctica, continuing until 2012, and increased loss rates in West Antarctica (Figure 2). The mass gain in East Antarctica is mainly a consequence of high precipitation events in Dronning Maud Land (Boening et al., 2012) (Supplementary Figure S6), cancelling out the increasing losses in West Antarctica after 2009 in the mass balance of the whole ice sheet. Apart from long-term trends, all regions exhibit substantial interannual fluctuations (Supplementary Figure S6). East Antarctica contributes the most to the interannual fluctuations of the entire ice sheet, which appear random, however, reflect the superposition of asynchronous mass signals arising from accumulation variability within smaller regions of this largest part of Antarctica, discussed later. In contrast, the mass balance in West Antarctica dominates the mass loss signal of the entire Antarctic Ice Sheet (-138 ± 23 Gt yr⁻¹ versus -131 ± 8 Gt yr⁻¹, respectively), with an additional contribution from the Antarctic Peninsula.

Next, we investigate the interannual variations of the mass changes M ( t ) from GRACE/GRACE-FO and S M B ( t ) , obtained as residual after removal of the fitted temporal model M ^ ( t ) shown in Eq. 2 (offset, linear and quadratic trends). For the entire Antarctic Ice Sheet M ( t ) − M ^ ( t ) from GRACE/GRACE-FO and the cumulative SMB show a Pearson correlation coefficient of 0.59, well above the 95% significance level of 0.30 accounting for auto-correlation, which reduces the effective degrees of freedom (Barker et al., 2016). Applying a three-month moving average filter reduces short-term GRACE/GRACE-FO monthly fluctuations, which are not present in the SMB time series and therefore likely associated with noise. This increases the correlation to 0.72 at a significance level of 0.40 (including autocorrelation). After application of the moving average, the temporal characteristics of M ( t ) from GRACE/GRACE-FO and S M B ( t ) are highly consistent (Supplementary Figure S6); therefore, we refrain from similar filtering of the SMB time series.

At the basin scale, including the 3-month moving average (Supplementary Figure S6), correlations lie above 0.8 for twelve of the 25 basins. These basins show a large magnitude of the SMB variability with respect to GRACE/GRACE-FO uncertainties. We identify two regions where de-trended GRACE/GRACE-FO mass changes correlate remarkably well with SMB fluctuations; basins 20 through 24 in the Amundsen and Bellingshausen Sea Embayment region (coefficient >0.92) and basins 4 through 8 in Dronning Maud Land (coefficient >0.85). Highly correlated, with coefficients above 0.83, are basins 11 through 13 in East Antarctica, while remaining basins show moderate correlation. For example coefficients for basins 14 through 17, adjacent to the Ross Sea, range between 0.38 and 0.64. A likely reason is that magnitudes of variation in SMB ( t ) are lower in these basins, leading to a smaller signal-to-noise ratio in the GRACE/GRACE-FO time series. Similar conditions may apply to basins 9 and 10 (Amery Ice Shelf region), even though here, correlations are typically greater than 0.62. We note that without filtering, correlations decrease on average by 0.07.

Rates of cumulative discharge, D ˙ , obtained for 2002–2017 are consistent between the indirect and the direct method at the regional scale (Figure 3); for the entire Antarctic Ice Sheet and four out of ten regions, values agree within their respective uncertainties. However, in six regions, deviations exceed the uncertainty ranges, possibly due to underestimated uncertainties of the SMB component in the indirect method. The median relative uncertainty for all regions is 5% for the SMB component, adding up to a relative uncertainty of the rates of cumulative discharge of 7% for the indirect method, compared to approximately 1% stated for the direct method.

For 2002–2020, we obtain an average annual surface mass balance, S M B ˙ based on ERA-5 snow accumulation minus snow melt and snow evaporation of 2228 ± 44 Gt yr⁻¹, comparing well with published values of 2098 ± 133 Gt yr⁻¹ for the time period 1979 to 2008 using a different regional climate model (Rignot et al., 2019a). Note that a two to three times larger relative uncertainty of SMB (10–15%; equal to ±88 to ±132 Gt yr⁻¹ for the entire ice sheet, respectively) would reconcile the indirect and direct discharge estimates within uncertainties for all regions, as would seven times larger discharge uncertainties of the direct method (7%; equal to ±62 Gt yr⁻¹).

In the following, we present basin-scale estimates of the acceleration of the mass change from GRACE/GRACE-FO, S M B ( t ) approximated by ERA-5 net snow accumulation, as well as D ( t ) calculated as difference of both observations for the time period 2002–2020. The spatial distribution of accelerations in each component ( M ¨ ,   S M B ¨ and D ¨ ) is shown in Figure 4 (2002–2020), displayed as average specific mass change within each basin (mass change divided by basin area), enabling a size-independent comparison of imbalances between basins.

Concerning rates of mass change, M ˙ , our estimates based the spectral inversion approach using combined GRACE/GRACE-FO solutions (Supplementary Figure S1) are in agreement with those based on published mass change products using different inversion approaches and input solutions. For example, the total mass balance of -131 ± 23 Gt yr⁻¹ is in good agreement with the independent estimates based on gridded products; -116 ± 23 Gt yr⁻¹ for CSR RL06 Mascons (Save et al., 2016) and -123 ± 23 Gt yr⁻¹ for TUD COST-G (Sasgen et al., 2019a). Furthermore, our mass balance falls within the range of gravimetric estimates from IMBIE2 (Shepherd et al., 2018) for the coincident time interval 2003–2016, for example -123 ± 23 Gt yr⁻¹ (this study) compared to -114 ± 33 Gt yr⁻¹ (IMBIE2) for the entire ice sheet (Supplementary Figure S4). Our results are more negative than a recent estimate of the mass balance of -99 ± 57 Gt yr⁻¹ for April 2002 to September 2019 (Velicogna et al., 2020), yet show agreement within uncertainties. Note that neither of these studies consider the regional GIA response in the Amundsen Sea Embayment region presented by Barletta et al. (2018) thus, possibly underestimating mass loss by 10 Gt yr⁻¹ to 14 Gt yr⁻¹.

To show the impact of atmospheric circulation patterns in the Southern Hemisphere on the interannual mass changes presented earlier, we now calculate temporal correlations with climate indices. In particular, we investigate the Southern Oscillation Index (SOI) (Parker, 1983) and the Southern Annular Mode (SAM) (Thompson, 2002), reflecting atmosphere pressure differences, as well as the NIÑO 3.4/4 indices, derived from gradients in sea-surface temperatures in the central Pacific (Trenberth, 1975). Previous studies have shown correlation of these indices with monthly rates of precipitation in Antarctica (Turner, 2004). As we analyze storage changes of the ice sheet represented by cumulative sums of monthly fluxes, we cumulate the indices analogous to the processing of the SMB data to enable a direct comparison.

The cumulative SMB, S M B ( t ) , and the cumulative indices exhibit groups of basins with correlated variations, similar to the mass changes between basins presented earlier (Supplementary Figure S7). Correlations cluster as blocks of negative and positive values. For example, the correlation with the SOI is negative for the Filchner-Ronne Ice Shelf region (basins 1 and 2), positive for the Dronning Maud Land region (basins 3 through 6), negative for the George V and Oates Land region (basins 13 and 14) and the Siple Coast eastward up to and including the Amundsen Sea Embayment region (basins 18 through 23). These correlations are above the 95% significance level (considering auto-correlation). The indices of NIÑO 3.4/4 based on sea-surface temperatures show largely the same pattern of covariance, however, with opposite sign owing to the strong anti-correlation of SOI and NIÑO 3.4/4 (Supplementary Figure S7) (Cullather et al., 1996). Fluctuations in relation to the SAM are visible mainly for central Dronning Maud Land, large parts of the ice sheet draining into the Amundsen Sea Embayment and the northern Antarctic Peninsula (Supplementary Figure S7). Although these correlations with individual indices are significant, we note that their superposition may describe the complex modes of atmospheric circulation and their effect on the SMB variability more accurately (Fogt and Bromwich, 2006; Mémin et al., 2015; Kaitheri et al., 2021). Furthermore, Mémin et al. (2015) and Kaitheri et al. (2021) suggest additional modulation of SMB anomalies related to the Antarctic Circumpolar Wave, propagating the climate anomalies eastward around Antarctica at a period of six to eight years.

The results show that interannual variations of S M B ( t ) and their correlations with atmospheric conditions are an important non-stationary temporal signal in the GRACE/GRACE-FO time series. In particular, the asynchronous modulation of snowfall anomalies between Dronning Maud Land and the Amundsen Sea Embayment region (Supplementary Figure S6 and Supplementary Figure S7) is part of global atmosphere variability closely linked to the El Niño Southern Oscillation (Zhang et al., 2021). For short time scales, changes in rates of mass change between East and West Antarctica are considerably driven by the SMB component (for example in 2009; Figure 2) and should not be misinterpreted as a persistent shift in ice dynamics.

For the entire Antarctic Ice Sheet, as well as for eight out of the ten regions, accelerations of cumulative discharge obtained in this study show agreement with the independent direct estimates within uncertainties between 2002 and 2017 (Rignot et al., 2019b) (Figure 5). Largest values of acceleration of D ( t ) (increasing mass loss) are detected for the Amundsen Sea Embayment region, amounting to a value of D ¨ of -5.0 ± 0.8 Gt yr⁻², well in line with the direct estimates. Exceptions from the overall agreement of direct and indirect estimates are Siple Coast and Edward VII Land, for which we estimate a deceleration D ¨ of +0.7 ± 0.3 Gt yr⁻², somewhat larger than that obtained from the direct approach (+0.3 ± 0.3 Gt yr⁻²). For the Filchner-Ronne Ice Shelf region, we detect an acceleration of c D ( t ) of -1.0 ± 0.5 Gt yr⁻² not present in the direct estimates (+0.1 ± 0.4 Gt yr⁻²). For both regions, cumulative SMB does not remove the decelerations and accelerations prominent in the GRACE data, leading to non-zero values for cumulative discharge in 2002–2017.

We further investigate the robustness of our approach to recover the quasi-stationary accelerations of D ( t ) by prolonging the time period to 2020, now including GRACE-FO data. Figure 5 illustrates that values of D ¨ only moderately change within their uncertainties, when including about three more years of data (2002–2020). In contrast to the stability of the discharge estimates, accelerations of GRACE/GRACE-FO mass change M ( t ) and, prominently, S M B ( t ) show considerable variations when using different time periods. For example, GRACE/GRACE-FO recover considerably different magnitudes of accelerating mass loss in the Amundsen Sea Embayment region in both periods. For the SMB components, accelerations estimates S M B ¨ even differ in sign. Similarly, for the entire ice sheet, acceleration of GRACE/GRACE-FO mass loss is larger in 2002–2017 compared to 2002–2020, while the estimated acceleration value D ¨ remains mostly unchanged, due to subtraction of different SMB contributions. This further supports our finding that it is possible to separate slow-varying discharge changes from SMB-induced fluctuations in the GRACE/GRACE-FO data.

In contrast to the uncertainties of the rates of D ( t ) discussed earlier (Figure 3), uncertainties of the acceleration estimates D ¨ are of similar magnitude for the direct and indirect methods. While possible biases in the absolute amount of snow accumulation limit the precision of the indirect discharge rates D ˙ , as they determine the mean rate of S M B ( t ) , accelerations are influenced by uncertainties in the temporal trends of snow accumulation. In addition, for the discharge rates D ˙ , the GIA correction of the GRACE/GRACE-FO is an important source of systematic uncertainty (Supplementary Figure S5), while it is much less important when considering accelerations D ¨ .

Our study locates the largest acceleration of D ( t ) (increase in mass loss) in the Amundsen Sea Embayment region (Figure 5), in accordance with prior studies (Sutterley et al., 2014; Gardner et al., 2018) and analysis using direct observations (Rignot et al., 2019a). In this region, various observations indicate that the ice sheet is in a state of transition, with the risk of exceeding the instability threshold related to the Marine Ice Sheet Instability (DeConto and Pollard, 2016; Sutter et al., 2016).

For many glacier systems, retreating and thinning ice shelves have triggered fundamental ice-dynamic responses. For example, the mass loss of the Getz Ice Shelf (adjacent to basin 20) was followed by an increase in the velocity of the tributary glaciers in the period of 2007–2013, as well as dynamic glacier thinning and grounding line retreat between 2010 and 2013 (Joughin et al., 2014; Rignot et al., 2014; Chuter et al., 2017; Gardner et al., 2018). An increase in ice discharge by 6% is determined for the Pine Island, Thwaites, Haynes, Pope, Smith and Kohler Glaciers between 2008 and 2015, which equals an acceleration of -2.4 Gt yr⁻² (Gardner et al., 2018), somewhat lower than in this study (basins 21 and 22; -3.1 ± 0.7 Gt yr⁻²). Overall, our estimates of cumulative discharge acceleration are in line with the increase of 77% in this sector from 1973 to 2013 (Rignot et al., 2019a).

In basins 1 and 2 of the Filchner-Ronne Ice Shelf region, we detect moderate acceleration of D ( t ) between 2002 and 2017 (-0.9 ± 0.5 Gt yr⁻²; Figure 5), in contrast to the direct estimate of +0.1 ± 0.4 Gt yr⁻² (Rignot et al., 2019a). Although accelerations remain negative using GRACE/GRACE-FO data of CSR RL06 Mascons or TUD COST-G, we cannot exclude possible biases in the gravimetry data caused by long-term tidal-aliasing signals, known to be a source of uncertainty in the area of the large ice shelves in Antarctica (Dobslaw et al., 2017).

For the Antarctic Peninsula, our results show no significant change in ice-dynamic discharge, which is in contrast to the direct estimates, indicating moderate acceleration (Figure 5). At the basin scale, we identify discharge deceleration (+1.3 ± 0.4 Gt yr⁻²) in northern and acceleration (-1.2 ± 0.7 Gt yr⁻²) in the southern part of the Antarctic Peninsula (Figure 4). Evaluating GRACE/GRACE-FO data over this elongated and narrow region is challenging and may be subject to biases introduced by processing choices, even though our analysis using published mass change products CSR RL06 Mascons and TUD COST-G yields consistent results. In addition, estimates of the SMB component could be biased due to the limited resolution of the ERA-5 re-analysis data. This could be clarified in future studies using higher resolution products. However, projections suggest that glacier discharge is substantially less sensitive to collapse of the Larsen C Ice Shelf along the northern Peninsula (basin 25) than the southern George IV Ice Shelf (basin 24) (Schannwell et al., 2018). Less frequent events of ice shelf disintegration compared to the 1990s are possibly related to incipient near-surface temperature cooling (Turner et al., 2016), even though most observations indicate increasing discharge for both basins (Wouters et al., 2015; Cook et al., 2016; Hogg et al., 2017).

In contrast to most of West Antarctica, cumulative discharge decelerates along the Siple Coast and Edward VII Land (basins 18 and 19; +0.7 ± 0.3 Gt yr⁻²) feeding the Ross Ice Shelf. During the last 4 decades, glaciers along the Siple Coast showed reduced ice flow (Rignot et al., 2019a) and experienced mass gain between 1992 and 2017 (Schröder et al., 2019). The Whillans Ice Stream is anticipated to slow down (Joughin et al., 2005), and stagnate by 2070 (Thomas et al., 2013), in-line with a discharge deceleration of 25.3% during 1997 and 2009, when including the adjacent Mercer Ice Stream (Gardner et al., 2018). The reduction in ice velocity is not only present near the grounding line, where the mass loss finally occurs, but reaches a few hundred kilometers inland (Scheuchl et al., 2012).

The Ross Ice Shelf region in East Antarctica (basins 15 through 17) include complex terrain of the Transantarctic Mountains (Fretwell et al., 2013). For this region, our estimates of acceleration D ¨ are close to zero, indicating ice-dynamic stability, enabled by a topography of broad and stabilizing ridges (Morlighem et al., 2020).

George V and Oates Land (basins 13 and 14) show moderate dynamic deceleration, D ¨ , after removal of large SMB-driven accelerations from the GRACE/GRACE-FO data (+1.0 ± 0.9 Gt yr⁻²; Figure 5). Our estimate are in line with previous studies, finding no pronounced glacier retreat between 1972 and 2013 for both basins (Lovell et al., 2017).

For Wilkes Land (basins 11 and 12), East Antarctica, we identify no significant acceleration of D ( t ) , despite a mass loss increase of -5.8 ± 0.2 Gt yr⁻² observed by GRACE/GRACE-FO and mean rates of cumulative discharge being among the largest in Antarctica (2002–2020: -465 ± 22 Gt yr⁻¹; Figure 3). Mass and elevation change within basin 12 (Figure 1), including Totten Glacier, are well documented in previous studies (Schröder et al., 2019), however reflect spatially confined changes close to the coast (Shepherd et al., 2019). Acceleration of M ( t ) (increase in mass loss) detected with GRACE/GRACE-FO is almost completely caused by a decline in snow accumulation ( S M B ¨ is -6.3 ± 1.2 Gt yr⁻² for 2002–2020), as part of a strong interannual variation in this region. In contrast to our estimate, studies have shown a 40% increase of dynamic discharge in 2017 compared to 1979–2003 (Rignot et al., 2019a), as well as a loss of grounded ice (Smith et al., 2020). The future evolution of this part of the ice sheet critically depends on the evolution of the Totten Glacier, due to its large catchment area and retrograde bedrock slope, promoting instability similar to West Antarctica (Gomez et al., 2010; Durand and Pattyn, 2015; Greenbaum et al., 2015). Comparable conditions apply to the Denman Glacier (Brancato et al., 2020; Smith et al., 2020; Miles et al., 2021). It is the consensus that the evolution of heat transport due to changing ocean currents and ocean temperatures will have a decisive impact on this marine-based part of the ice sheet (Rintoul et al., 2016; Levermann et al., 2020). It will be important to reconcile the observational data sets in order to reliably quantify the dynamic changes in this region, holding high potential of sea-level rise.

The Amery Ice Shelf region exhibits only small values of dynamic acceleration D ¨ (+0.5 ± 0.2 Gt yr⁻²), but also of S M B ¨ (-0.7 ± 0.2 Gt yr⁻²), indicating relative dynamic stability and a low influence of mass change by atmospheric circulation, particularly in the interior part (basin 10). Our finding is supported by unchanging ice velocities since 1970 (Rignot et al., 2019a) and observation of slight deceleration in ice flow from 1968 to 2014 (Tong et al., 2018).

Enderby Land and Kemp Land (basins 7 and 8) show no sign of acceleration of cumulative discharge in the investigation period, in agreement with direct observations (Figure 5). However, for Dronning Maud Land (basins 3 through 6) we isolate an acceleration of -1.5 ± 0.8 Gt yr⁻² for D ( t ) , which is similar for both time periods investigated, despite a large variation in the SMB component (Figure 5). We note that mass changes obtained from altimetry (Schröder et al., 2019) are consistent with GRACE/GRACE-FO. Therefore, we expect accelerations to be of similar magnitude and sign using altimetry instead of gravimetry data.

Next we extrapolate the present-day mass balance (-131 ± 23 Gt yr⁻¹) together with the acceleration of D ( t ) (-5.3 ± 2.2 Gt yr⁻²) for a time period of 80 years, to explore the relevance of the quadratic term of today’s mass balance for the future sea-level rise contribution from Antarctica. Assuming 360 Gt of mass loss to be equivalent to 1 mm of sea-level rise and applying Eq. 3, we obtain a sea-level contribution of 7.6 ± 2.9 cm by the year 2100. Note that we have included an additional uncertainty of ±0.4 Gt yr⁻² associated with the stochastic variably of ice discharge, as inferred by Wouters et al. (2013). This sea-level contribution is more than two times the amount obtained by extrapolating only the linear trend of the mass balance (2.9 ± 0.5 cm), highlighting the sensitivity of the future sea-level contributions to ice-dynamic accelerations. Adopting non-linear GIA correction for the GRACE/GRACE-FO data of +0.6 Gt yr⁻² (Barletta et al., 2018) would further increase the acceleration of cumulative discharge to -6.0 ± 2.2 Gt yr⁻², resulting in a sea-level contribution of 8.0 ± 2.9 cm by 2100.

The extrapolation including accelerations yields with 7.6 ± 2.9 cm a sea-level contribution that lies within the medium range of -8 and 30 cm of CMIP5/ISMIP6 projections for unmitigated climate change (Seroussi et al., 2020) and the CMIP6/ISMIP6 ensemble evaluation of Edwards et al. (2021), projecting likely contribution between -5 and 14 cm by year 2100. Lowry et al. (2021) show that limiting projections to trajectories satisfying present-day observations considerably reduces the ensemble spread caused by uncertainties in model parameters, for example related to ice rheology and sliding. As the uncertainty of the present-day observations expressed as sea-level rise in 2100 (±2.9 cm) is considerably lower than the ensemble spread, we expect that evaluating the projection ensembles with observed trends and accelerations will further confine the range of likely trajectories.

We emphasize the severe limitation that our extrapolation assumes constant rate of S M B ( t ) until the year 2100. Projections have shown that long-term changes in snow accumulation related to a higher moisture content in the warming atmosphere may compensate dynamic losses, potentially leading to net mass gain of the ice sheet and an associated fall in sea level (Seroussi et al., 2020). However, the spread in projected SMB is large due to differences in forcing, models and downscaling (Payne et al., 2021), and further research is needed to quantify the relative importance of SMB, ice-dynamic changes and their interaction in projections of sea-level rise. Our results show that it is possible to separate superimposed substantial multi-year fluctuations in the SMB component from ice-dynamic changes in GRACE/GRACE-FO data. We consider such an isolation in current satellite data an essential first step for model validation and calibration, in order to reduce uncertainties in projections of the future sea-level contribution from the Antarctic Ice Sheet.