# Exploring the uncertainty in GRACE estimates of the mass redistributions at the Earth surface: implications for the global water and sea level budgets

Alejandro Blazquez, B. Meyssignac, Jm Lemoine, Etienne Berthier, A Ribes, A. Cazenave

► To cite this version:

Alejandro Blazquez, B. Meyssignac, Jm Lemoine, Etienne Berthier, A Ribes, et al.. Exploring the uncertainty in GRACE estimates of the mass redistributions at the Earth surface: implications for the global water and sea level budgets. Geophysical Journal International, 2018, 215 (1), pp.415-430. 10.1093/gji/ggy293. hal-01952652

HAL Id: hal-01952652

https://hal.science/hal-01952652

Submitted on 25 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Exploring the uncertainty in GRACE estimates of the mass redistributions at the Earth surface: implications for the global water and sea level budgets

A. Blazquez,<fnref n="1" />,<fnref n="2" /> B. Meyssignac,<fnref n="1" />,<fnref n="2" /> J.M. Lemoine,<fnref n="2" />,<fnref n="3" />,<fnref n="4" /> E. Berthier,<fnref n="1" /> A. Ribes<fnref n="5" /> and A. Cazenave<fnref n="1" />

Accepted 2018 July 19. Received 2018 July 10; in original form 2017 November 08

## SUMMARY

Observations from the Gravity Recovery and Climate Experiment (GRACE) satellite mission provide quantitative estimates of the global water budget components. However, these estimates are uncertain as they show discrepancies when different parameters are used in the processing of the GRACE data. We examine trends in ocean mass, ice loss from Antarctica, Greenland, arctic islands and trends in water storage over land and glaciers from GRACE data (2005–2015) and explore the associated uncertainty. We consider variations in six different GRACE processing parameters, namely the processing centre of the raw GRACE solutions, the geocentre motion, the Earth oblateness, the filtering, the leakage correction and the glacial isostatic adjustment (GIA). Considering all possible combinations of the different processing parameters leads to an ensemble of 1500 post-processed GRACE solutions, which is assumed to cover a significant part of the uncertainty range of GRACE estimates. The ensemble-mean trend in all global water budget components agree within uncertainties with previous estimates based on different sources of observations. The uncertainty in the global water budget is ±0.27 mm yr⁻<fnref n="1" /> [at the 90 per cent confidence level (CL)] over 2005–2015. We find that the uncertainty in the geocentre motion and GIA corrections dominate the uncertainty in GRACE estimate of the global water budget. Their contribution to the uncertainty in GRACE estimate is respectively ±0.21 and ±0.12 mm yr⁻<fnref n="1" /> (90 per cent CL). This uncertainty in GRACE estimate implies an uncertainty in the net warming of the ocean and the Earth energy budget of ±0.25 W m⁻<fnref n="2" /> (90 per cent CL) when inferred using the sea level budget approach.

**Key words:** Global change from geodesy; Satellite geodesy; Sea level change; Time variable gravity.

## 1 INTRODUCTION

The Gravity Recovery and Climate Experiment (GRACE) mission sponsored by the National Aeronautics and Space Administration and the Deutsches Zentrum fur Luft-und Raumfahrt has been providing precise, time-varying measurements of the Earth’s gravitational field since 2002 (Tapley et al. 2004). Changes in the Earth gravitational field are caused by changes in the mass distribution in the solid Earth and at the Earth’s surface, i.e. in the ocean, the atmosphere and on land. The solid Earth processes that can cause significant variations in the Earth gravity field at interannual to decadal timescales include essentially the Earth and Ocean-load tides, the solid Earth pole tide, the glacial isostatic adjustment (GIA) and

the earthquakes (at small spatial scale). Once corrected for these effects, the variations in the time-varying gravity field represent the water mass exchanges at the Earth surface within and among the ocean, the atmosphere and the water/snow/ice storage on land. By providing the first global, satellite-based, accurate measurement of these exchanges on a monthly basis, the GRACE mission has given unprecedented insights on the Earth water cycle changes in response to the current climate change (Wouters et al. 2014).

In particular GRACE measurements have provided since 2002 estimates of many components of the global water budget. These components are the ice loss from glaciers (Jacob et al. 2012; Gardner et al. 2013; Schrama et al. 2014), ice sheets (Shepherd et al. 2012; Velicogna et al. 2014), the water storage changes in major river

<footnote n="1">¹LEGOS, Université de Toulouse, CNES, CNRS, UPS, IRD, 31400 Toulouse, France. E-mail: alejandro.blazquez@cnes.fr</footnote>

<footnote n="2">²CNES, 31400 Toulouse, France</footnote>

<footnote n="3">³GET, Université de Toulouse, CNRS, IRD, UPS, 31400 Toulouse, France</footnote>

<footnote n="4">⁴GRGS, 31400 Toulouse, France</footnote>

<footnote n="5">⁵CNRM, Météo France – CNRS, 31100 Toulouse, France</footnote>
basins (e.g. Llovel *et al.* 2011; Chen *et al.* 2015a; Reager *et al.* 2016) and the variations of the ocean mass due to transfer of water between continents and the ocean (Chambers 2009; Leuliette & Willis 2011). GRACE estimates further helped in identifying and unraveling the role of each of the global water budget components in the contemporary sea level rise (Yi *et al.* 2015; Dieng *et al.* 2015a; Reager *et al.* 2016; Rietbroek *et al.* 2016).

An interesting feature of GRACE-based estimates of the different components of the global water budget is that they are derived from the same single observing system and thus show a high level of consistency with each other when computed with the same global GRACE solution. For this reason, GRACE solutions have provided essential and critical observations to analyse and test the closure of the global water budget (Church *et al.* 2013; Llovel *et al.* 2014; Yi *et al.* 2015; Reager *et al.* 2016; Rietbroek *et al.* 2016; Dieng *et al.* 2017). However, previous studies have noted significant differences in the water budget components estimates when different GRACE solutions from different data processing centres are used or when different post-processing is applied to the data (e.g. Gardner *et al.* 2013, for glacier mass changes; Barletta *et al.* 2013, Velicogna & Wahr 2013 for ice sheets mass changes; Quinn & Ponte 2010, Chambers & Bonin 2012, for ocean mass changes and Reager *et al.* 2016; for the terrestrial water storage changes). These results indicate that each single post-processed GRACE solution provides consistent estimates of the global water budget components (in the sense that they close the global water budget), but they are potentially biased. Thus if GRACE solutions are considered alone, they likely underestimate the true uncertainty in global water budget components.

In this study, we propose to estimate the components of the global water budget from GRACE in terms of trends over the period of January 2005 to December 2015 and to evaluate the associated uncertainty using an ensemble of global GRACE solutions and an ensemble of post-processing parameters. We break down the uncertainty into uncertainties associated to each post-processing parameter. Here, we consider the different component of the global water budget as global ocean mass, Greenland (ice sheet + peripheral glaciers), Antarctica (ice sheet + peripheral glaciers), glaciers on arctic islands while others glacier are lumped with terrestrial water storage. Glaciers on arctic islands are separated from other glaciers as they represent a coherent component for which several mass change estimates have been published in the recent past.

We analyse five GRACE spherical harmonics solutions from five different processing centres: the Center for Space Research (CSR), the Jet Propulsion Laboratory (JPL), the Deutsches Geo-ForschungsZentrum (GFZ), the Technische Universität Graz (TUG) and the Groupe de Recherche en Géodésie Spatiale (GRGS). These solutions cannot be directly used to evaluate the global water budget components; they need to be post-processed first (Wahr *et al.* 1998). The post-processing parameters include (i) the addition of independent estimates of the geocentre motion as these harmonics are not observable by GRACE (ii) the substitution of the Earth oblateness by independent estimates as this harmonic is poorly observed by GRACE (iii) a filtering for correlated errors that map into characteristic north-south stripes, (iv) a correction for the large land signals (from hydrology, glaciers and ice sheet) that can ‘leak’ into the ocean because of the coarse spatial resolution of GRACE, and (v) a correction for GIA. For each GRACE solution from our ensemble we test a range of post-processing parameters to get a spread of estimates of the ice sheet mass changes, the continental water storage changes and the ocean mass changes. Our selection of GRACE solutions and post-processing parameters is not exhaustive

and thus does not cover the whole range of uncertainty of GRACE. However, because we choose the most up-to-date and the largest possible number of parameters, the spread of our ensemble should approach the real underlying uncertainty. The implications on the closure of the global water budget, the sea level budget and the Earth energy budgets are further analysed.

Previous studies have used a similar ensemble approach based on different GRACE solutions and different post-processing to assess the uncertainty in GRACE estimates but they all focused on a single component of the water budget such as the ocean mass change (Quinn & Ponte 2010; Chambers & Bonin 2012), the ice sheet mass changes (Velicogna & Wahr 2013; Chen *et al.* 2015b) or the glacier mass changes (Gardner *et al.* 2013). Here we analyse all the components together in a consistent way. This novel approach enables us to explore whether the uncertainty in the different components of the global water budget are correlated (or not) when assessed with GRACE measurements. This issue is essential when assessing the closure of the global water budget and the sea level budget (the latter includes all components of the water budget plus the thermal expansion of the ocean).

In the second section of this paper we present the five GRACE spherical harmonic solutions that are used in this study and we briefly explain the main differences among them (Section 2.1). We also describe the range of post-processing parameters that are applied to the GRACE solutions (Sections 2.2–2.5) and the statistical framework used to evaluate the uncertainties (Section 2.6). In Section 3 we present an ensemble of GRACE post-processed solutions based on all possible combinations from the set of post-processing parameters and the set of GRACE solutions. From this ensemble, we compute an estimate of the global water budget components and test its sensitivity to the post-processing parameters and GRACE solutions. In Section 4 we compare the ensemble of solutions with the three mascons solutions available from CSR, JPL and GSFC. We propose a new estimate of the uncertainties associated to the water budget components. We also discuss the implications on the closure of the global water budget and sea level budget and on the indirect estimate of the ocean warming and Earth energy imbalance through the sea level budget approach.

In addition to the spherical harmonic solutions, we analyse three mascons solutions from JPL, CSR and the Goddard Space Flight Center (GSFC). These solutions are already post-processed to add the geocentre and the earth oblateness and to correct for the GIA signal. This set of post-processing parameters is the same for the three mascon solutions. We decided not to add these three mascon solutions to our ensemble because it would have resulted in arbitrary larger weight for this set of post-processing parameters when we compute ensemble means. We prefer to use the mascon solutions for comparison (see Section 4.3).

# 2 DATA AND METHODS

## 2.1 GRACE data

We focus on global solutions that are provided to users in the form of spherical harmonic gravity coefficients (Stokes coefficients). Five global solutions are obtained from five different processing centres: CSR, GFZ, GRGS, JPL and TUG. We use the release five of CSR, GFZ and JPL solutions and the release ITSG 2016 of the TUG solution (Klinger *et al.* 2016) from the International center for global Earth models (ICGEM). We use the release 3.3 of the GRGS solution (Lemoine *et al.* 2016). All centres process the same
raw data (level L1) from the GRACE mission which includes the intersatellite range and range rates measured by the K-Band Range (KBR), the range and phase measurements between satellites and the GPS constellation, the accelerometers and the star trackers data. The GRGS centre also uses the satellite laser ranging (SLR) data from LAGEOS1-2, Starlette and Stella in a joint inversion with the GRACE data. To perform the inversion, which yields the gravity field, processing centres use as initial guess an *a priori* model of the gravity field (which is usually a mean field with time-variable components in the form of drift and periodic terms, based on an earlier release of their monthly solutions). In order to reduce the noise, they use ocean and atmosphere models (and sometimes hydrological models) to account for sub-month anomalies that would otherwise alias into the GRACE solutions. They may also modify the weighting of the GPS and the misalignment of the KBR antennas. During the inversion procedure different processing centres use different assumptions, different parameters and different models, which lead to differences among the solutions. Of particular note, the JPL version was considered as a validation product and is by design processed differently than the official CSR and GFZ products. Except for the GRGS product, all processing centres deliver unconstrained solutions, which are not usable without a filtering process (see in Section 2.3). The GRGS product uses an inversion scheme that allows to control the noise of the solutions at small spatial scales (by a normalization of the higher degrees) and therefore does not need any *a posteriori* filtering (Lemoine et al. 2016). The TUG product uses an empirical covariance function of the KBR range-rate data in order to better decorrelate the KBR measurements.

For all GRACE solutions we use the spherical harmonic coefficients up to degree 60 for the conversion to gridded mass anomalies. CSR and TUG solutions provide directly a 60-degree solution, which we use here. For the solutions from other centres, we truncate at the degree 60 to keep the comparison consistent among solutions. To account for the full-mass variability estimated by GRACE, the ocean and atmosphere background models (i.e. initial guess *a priori* models of the gravity field) has to be restored to GRACE solutions. Here we are not interested in the atmospheric mass variability so we only restore the ocean background model and the atmospheric load over the ocean (GAD products, Flechtner et al. 2015) so that we get a full-mass variability estimated by GRACE without the atmospheric variations over land. Over the ocean, the atmospheric model is restored in order to correct the ocean mass variations for the inverse barometer effect and make it comparable with satellite altimetry (of course, for this purpose we remove the time-variable spatial mean of the atmospheric load over the ocean before restoring the atmospheric model). A problem with such an approach is that the atmospheric and oceanic background models are not the same for every GRACE solution (GRGS uses TUGO for the ocean and ERA Interim for the atmosphere while the other centres use OMCT for the ocean and ECMWF's IFS for the atmosphere). To address this issue we use an alternative approach for the GRGS solution: we first restore its own ocean and atmosphere background models to get the estimate of the full mass variability; then we remove the atmospheric load over land and the mean of the atmospheric load over the ocean using the ECMWF's IFS model.

GRACE solutions must be corrected for the pole tide (the solid Earth and ocean response to the polar motion) to get estimates of the surface mass variations. All processing centres except GFZ, use the International Earth Rotation and Reference System Service (IERS) recommended pole tide correction (Altamimi et al. 2016).

But this correction accounts only for the timescales at the Chandler period, since it involves the offset between the instantaneous pole and the mean pole. It does not include the long-period pole tide signals (interannual and longer periods, Wahr et al. 2015). To remove these long-period signals that contaminate GRACE surface mass estimates, we apply to all GRACE solutions the correction proposed by Wahr et al. (2015). For the GFZ solution we also apply the IERS-recommended pole tide correction (which is not applied *a priori*) to get a full pole tide correction and be consistent with other solutions. The mascon solutions that are additionally used in this study are described in Section 4.3. A time mean over 2005–2015 is removed from all GRACE solutions to compute anomalies.

## 2.2 Geocentre motion and Earth oblateness

Several corrections need to be made to the GRACE data in order to be usable over land and over the ocean. The first correction is to include estimates of the degree 1 spherical harmonic coefficients. The degree 1 terms are proportional to the position of the geocentre, defined as the position of the centre of mass of the Earth (CM), relative to the centre of figure (CF) of the Earth's outer solid surface (i.e. the Earth surface over land). As the GRACE mission orbits around the CM and measures the intersatellite range and range rates it cannot observe the degree 1 term. Previous studies based on SLR (Watkins & Eanes 1997; Cretaux et al. 2002; Cheng et al. 2013b; Lemoine & Reinquin 2017), GPS (Fritsche et al. 2010; Wu Xiaoping 2010), ocean models (Swenson et al. 2008) and joint inversions using GPS, GRACE and models (Rietbroek et al. 2012; Wu et al. 2017) provide some estimates of the degree 1 correction. In this study, we consider the four most recent estimates of the degree 1 (Rietbroek et al. 2012; Cheng et al. 2013b; Lemoine & Reinquin 2017; Wu et al. 2017) and we consider also the degree 1 correction from Swenson et al. (2008) which is widely used in the GRACE community (Fig. 1). Lemoine & Reinquin (2017) and Cheng et al. (2013b) estimates are both based on SLR data but they differ in the inversion technique. The Swenson et al. (2008) solution is based on a combination of an ocean model with estimates of degree 2 and higher from GRACE. The Rietbroek et al. (2012) and the Wu et al. (2017) solutions use a combination of observations from GPS and GRACE. All solutions compare well in terms of annual cycle but they show different trends (see Fig. 1). The estimates from Cheng et al. (2013b), Lemoine & Reinquin (2017) and Wu et al. (2017) show a trend in the Z component that is positive (+0.12, +0.18 and +0.24 mm yr⁻¹, respectively) while estimates from Swenson et al. (2008) and Rietbroek et al. (2012) show a negative trend (−0.05 and −0.04 mm yr⁻¹, respectively). Past studies have focused on estimating the annual cycle of the geocentre. But here we are interested in the trends. The trends are more difficult to retrieve than the annual cycle because they are at the limit of what is achievable with some of the approaches developed so far to estimate the geocentre position (Wu et al. 2012). A recent study from Riddell et al. (2017) suggests that the uncertainty in the trend of the geocentre could be up to ±0.21, ±0.28 and ±0.54 mm yr⁻¹ in X, Y and Z, respectively (at 1.65σ level). We adopt here a conservative approach, considering the most recent state-of-the-art published geocentre motions available at the time of writing and use the spread in their trend as an estimate of the uncertainty in the trend. Here, the spread in trend between different geocentre is slightly smaller (±0.30 mm yr⁻¹ on the Z-axis) than in Riddell et al. (2017; see Fig. 1).

The second correction consists in including an estimate of the Earth oblateness which corresponds to the degree 2 order 0 zonal
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Geocentre motion in the (X,Y,Z). Time-series (a) and associated trends over the period January 2005–December 2015 (b) from Swenson *et al.* (2008), Rietbroek *et al.* (2012), Cheng *et al.* (2013b), Lemoine & Reinquin (2017) and Wu *et al.* (2017). A low-pass filter to remove sub-annual frequencies have been applied to the time-series. The trends are represented by error bars considering the 1.65 std from the individual errors.</figcaption></figure>

spherical harmonic coefficient ($C_{2,0}$). In principle, it can be determined from the nodal precession acceleration of GRACE. But GRACE is in a near-polar orbit which makes the satellite-to-satellite tracking technique weakly sensitive to the nodal precession acceleration. As a consequence GRACE estimates of the $C_{2,0}$ coefficient are inaccurate (Wahr *et al.* 2006; Cheng & Ries 2017), and have to be replaced by estimates based on SLR data (Chambers 2006). In this study we consider two recent estimates of $C_{2,0}$ from Cheng *et al.* (2013a) and Lemoine & Reinquin (2017; see Fig. 2). Cheng *et al.* (2013a) estimate is obtained from the analysis of SLR data only, using five geodetic satellites: LAGEOS 1 and 2, Starlette, Stella and Ajisai. Lemoine & Reinquin (2017) use the same SLR raw data (except the Ajisai data) but obtain the $C_{2,0}$ estimate through a joint inversion of the SLR data with GRACE data. Both $C_{2,0}$ estimates are very close (see Fig. 2) in terms of annual to interannual variability and also in terms of trends over the period 2005–2015. This is because the weighting scheme in the joint inversion of Lemoine & Reinquin (2017) makes the calculation of $C_{2,0}$ heavily dependent on SLR data.

## 2.3 Filtering

GRACE solutions show correlated errors in the high degrees and orders coefficients (Swenson & Wahr 2006). This noise, which is more pronounced on the sectorial coefficients than on the tesseral and zonal coefficients, propagates into a north-south striping pattern when the spherical harmonics are converted into equivalent water height grids (EWH). The origin of these stripes lays in the near-polar orbit of GRACE and the associated weak sensitivity of the satellite-to-satellite tracking technique to East-West gravity gradients. Instrument errors, background model inaccuracies and processing errors are minimized in the North-South direction and tend to end up in the East-West gravity gradients where the constraint from observations is weaker. These stripes need to be corrected for, if we

want to retrieve accurate estimates of the surface mass variations at small spatial scales (300–500 km). There are several methods to reduce this noise. The most used are the simple Gaussian smoothing of Jekeli (1981) adapted by Wahr *et al.* (1998), the empirical destriping (Chambers 2006; Swenson & Wahr 2006), the Wiener filtering (Sasgen *et al.* 2006), the empirical orthogonal function filtering (Schrama *et al.* 2007), the two-step destriping method (Swenson & Wahr; Chen *et al.* 2007) and the smoothing with order convolution filters (Kusche *et al.* 2009). Every method carries its limitations. Here we choose to test two methods: the two steps method P3M6 from Chen *et al.* (2007) and the order convolution filters from Kusche *et al.* (2009). The P3M6 method consists in smoothing all the coefficients above the sixth order by computing a polynomial fit of order 3 to the odd and even coefficients of a given order. For the order convolution filter, ICGEM provides eight DDK filters with different degrees of smoothing referred to as DDK1 to DDK8. The DDK1 filter applies the strongest level of smoothing approximately corresponding to a Gaussian smoothing radius of 530 km along the East-West direction. The level of smoothing decreases to the DDK8 filter, which applies the weakest level of smoothing corresponding to a Gaussian radius of 200 km (more information on http://icgem.gfz-potsdam.de/Visualized_DDKkernels.pdf). The large smoothing applied by the DDK1, DDK2 filters and Gaussian filter results in an increased leakage of the ice sheet, and coastal glacier signals into the ocean and a dampened signal on land while the weak smoothing of the DDK7 and DDK8 filters make the destriping inefficient in the ocean (see Fig. S1 in the Supporting Information). For these reasons, these DDK filters are discarded here and we only test four DDK filters, namely DDK3, DDK4, DDK5 and DDK6. The filtering is applied to all GRACE solutions except the GRGS solution which does not need a *posteriori* filtering. Indeed, the GRGS centre uses a truncated single value decomposition scheme (leaving out less than 12 per cent of the total variance) for the inversion instead of a classical Choleski inversion. This allows the resolution of the better-determined linear combinations of
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Time-series of the $C_{20}$ coefficient from Lemoine & Reinquin (2017) and Cheng et al. (2013a). A low-pass filter to remove sub-annual frequencies have been applied to the time-series. Time-series are centred around -4 841 653.22 × 10⁻¹⁰. The trends are computed for the period January 2005–December 2015. $C_{20}$ is adimensional and its trend is expressed in yr⁻¹.</figcaption></figure>

the gravity coefficients and prevents the resolution of the most ill-determined ones, responsible for the stripes in the classical solution.

## 2.4 Leakage correction

At monthly and longer periods, mass variations are in general much larger over land than over the ocean. Close to the coast the large mass signal from land tends to spread out into the nearby ocean because of the limited GRACE resolution. This effect, called ‘leakage’, generates spurious signal on the coastal ocean and an underestimation of the signal on land in the coastal zone. Based on a model approach (Landerer & Swenson 2012), we correct for this land leakage signal by using an independent reference estimate of the ocean mass variations. The reference ocean mass estimate is based on ocean reanalyses that assimilates observations from satellite altimetry and *in situ* temperature and salinity profiles. From the reference ocean mass estimate, we compute for each gridpoint in the coastal ocean, monthly mass anomalies (with respect to the local mean value over 2005–2015) in EWH. We assume these monthly mass anomalies from the reanalyses as references and we consider that any difference between these anomalies and GRACE anomalies is land leakage signal in GRACE solutions. We correct GRACE solutions for this land leakage signal by transferring it from the coastal ocean gridpoints to the closest land gridpoint. We apply this correction to all ocean gridpoints located in the coastal ocean defined as the ocean within 300 km from the coast. Once corrected, we verify that the new mass distribution is compliant with GRACE raw uncertainty (see Text S1 and Fig. S2 in the Supporting Information for the GRACE raw uncertainty and Text S2, Figs S3 and S4 in the Supporting Information for further details on the leakage correction). This leakage correction depends on the reference ocean mass estimate. We test different values for this parameter.

Two different reference ocean mass estimates are considered. They are computed from two ocean reanalysis namely ORA S4 (Balmaseda et al. 2013) and GLORYS2V4 (Garric et al. 2018). These estimates cover the whole ocean including high latitudes

and the coastal zones. They assimilate observations from satellite altimetry and temperature and salinity profiles worldwide. The dynamical core of the ORA S4 and GLORYS2V4 models enable to interpolate the temperature and salinity fields in a physical consistent way in regions where observations are scarce like coastal regions.

We tested different methods to allocate the leakage signal on land (distributing the signal on all land gridpoints with an amplitude inversely proportional to the distance instead of distributing the signal on the closest land point) but it did not yield significant differences in the results (not shown). We tested the boundary conditions of the coastal ocean for the leakage correction and we choose 300 km (see Fig. S5 in the Supporting Information). We also tested solutions for which the destriping was not done before the leakage correction. It yields differences in ocean mass trend (and other components of the global water budget) below 0.03 mm yr⁻¹ (see Fig. S1 in the Supporting Information). Another issue is the ‘ringing’ that the truncation of the spherical harmonics generates. Note that the leakage correction does not prevent for this effect. However, although the ringing presents an important local effect it does not affect the global budgets so we neglect this effect here.

## 2.5 GIA correction

Several solid Earth processes generate significant variations in the Earth gravity field at interannual to decadal timescales and can blur the gravity changes associated with water mass redistributions at the Earth surface. These processes are essentially the solid Earth load tides, the solid Earth pole tides, the earthquakes and the GIA. The solid Earth load tides and pole tides are already corrected for in the solutions provided by the processing centres. Earthquakes generate local mass redistributions on short timescales (over a few days). They can also impact mass redistributions at interannual to multidecadal timescales through long-term post-seismic adjustment. However their effect is local (typically over a few 100 km) and their impact is small at global scale (Reager et al. 2016). In this
study we neglect their effect. In contrast, the GIA signal induces significant trends in GRACE solutions that must be removed. The only way to correct for this signal is to use GIA models. GIA models primarily depend on models of the deglaciation history and the mantle viscosity profile. Here we consider three different models namely AG 2013 (Geruo et al. 2013), ICE 6G_D 2017 (Peltier et al. 2017) and ICE6G_ANU_D 2018 (Purcell et al. 2018). AG 2013 is based on the ICE5G deglaciation model and the mantle viscosity is computed numerically via a finite element method considering a viscoelastic profile that varies laterally. ICE6G_D 2017 is an update of ICE5G_C which was developed simultaneously to the deglaciation model (ICE 6 G; Stuhne & Peltier 2015) and the viscosity mantle model (VM5a). ICE6G_ANU_D is an alternative version of ICE6G_D that is corrected for anomalously large uplift signals over regions where ice has been grounded below sea level at or since the Last Glacial Maximum (Purcell et al. 2016).

Note that regional GIA models exist to correct GRACE in specific regions. These models often use regionally refined sets of observations of the vertical land movement and refined sets of local proxies of the past ice extent in their inversion scheme. Regional GIA models show results that can be different from global models locally. They yield to significantly different ice mass change estimates when they are used to correct local GRACE solutions in particular in Antarctica (Shepherd et al. 2012; Whitehouse et al. 2012; Ivins et al. 2013; Martín-Español et al. 2016) and Greenland (Khan et al. 2016). Here we do not use regional GIA models to correct global GRACE solutions as they give spurious estimates of the GIA signal out of the specific regions for which they have been designed (Whitehouse et al. 2012). They would lead to inconsistent estimates of mass changes at global scale. We only use two regional GIA models over Antarctica (Whitehouse et al. 2012) and over Greenland (Khan et al. 2016) as references to quantify the impact of using regionally refined GIA models rather than global GIA models on GRACE mass trends (see Table S1 in the Supporting Information).

## 2.6 Statistical method to evaluate GRACE uncertainty

Variations in processing and post-processing parameters affect GRACE estimates of the mass redistributions. In order to determine to which extent these variations can bias GRACE estimates, we consider the range of the six processing parameters described in previous sections (i.e. processing centre, geocentre motion, $C_{2,0}$, the filtering parameter, the leakage parameter and the GIA correction) and we compute an ensemble of GRACE post-processed solutions. The range of processing includes five GRACE solutions from five processing centres, five geocentre motion corrections, two $C_{2,0}$ corrections, five filtering techniques, two leakage corrections and three GIA corrections which lead to 1500 possible combinations to generate a post-processed GRACE solution. We thus build an ensemble of 1500 post-processed solutions and use it to estimate the global water budget.

The global water budget is broken down into five components: the ocean mass change, the Antarctica mass change, the Greenland mass change, the arctic islands mass change and remaining glaciers plus terrestrial water storage (TWS) mass change (glaciers & TWS). The ocean component is further broken down into the coastal ocean mass change (within 300 km from the coast, which represents 17 per cent of the ocean surface) and the open ocean mass change because it is essentially the coastal ocean that is affected by variations in

processing techniques (see below). Greenland and Antarctica components include the ice sheet and the peripheral glaciers. Antarctica is defined using the bedmap2 mask (Fretwell et al. 2013). The arctic islands component is defined using the version 4.0 of the Randolph glaciers inventory (RGI) from Pfeffer et al. (2014). It includes islands in the northern and southern arctic Canada (region 3 and 4 in RGI), in Iceland (region 6 in RGI), in Svalbard (region 7 in RGI) and in the Russian arctic (region 9 in RGI). The glaciers & TWS component includes all other regions on Earth (see Fig. S6 in the Supporting Information). All contributions are expressed in mm Sea Level Equivalent (mm SLE is defined as the mass change of a water budget component normalized by the total global ocean area). Note that because water mass is conserved in the Earth system and in GRACE observations (GRACE solutions do not include any degree 0), GRACE estimate of the ocean mass change is equal to the sum of GRACE estimate of the other water budget components.

In this article, we focus on the trends of the global water budget components over the period January 2005–December 2015. January 2005 corresponds to the full deployment of the ARGO profiling floats in the ocean (Roemmich 2009), increasing the performance of ocean reanalysis. However, the ensemble solution is provided from August 2002, which is the first common month from the five chosen processing centres. Note that all the trends in this study are computed by least-squares fitting a first-order polynomial after removing the annual and semi-annual cycles. This method provides an estimate of the formal error associated to the trend estimate. This formal error is very small ($<0.001 \text{ mm yr}^{-1}$) because we use monthly time-series that cover a decade. We neglect this source of uncertainty.

The ensemble is built with the intent to consider the largest number of state-of-the-art processing parameters as possible. Any processing parameter that has been shown to have deficiencies has been discarded (e.g. the filters ddk1 or ddk2). All the other parameters are weighted equally in the ensemble because we could not find any reason in the literature to favour a parameter over another. With this selection, we expect that our ensemble actually gathers most of the current state-of-the-art GRACE solutions. As such, we expect that the variance of our ensemble can provide insights on GRACE uncertainty in global water budget component trends. That said, we acknowledge that our ensemble is not exhaustive. It builds on prior work to develop and refine the best processing parameters and thus it does not include all GRACE estimates. In addition, it does not account for any unknown systematic error that could affect all GRACE solutions. As such, the spread of our ensemble do certainly not represent the true uncertainty in GRACE. But we believe it gives interesting insights on the sources of uncertainty in state-of-the-art GRACE solutions.

To evaluate the role of the variations in a given processing parameter on GRACE uncertainty, we first average the ensemble across the five other processing parameters. Then we compute the variance of the resulting single-parameter ensemble (see Sections 3.1–3.6). To evaluate the total GRACE uncertainty, we compute the variance of the whole 1500-member ensemble (see Section 3.7). This variance can be expressed as the sum of the variance of each single-parameter ensemble plus an interaction term. This method is similar to the ANanalysis Of VAriance procedure – ANOVA e.g. (Fisher 1925) – except that we do not have residual terms here. This method was used by Geoffroy et al. (2012) to assess the contribution of different parameters to the variations of climate sensitivity, and we refer to that study for details on the method. To summarize, the total variance of GRACE reconstruction is decomposed as follows:
Eq. (1) computation of the variance as the sum of the variance of each single-parameter ensemble plus an interaction term

$$ 
\begin{aligned}
\operatorname{var}(\operatorname{Tr}_i(\mathrm{pc}, \mathrm{GM}, C_{2,0}, \mathrm{Fp}, \mathrm{Lp}, \mathrm{GIA})) \\
&= \operatorname{var}\left(\left(\operatorname{Tr}_i(\mathrm{pc})\right)_{\mathrm{GM}, C_{2,0}, \mathrm{Fp}, \mathrm{Lp}, \mathrm{GIA}}\right) \\
&\quad+\operatorname{var}\left(\left(\operatorname{Tr}_i(\mathrm{GM})\right)_{\mathrm{pc}, C_{2,0}, \mathrm{Fp}, \mathrm{Lp}, \mathrm{GIA}}\right) \\
&\quad+\operatorname{var}\left(\left(\operatorname{Tr}_i(C_{2,0})\right)_{\mathrm{pc}, \mathrm{GM}, \mathrm{Fp}, \mathrm{Lp}, \mathrm{GIA}}\right) \\
&\quad+\operatorname{var}\left(\left(\operatorname{Tr}_i(\mathrm{Fp})\right)_{\mathrm{pc}, \mathrm{GM}, C_{2,0}, \mathrm{Fp}, \mathrm{GIA}}\right) \\
&\quad+\operatorname{var}\left(\left(\operatorname{Tr}_i(\mathrm{Lp})\right)_{\mathrm{pc}, \mathrm{GM}, C_{2,0}, \mathrm{Fp}, \mathrm{GIA}}\right) \\
&\quad+\operatorname{var}\left(\left(\operatorname{Tr}_i(\mathrm{GIA})\right)_{\mathrm{pc}, \mathrm{GM}, C_{2,0}, \mathrm{Fp}, \mathrm{Lp}}\right) \\
&\quad+\operatorname{var}(\mathrm{I}(\mathrm{pc}, \mathrm{GM}, C_{2,0}, \mathrm{Fp}, \mathrm{Lp}, \mathrm{GIA})),
\end{aligned}
\quad (1) $$

where the brackets $\langle \rangle_{p_1, \dots, p_5}$ indicate the mean across the five parameters $p_1, \dots, p_5$. $\operatorname{Tr}_i$ indicates GRACE estimate of the trend in the ith global water budget component (ocean mass, Greenland, Antarctica, Artic islands and glaciers & TWS). $\mathrm{pc}, \mathrm{GM}, C_{2,0}, \mathrm{Fp}, \mathrm{Lp}, \mathrm{GIA}$ indicate the six processing parameters (namely the processing centre, the geocentre motion, the $C_{2,0}$, the filtering parameter, the leakage parameter and the GIA parameter) and I indicates the interaction term. Note that this decomposition of the variance of the ensemble is exact. There is no approximation. The interaction term corresponds to the non-additive interaction of processing parameters in the GRACE post-processed solution (i.e. the $\operatorname{Tr}_i$ function) on one or several processing parameters. We evaluate the interaction term in Section 3.7 and discuss the consequences for GRACE uncertainty.

# 3 RESULTS: UNCERTAINTY IN GRACE ESTIMATE OF THE WATER BUDGET COMPONENTS

## 3.1 Uncertainty related to the GRACE processing centre

Fig. 3(a) shows the effect of varying the processing centre parameter on GRACE estimates of the global water budget trends for the period 2005–2015. As described in eq. (1), we compute the trend for each water budget component as the mean of the trends of the sub-ensemble defined by the processing centre parameter. For all water budget components, the trend varies by less than 0.20 mm yr⁻¹ SLE among the different processing centres. The largest differences are obtained for the open ocean mass and glaciers & TWS trends for which largest differences reach 0.20 mm yr⁻¹ SLE. For Antarctica, differences in trends estimates are below 0.09 mm yr⁻¹ SLE and for Greenland and the arctic Islands the differences are smaller (0.03 mm yr⁻¹ SLE).

## 3.2 Uncertainty related to the geocentre motion correction

The geocentre motion parameter has a small impact on the estimate of the Greenland and arctic islands mass trends (smaller than 0.03 mm yr⁻¹ SLE). The Antarctica mass trend differences are below 0.08 mm yr⁻¹ SLE. However ocean and glaciers & TWS mass trends are more sensitive to the geocentre parameter (differences up to 0.42 mm yr⁻¹ SLE). The use of Lemoine & Reinquin (2017) or Wu et al. (2017) geocentre motion parameters yield similar results. They lead to small ocean mass trend (1.49–1.51 mm yr⁻¹ SLE) and small glaciers & TWS mass trends (−0.24

and −0.21 mm yr⁻¹ SLE). The differences between both solutions are below 0.02 mm yr⁻¹ SLE for all components of the global water budget. Using Cheng et al. (2013b) or Rietbroek et al. (2012) yield also similar results but they lead to a different picture (Fig. 3b) compared to Lemoine & Reinquin (2017) or Wu et al. (2017). In the case of Cheng et al. (2013b) and Rietbroek et al. (2012), the ocean mass and glaciers & TWS trends are both significantly larger (1.63–1.68 mm yr⁻¹ SLE) than when we consider the Lemoine & Reinquin (2017) and Wu et al. (2017) parameter. The use of the Swenson et al. (2008) geocentre parameter leads to the largest ocean mass trend (+1.84 mm yr⁻¹) and the largest glaciers & TWS mass trend (+0.14 mm yr⁻¹). These large differences in ocean and glaciers & TWS mass trend among GRACE solutions are essentially caused by differences in trends of the Z component of the geocentre motion parameter (see Barletta et al. 2013 and Table S2 in the Supporting Information for the impact of the x, y and z component of the geocentre motion parameter on the trends of the global water budget).

## 3.3 Uncertainty related to the $C_{2,0}$ correction

The use of different $C_{2,0}$ corrections generate small differences in the water budget component trends. These differences are smaller than 0.03 mm yr⁻¹ for all components (Fig. 3c). The largest effect is on the ocean and Antarctica mass trends.

## 3.4 Uncertainty related to the filtering

The filtering parameter has a small effect. This is for three reasons. First, because we selected here only up-to-date state-of-the-art filters and these filters tend to yield similar results. Second, because we consider separately the leakage effect and the filtering (unlike many previous studies). And third, because we only analyse large areas (see Fig. 4a). We probably would find significant differences among filters at basin or glacier scale. Here we find that the differences in the water budget trends generated by the different filtering are below 0.04 mm yr⁻¹ SLE for all components and all filters.

## 3.5 Uncertainty related to the leakage correction

As explained in Section 2.4, the leakage correction is sensitive to the choice of the reference ocean mass used to correct the coastal ocean. Here the two different reference ocean masses used in the leakage correction generates differences in GRACE estimates of the trends in the water budget components that are below 0.11 mm yr⁻¹ SLE (see Fig. 4b). Because of the localization of the leakage on the coast, only the coastal ocean mass and land components are affected. The effect of this parameter is important on the glaciers & TWS component. We find that it is small on Greenland, arctic Islands and Antarctica mass trends.

## 3.6 Uncertainty related to the GIA correction

There is a good agreement between GRACE estimates of Greenland, the arctic islands and glaciers & TWS mass trends when any of the three GIA model is considered. With ICE-6G_D 2017 and ICE6G_ANU_D 2018 models, similar Antarctica mass trends are obtained. However, the ocean mass trends are different by up to 0.08 mm yr⁻¹ SLE. When AG 2013 is used, the picture is different with a global ocean mass trend and Antarctica mass trend higher (1.72 and 0.73 mm yr⁻¹ SLE) than the two other global GIA models.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Trends in global mass budget (in mm yr⁻¹ SLE over 2005–2015) for all components of the global water budget. Comparison of the mean values over the sub-ensembles according to (a) processing centre, (b) Geocentre motion and (c) C₂₀.</figcaption></figure>

ICE-6G_D 2017 and ICE6G_ANU_D 2018 models lead to similar Antarctica mass trend (Whitehouse et al. 2012; Ivins et al. 2013) while AG 2013 present the highest value of all GIA models available for Antartica (Martín-Español et al. 2016). The reason of the smaller difference comes from the fact that here two of the three global GIA models have been recently updated considering the same BEDMAP2 bathymetry for the Southern Ocean (Peltier et al. 2017; Purcell et al. 2018). For Greenland, the regional model of Khan et al. (2016) yields to a local GIA apparent mass higher than global models, leading to a Greenland SLE mass loss which is 0.06 mm yr⁻¹ SLE smaller than when using global models (see Table S1 in the Supporting Information).

## 3.7 Total uncertainty

Table 1 summarizes the sources of uncertainty in GRACE estimate of the global water budget components associated to each processing parameter. The uncertainty is expressed as 1.65 standard deviation of each single-parameter ensemble corresponding to the 90 per cent

confidence level (CL) when assuming a Gaussian distribution. The uncertainty of the whole ensemble of GRACE solutions (referred to as the ‘total uncertainty’) is computed also as the 1.65 standard deviation of the whole ensemble. We find that the total uncertainty in GRACE estimates of the trends in Greenland and arctic islands mass is below 0.04 mm yr⁻¹. It represents respectively 5 and 13 per cent of the mass trend signal in these regions confirming the accuracy of GRACE estimates of the current Greenland and arctic islands mass loss. For the global ocean, Antarctica and glaciers & TWS mass changes, the total GRACE uncertainty is about one order of magnitude larger than for Greenland and the arctic islands. It reaches respectively 0.27 mm yr⁻¹ SLE, 0.15 mm yr⁻¹ SLE and 0.27 mm yr⁻¹ SLE which represent 16, 23 and 300 per cent of the average mass change signal in these regions confirming earlier studies which point out the large uncertainty in GRACE estimates of the current changes in Glaciers & TWS (Yi et al. 2015; Reager et al. 2016; Rietbroek et al. 2016; Dieng et al. 2017).

Table 1 also shows the interaction term, which is computed as the difference between the total uncertainty and the quadratic sum
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Trends in global mass budget (in mm yr⁻¹ SLE over 2005–2015) for all components of the global water budget. Comparison of the mean values over the sub-ensembles according to (a) filtering parameter, (b) leakage correction and (c) GIA correction.</figcaption></figure>

of the uncertainty associated to each single processing parameter. As explained in Section 2.6 the interaction term arises from the non-additive interaction of parameters in the ensemble of GRACE post-processed solutions. Here it represents less than 1 per cent of the total GRACE uncertainty for all global water budget components. The reason for this low interaction term is that here all estimates of the processing parameters are sampled independently and for this reason they are largely uncorrelated. This approach is partly biased because some estimates of the different processing parameters rely on a common reference frame (e.g. estimates of the geocentre motion, C₂₀ and GIA) that could impact several parameters simultaneously, and thus generate covariance among the processing parameter estimates. So the interaction term here is likely underestimated. However the underestimation is probably small as the biases of the International Terrestrial Reference Frame

are likely small (of the order of a few tenth of mm yr⁻¹; Altamimi et al. 2016).

At interannual timescales, we estimate the uncertainty in GRACE solutions by computing for each year, the standard deviation of the ensemble of all GRACE annual estimates. This uncertainty reaches up to ±4.23 mm SLE for the most uncertain years (at 1.65σ, i.e. 90 per cent CL assuming a Gaussian distribution).

Fig. 5(c) shows the ensemble-mean time-series of the global water budget components. The uncertainty in glaciers & TWS mass variations at interannual timescales reaches up to ±4.28 mm SLE (90 per cent CL). The main contributor to this uncertainty is the uncertainty in geocentre motion (Fig. 6b). For Antarctica, the uncertainty in mass variations at interannual timescales reaches up to ±1.24 mm SLE (90 per cent CL) while for Greenland and arctic islands the uncertainty is smaller (up to ±0.37 mm SLE).
<table><caption>**Table 1.** Uncertainties in trend over 2005–2015 of GRACE estimates of the global water budget components (mm yr⁻¹).</caption><thead><tr><th>(mm yr<sup>-1</sup> SLE)</th><th>Ocean mass</th><th>Greenland</th><th>Antarctica</th><th>Arctic islands</th><th>Glacier &amp; TWS</th></tr></thead><tbody><tr><td>Mean trend</td><td>1.63</td><td>0.80</td><td>0.63</td><td>0.29</td><td>-0.09</td></tr><tr><td>Processing centre</td><td>0.09</td><td>0.02</td><td>0.06</td><td>0.01</td><td>0.12</td></tr><tr><td>Geocentre motion</td><td>0.21</td><td>0.01</td><td>0.05</td><td>0.02</td><td>0.23</td></tr><tr><td>C<sub>2,0</sub></td><td>0.02</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td></tr><tr><td>Filtering</td><td>0.02</td><td>0.03</td><td>0.01</td><td>0.02</td><td>&lt;0.01</td></tr><tr><td>Leakage correction</td><td>0.08</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td><td>0.09</td></tr><tr><td>GIA</td><td>0.12</td><td>&lt;0.01</td><td>0.12</td><td>0.01</td><td>0.03</td></tr><tr><td>Total uncertainty</td><td>0.27</td><td>0.04</td><td>0.15</td><td>0.04</td><td>0.27</td></tr><tr><td>Quadratic sum of individual uncertainties</td><td>0.27</td><td>0.04</td><td>0.15</td><td>0.04</td><td>0.27</td></tr><tr><td>Interaction</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td></tr></tbody></table>

All uncertainties are given at 1.65σ (i.e. 90 per cent CL assuming a Gaussian distribution). Relative uncertainties are provided in per cent and calculated using the mean trend value. Total uncertainty is computed as the 1.65 × RMS of the whole ensemble.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Time-series of the global water budget components after applying a low-pass filter to remove sub-annual frequencies (left) and their trends distributions for 2005–2015 (right). In the right-hand panel, the shaded areas indicate the distribution of trends, while the red lines indicate the distribution that would be obtained under the assumption of a Gaussian distribution. For comparison, estimates have been offset in the left-hand panel.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Impact of the geocentre parameter in the components of the global water budget. The small boxes correspond to the distribution of trends of each sub-ensembles for each of the four geocentre series. A low-pass filter to remove sub-annual frequencies have been applied to the time-series.</figcaption></figure>

<table><caption>**Table 2.** Comparison with previous global water budget estimates in mm yr⁻¹.</caption><thead><tr><th>Period</th><th>Source</th><th>Ocean mass<br>(mm yr<sup>-1</sup> SLE)</th><th>Greenland<br>(mm yr<sup>-1</sup> SLE)</th><th>Antarctica<br>(mm yr<sup>-1</sup> SLEr)</th><th>Rest of the<br>land (mm yr<sup>-1</sup> SLE)</th></tr></thead><tbody><tr><td>Apr 2002–Dec 2014</td><td>Reager et al. 2016<sup>a</sup></td><td>1.58 ± 0.43</td><td>0.77 ± 0.13</td><td>0.49 ± 0.53</td><td>0.32 ± 0.21</td></tr><tr><td>Aug 2002–Dec 2014</td><td>This study</td><td><strong>1.32 ± 0.21</strong></td><td><strong>0.77 ± 0.03</strong></td><td><strong>0.49 ± 0.13</strong></td><td><strong>0.07 ± 0.20</strong></td></tr><tr><td>Apr 2002–Jun 2014</td><td>Rietbroek et al. 2016<sup>a</sup></td><td>1.36 ± 0.69<sup>b</sup></td><td>0.73 ± 0.03</td><td>0.26 ± 0.07</td><td>0.09 ± 0.27<sup>c</sup></td></tr><tr><td>Aug 2002–Jun 2014</td><td>This study</td><td><strong>1.31 ± 0.23</strong></td><td><strong>0.77 ± 0.03</strong></td><td><strong>0.46 ± 0.13</strong></td><td><strong>0.07 ± 0.20</strong></td></tr><tr><td>Jan 2005–Jul 2014</td><td>(Yi et al. 2015)<sup>d</sup></td><td><strong>2.03 ± 0.21</strong></td><td>0.77 ± 0.05</td><td>0.60 ± 0.18</td><td><strong>0.65 ± 0.28</strong><sup>e</sup></td></tr><tr><td>Jan 2005–Jul 2014</td><td>This study</td><td><strong>1.40 ± 0.22</strong></td><td><strong>0.83 ± 0.03</strong></td><td><strong>0.57 ± 0.14</strong></td><td><strong>0.00 ± 0.21</strong></td></tr><tr><td>Jan 2004–Dec 2015</td><td>(Dieng et al. 2017)<sup>a</sup></td><td><strong>2.24 ± 0.16</strong><sup>f</sup></td><td>0.82 ± 0.10</td><td>0.33 ± 0.10</td><td><strong>1.03 ± 0.17</strong></td></tr><tr><td>—</td><td>—</td><td><strong>2.35 ± 0.27</strong><sup>b</sup></td><td>—</td><td>—</td><td>—</td></tr><tr><td>Jan 2004–Dec 2015</td><td>This study</td><td><strong>1.57 ± 0.21</strong></td><td><strong>0.78 ± 0.03</strong></td><td><strong>0.57 ± 0.14</strong></td><td><strong>0.21 ± 0.20</strong></td></tr></tbody></table>

Our GRACE ensemble starts in August 2002 (and not April 2002) because there are some missing months in one of the processing centres. 'Rest of the land' corresponds to the sum of arctic islands and glaciers & TWS from our study. Uncertainty is expressed in 1.65σ (90 per cent of CL). Red values do not overlap within 1.65σ (90 per cent of CL).

# 4 DISCUSSION

## 4.1 Impact of the geocentre motion in the global water budget

As shown by previous studies (Barletta et al. 2013; Reager et al. 2016) most of the GRACE uncertainty related to the geocentre is due to the uncertainty in the geocentre velocity along the Z-axis. We

tested this hypothesis by detrending the z component of geocentre solutions and recomputing GRACE solutions. This experience resulted in a significant reduction of the spread in GRACE estimates of the ocean mass trend confirming previous studies (the spread among the trends dropped from ±0.35 to ±0.17 mm yr⁻¹ SLE, see Fig. S7 in the Supporting Information). Here the five different estimates of the geocentre motion (Swenson et al. 2008; Rietbroek

<footnote><sup>a</sup>Values in the article are expressed in 1σ</footnote>

<footnote><sup>b</sup>Ocean mass computed as altimetry (total) minus steric from the reference</footnote>

<footnote><sup>c</sup>Rest of land computed as sum of glaciers and hydro from the reference</footnote>

<footnote><sup>d</sup>Values in the article are expressed in 1.96σ (95 per cent of CL)</footnote>

<footnote><sup>e</sup>Rest of land computed as glacier, ice caps and land water from the reference</footnote>

<footnote><sup>f</sup>Ocean mass is based on Chambers & Bonin (2012)</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison of our ensemble with previous estimates of some of the water budget components trends. Vertical lines represent ±1.65 standard deviation around the mean values and the grey shaded areas the distribution of trends from our ensemble. The estimates from Table 2 are included to complete the comparison.</figcaption></figure>

et al. 2012; Cheng et al. 2013b; Lemoine & Reinquin 2017; Wu et al. 2017) show a geocentre velocity along Z ranging from -0.05 and +0.24 mm yr⁻¹, and a dispersion around the ensemble mean in general agreement with recent studies (Riddell et al. 2017 – shows an uncertainty in the trend of the z component of the geocentre of ±0.54 mm yr⁻¹). However our range may underestimate the real range as it is based on only five geocentre solutions and it does not take into account any potential sources of systematic bias.

## 4.2 Comparison with previous estimates

Table 2 and Fig. 7 show the comparison of the trends estimates of the global water budget components from our ensemble with previous

estimates over similar time period (Shepherd et al. 2012; Gardner et al. 2013; McMillan et al. 2014; Mémin et al. 2014; Purkey et al. 2014; Schrama et al. 2014; Velicogna et al. 2014; Nilsson et al. 2015; Yi et al. 2015; Reager et al. 2016; Rietbroek et al. 2016; Dieng et al. 2017). Overall, the general agreement (within uncertainties at 1.65σ) of our solution with previous published estimates of the trends in the water budget components gives confidence in the choice of post-processing parameters. There are only two cases for which our estimate significantly differs from previous estimates: the ocean mass and glaciers & TWS trends estimates from Yi et al. (2015) and from Dieng et al. (2017). Dieng et al. (2017), is based on older data sets which tend to be higher than what is estimated here.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison between mascon solutions and a subset of our ensemble for the (a) global ocean mass evolution, (b) global ocean mass trends for 2005–2015 and (c) trends in global mass budget (in mm yr⁻¹ SLE over 2005–2015) for all components of the global water budget. A low-pass filter to remove sub-annual frequencies has been applied to the time-series.</figcaption></figure>

## 4.3 Comparison with mascons solutions

In this section, we compare our ensemble with the mascon solutions from JPL (RL05M_1.MSCNv02CRIv02 solution with Coastal Regional Improvement filter applied; Watkins et al. 2015), CSR (the CSR RL05_Mascons_v01; Save et al. 2016) and GFSC (the Global Solution v2.3; Luthcke et al. 2013). For a rigorous comparison, we select the subset of our ensemble that is consistent with the mascon solutions in terms of post-processing. The subset includes all solutions that use the Swenson et al. (2008) geocentre motion correction, the Cheng et al. (2013a) C₂₀ correction and the Geruo et al. (2013) GIA model. JPL and CSR solutions are released with the ocean model and the atmosphere model restored over the ocean. To ensure that these solutions are consistent with the spherical harmonic solution, we restore the mean atmosphere load over the ocean (Landerer, personal communication 2017). We use the land/ocean mask of each mascon centre to compute the trends in the water budget components.

Regarding the trend in global ocean mass (Fig. 8b), there is a good agreement between our subset and the JPL mascon estimate while the CSR mascon estimate is significantly smaller (by 0.09 mm yr⁻¹, 90 per cent CL) and the GSFC estimate is significantly larger (by 0.14 mm yr⁻¹ at the 90 per cent CL). At interannual timescales (Fig. 8a), our sub-ensemble is in general agreement with mascons solutions except during La Nina 2011 where it differs significantly from all mascons solutions. Regarding the glaciers & TWS trend

component our sub-ensemble differs with the three mascons solutions but it is not an outlier (0.13 mm yr⁻¹ SLE versus 0.35, 0.23 and 0.02 mm yr⁻¹ SLE). Regarding the Antarctica trend our sub-ensemble presents the highest mass loss (0.68 mm yr⁻¹ SLE versus 0.64, 0.64 and 0.62 mm yr⁻¹ SLE). In terms of the other global water budget components, there is a good agreement (±0.05 mm yr⁻¹ SLE). Note that CSR solution do not close the global water budget because of some interpolation issues between the hexagonal computation grid and the final product that is in a rectangular grid (Save, personal communication 2017). The JPL solution closes the mass budget with a precision of 0.01 mm yr⁻¹ SLE. The small non-closure of 0.01 mm yr⁻¹ probably comes from a residual trend of 0.01 mm yr⁻¹ in the atmospheric background model (over land) which has not been restored here.

## 4.4 Implications for the global water budget, the sea level budget and the Earth energy imbalance

Based on our ensemble, we find a total uncertainty in global ocean mass from GRACE of ±0.27 mm yr⁻¹ over 2005–2015 (SLE at the 90 per cent CL). This uncertainty is comparable to the uncertainty in thermal expansion [±0.15 mm yr⁻¹, Dieng et al. (2017) and ±0.27 mm yr⁻¹, Desbruyères et al. (2016)] and the uncertainty in global mean sea level from satellite altimetry over the same period (which is of ±0.33 mm yr⁻¹, Ablain et al. 2015).

When GRACE observations of the ocean mass changes are combined with satellite altimetry, it gives an alternative estimate of
the thermal expansion with an uncertainty of ±0.30 mm yr⁻¹ over 2005–2015 (at the 90 per cent CL assuming that the uncertainty from GRACE and satellite altimetry are independent). Assuming a global expansion efficiency of heat of 0.12 ± 0.01 m YJ⁻¹ (Levitus et al. 2012) we find that this uncertainty implies an uncertainty in the net warming of the ocean and the Earth energy budget of ±0.25 W m⁻² over 2005–2015 (at 90 per cent CL assuming Gaussian distribution) when inferred using the sea level budget approach.

## 5 CONCLUSION

We compared trends in GRACE estimates of the global water budget components over 2005–2015 for different set of processing parameters. We considered variations in six different processing parameters namely the processing centre of the raw GRACE solutions, the geocentre motion and C₂₀ corrections, the filtering, the leakage correction and the GIA correction. With all possible combinations of the different processing parameters we computed an ensemble of 1500 post-processed GRACE solutions from which we evaluated the trends in global water budget components and their associated uncertainty. As an estimate of the uncertainty we consider the spread around the ensemble mean. This approach has limitations and may not reveal the true uncertainty of the GRACE estimates but it gives interesting insights on the sources of GRACE errors as it enables to explore the dependency of GRACE uncertainty to the different processing parameters.

The analysis of the uncertainties associated with GRACE solutions shows that they are dominated by the uncertainty in the geocentre motion correction and the GIA correction. The uncertainty in the geocentre motion generates an uncertainty in GRACE estimates of ±0.21 mm yr⁻¹ on the global ocean mass, ±0.01 mm yr⁻¹ SLE on Greenland mass loss, ±0.05 mm yr⁻¹ SLE on Antarctica mass loss, ±0.02 mm yr⁻¹ SLE on arctic islands mass loss and ±0.23 mm yr⁻¹ SLE on glaciers & TWS changes over 2005–2015 (at 1.65σ, i.e. 90 per cent CL).

The uncertainty in the GIA correction generates an uncertainty in GRACE estimates of ±0.12 mm yr⁻¹ on the global ocean mass, ±0.01 mm yr⁻¹ SLE on the Greenland mass loss, ±0.12 mm yr⁻¹ SLE on the Antarctica mass loss, ±0.01 mm yr⁻¹ SLE on arctic islands mass loss and of ±0.03 mm yr⁻¹ SLE on the glaciers & TWS changes over 2005–2015 (at 1.65σ, i.e. 90 per cent CL). This uncertainty in GIA is likely underestimated because it does not take into account the total uncertainty in the ice history (we only considered here GIA solutions based on the ICE5G and ICE6G models) and two of the three models have updated and converged in the estimate of the Antarctica GIA mass apparent trend with the regional GIA models. This is evidenced in Greenland where the use of a regional GIA model yield a mass trend that is significantly smaller (by 0.06 mm yr⁻¹) than our ensemble estimate.

There are other sources of uncertainty in GRACE solutions that have not been taken into account here like the effect of earthquakes or the uncertainties in the Earth reference frame that are not related to the geocentre motion. But they are likely small compared to the GIA and the geocentre motion effects (Reager et al. 2016).

A more general caveat concerning our uncertainty estimate is that our calculations are based on an ensemble of products which is limited (because of a limited available sampling of the parameters), and which does not explore all sources of uncertainties (because some uncertainty sources such as earthquakes, are not represented). These

two characteristics make us think that our uncertainty estimates are likely underestimated.

Compared with mascons solutions, we find a good agreement in terms of interannual variability and trends with JPL solution. However, with the CSR and GSFC solutions, we find significant discrepancies in terms of trends.

Compared to the other GRACE-based studies (listed in Table 2 and Fig. 7), our estimates of the uncertainties are in general agreement for all global water budget components

When the uncertainty on the geocentre motion is included, the total uncertainty of GRACE on the global water and sea level budget is ±0.27 mm yr⁻¹ SLE (at 1.65σ, i.e 90 per cent CL; Table 1). This uncertainty is similar to the uncertainty in sea level estimate. It fixes a limit on the constraints that the budget approach provides on missing or poorly known contributions to sea level rise (such as the deep ocean contribution; see Dieng et al. 2015b) or on related essential climate variable such as the total ocean heat content from which the Earth’s energy imbalance can be deduced (Llovel et al. 2014; Dieng et al. 2015b). These results call for more research to refine our estimates of the geocentre motion and the geoid response to GIA.

## ACKNOWLEDGEMENTS

The Grace Level 3 data were obtained from ICGEM (www.icgemanz-potsdam.de/ICGEM/) and GRGS (grgs.obs-mip.fr/grace). Ocean reanalysis data from www.ecmwf.int/en/research/climate-reanalyses/ocean-reanalysis and www.marine.copernicus.eu. We acknowledge support from CNES.

The ensemble of post-processed GRACE solutions created for this study is available online at ftp.legos OBSMIP.FR/PUB/SOA/gravimetrie/grace_legos

We want to thank Dr A, Dr Chen, Dr Cheng, Dr Li, Dr Rietbroek, Dr Save and Dr Wu for providing data. We also want to thank Dr Chambers, Dr Horwath, Dr Kusche, Dr Landerer and an anonymous reviewer for their constructive remarks, corrections and suggestions that improved this article significantly.

*Author contributions:* BM and AB designed the study. AB made the data analysis. AB and BM lead the writing of the paper. Others co-authors discussed the results and contributed to the writing. CNES supported the work of BM, AB and EB.

## REFERENCES

Ablain, M. et al., 2015. Improved sea level record over the satellite altimetry era (1993–2010) from the Climate Change Initiative project, *Ocean Sci.*, **11**(1), 67–82.

Altamimi, Z., Rebischung, P., Métivier, L. & Collilieux, X., 2016. ITRF2014: a new release of the International Terrestrial Reference Frame modeling nonlinear station motions, *J. geophys. Res.*, **121**(8).

Balmaseda, M.A., Mogensen, K. & Weaver, A.T., 2013. Evaluation of the ECMWF ocean reanalysis system ORAS4, *Q. J. R. Meteorol. Soc.*, **139**(674), 1132–1161.

Barletta, V.R., Sørensen, L.S. & Forsberg, R., 2013. Scatter of mass changes estimates at basin scale for Greenland and Antarctica, *Cryosphere*, **7**(5), 1411–1432.

Chambers, D.P., 2006. Evaluation of new GRACE time-variable gravity data over the ocean, *Geophys. Res. Lett.*, **33**(17).

Chambers, D.P., 2009. Calculating trends from GRACE in the presence of large changes in continental ice storage and ocean mass, *Geophys. J. Int.*, **176**(2), 415–419.
Chambers, D.P. & Bonin, J.A., 2012. Evaluation of Release-05 GRACE time-variable gravity coefficients over the ocean, *Ocean Sci.*, **8**(5), 859–868.

Cheng, M. & Ries, J., 2017. The unexpected signal in GRACE estimates of C20, *J. Geod.*, **91**(8), 897–914.

Cheng, M., Tapley, B.D. & Ries, J.C., 2013a. Deceleration in the Earth's oblateness, *J. geophys. Res.*, **118**(2), 740–747.

Cheng, M.K., Ries, J.C. & Tapley, B.D., 2013b. Geocenter variations from analysis of SLR data. In *Reference Frames for Applications in Geosciences*, pp. 19–25, eds Altamimi, Z. & Collilieux, X., Springer.

Chen, J., Famiglietti, J.S., Scanlon, B.R. & Rodell, M., 2015a. Groundwater storage changes: present status from GRACE observations, *Surv. Geophys.*, **37**(2), 397–417.

Chen, J.L., Wilson, C.R., Li, J. & Zhang, Z., 2015b. Reducing leakage error in GRACE-observed long-term ice mass change: a case study in West Antarctica, *J. Geod.*, **89**(9), 925–940.

Chen, J.L., Wilson, C.R., Tapley, B.D. & Grand, S., 2007. GRACE detects coseismic and postseismic deformation from the Sumatra-Andaman earthquake, *Geophys. Res. Lett.*, **34**(13), L13302.

Church, J.A. et al., 2013. Sea level change. In *Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*, eds Stocker, T.F. et al. Cambridge University Press.

Cretaux, J.F., Soudarin, L., Davidson, F.J.M., Gennero, M.C., Berge-Nguyen, M. & Cazenave, A., 2002. Seasonal and interannual geocenter motion from SLR and DORIS measurements: comparison with surface loading data, *J. geophys. Res.*, **107**(B12), ETG 16–1–ETG 16–9.

Desbruyères, D.G., Purkey, S.G., McDonagh, E.L., Johnson, G.C. & King, B.A., 2016. Deep and abyssal ocean warming from 35 years of repeat hydrography, *Geophys. Res. Lett.*, **43**(19).

Dieng, H.B., Cazenave, A., Meyssignac, B. & Ablain, M., 2017. New estimate of the current rate of sea level rise from a sea level budget approach, *Geophys. Res. Lett.*, **44**(8).

Dieng, H.B., Champollion, N., Cazenave, A., Wada, Y., Schrama, E. & Meyssignac, B., 2015a. Total land water storage change over 2003–2013 estimated from a global mass budget approach, *Environ. Res. Lett.*, **10**(12).

Dieng, H.B., Palanisamy, H., Cazenave, A., Meyssignac, B. & von Schuckmann, K., 2015b. The sea level budget since 2003: inference on the deep ocean heat content, *Surv. Geophys.*, **36**(2), 209–229.

Fisher, R.A., 1925. *Statistical Methods for Research Workers*, Oliver & Boyd.

Flechtner, F., Dobslaw, H. & Fagiolini, E., 2015. *AOD1B Product Description Document for Product Release 05 Rev4.4 GRACE 327–750*, Geo Forschungszentrum Potsdam, Potsdam, Germany.

Fretwell, P. et al., 2013. Bedmap2: improved ice bed, surface and thickness datasets for Antarctica, *Cryosphere*, **7**(1), 375–393.

Fritsche, M., Dietrich, R., Rülke, A., Rothacher, M. & Steigenberger, P., 2010. Low-degree earth deformation from reprocessed GPS observations, *GPS Solutions*, **14**(2), 165–175.

Gardner, A.S. et al., 2013. A reconciled estimate of glacier contributions to sea level rise: 2003 to, 2009, *Science*, **340**(6134), 852–857.

Garric, G. et al., 2018. Performance and quality assessment of the global ocean eddy-permitting physical reanalysis GLORYS2V4, in *Proceedings of the Eight EuroGOOS International Conference*, Paper Presented at the Operational Oceanography Serving Sustainable Marine Development, EuroGOOS, Brussels, Belgium. 2018: Bergen, Norway.

Geoffroy, O., Saint-Martin, D. & Ribes, A., 2012. Quantifying the sources of spread in climate change experiments, *Geophys. Res. Lett.*, **39**(24).

Geruo, A., Wahr, J. & Zhong, S., 2013. Computations of the viscoelastic response of a 3-D compressible Earth to surface loading: an application to Glacial Isostatic Adjustment in Antarctica and Canada, *Geophys. J. Int.*, **192**(2), 557–572.

Ivins, E.R., James, T.S., Wahr, J.O., Schrama, E.J., Landerer, F.W. & Simon, K.M., 2013. Antarctic contribution to sea level rise observed by GRACE with improved GIA correction, *J. geophys. Res.*, **118**(6), 3126–3141.

Jacob, T., Wahr, J., Pfeffer, W.T. & Swenson, S., 2012. Recent contributions of glaciers and ice caps to sea level rise, *Nature*, **482**(7386), 514–518.

Jekeli, C., 1981. Modifying Stokes' function to reduce the error of geoid undulation computations, *J. Geophys. Res.*, **86**(B8), 6985–6990, doi:10.1029/JB086iB08p06985.

Khan, S.A. et al., 2016. Geodetic measurements reveal similarities between post-Last Glacial Maximum and present-day mass loss from the Greenland ice sheet, *Sci. Adv.*, **2**(9).

Klinger, B., Mayer-Gürr, T., Behzadpour, S., Ellmer, M., Kvas, A. & Zehentner, N., 2016. The new ITSG-Grace2016 release, in *EGU 2016 General Assembly*, 2016, Vienna, Austria.

Kusche, J., Schmidt, R., Petrovic, S. & Rietbroek, R., 2009. Decorrelated GRACE time-variable gravity solutions by GFZ, and their validation using a hydrological model, *J. Geod.*, **83**(10), 903–913.

Landerer, F.W. & Swenson, S.C., 2012. Accuracy of scaled GRACE terrestrial water storage estimates, *Water Resour. Res.*, **48**(4), W04531.

Lemoine, J.-M., Bourgogne, S., Biancale, R., Bruinsma, S. & Gégout, P., 2016. CNES/GRGS solutions Focus on the inversion process, in *Paper presented at the GRACE Science Team Meeting*, A1-02. Berlin, Germany.

Lemoine, J.-M. & Reinquin, F., 2017. Processing of SLR observations at CNES, *Newsletter EGSIEM*, October, Page 3.

Leuliette, E.W. & Willis, J.K., 2011. Balancing the sea level budget, *Oceanography*, **24**(2), 122–129.

Levitus, S. et al., 2012. World ocean heat content and thermosteric sea level change (0-2000 m), 1955-2010: World ocean heat content, *Geophys. Res. Lett.*, **39**(10).

Llovel, W., Willis, J.K., Landerer, F.W. & Fukumori, I., 2014. Deep-ocean contribution to sea level and energy budget not detectable over the past decade, *Nature Clim. Change*, **4**(11), 1031–1035.

Llovel, W. et al., 2011. Terrestrial waters and sea level variations on inter-annual time scale, *Glob. Planet. Change*, **75**(1–2), 76–82.

Luthcke, S.B., Sabaka, T.J., Loomis, B.D., Arendt, A.A., McCarthy, J.J. & Camp, J., 2013. Antarctica, Greenland and Gulf of Alaska land-ice evolution from an iterated GRACE global mascon solution, *J. Glaciol.*, **59**(216), 613–631.

Martín-Español, A., King, M.A., Zammit-Mangion, A., Andrews, S.B., Moore, P. & Bamber, J.L., 2016. An assessment of forward and inverse GIA solutions for Antarctica, *J. geophys. Res.*, **121**(9).

McMillan, M., Shepherd, A., Sundal, A., Briggs, K., Muir, A., Ridout, A., Hogg, A. & Wingham, D., 2014. Increased ice losses from Antarctica detected by CryoSat-2, *Geophys. Res. Lett.*, **41**(11), 3899–3905.

Mémin, A., Flament, T., Rémy, F. & Llubes, M., 2014. Snow- and ice-height change in Antarctica from satellite gravimetry and altimetry data, *Earth planet. Sci. Lett.*, **404**, 344–353.

Nilsson, J., Sandberg Sørensen, L., Barletta, V.R. & Forsberg, R., 2015. Mass changes in Arctic ice caps and glaciers: implications of regionalizing elevation changes, *Cryosphere*, **9**(1), 139–150.

Peltier, W.R., Argus, D.F. & Drummond, R., 2017. Comment on “An Assessment of the ICE-6G-C (VM5a) Glacial Isostatic Adjustment Model” by Purcell et al, *J. geophys. Res.*, **123**(2), 2019–2028.

Pfeffer, W.T. et al., 2014. The Randolph Glacier Inventory: a globally complete inventory of glaciers, *J. Glaciol.*, **60**(221), 537–552.

Purcell, A., Tregoning, P. & Dehecq, A., 2016. An assessment of the ICE6G-C(VM5a) glacial isostatic adjustment model, *J. geophys. Res.*, **121**(5).

Purcell, A., Tregoning, P. & Dehecq, A., 2018. Reply to Comment by W. R. Peltier, D. F. Argus, and R. Drummond on “An Assessment of the ICE6G-C (VM5a) Glacial Isostatic Adjustment Model.”, *J. geophys. Res.*, **123**(2), 2029–2032.

Purkey, S.G., Johnson, G.C. & Chambers, D.P., 2014. Relative contributions of ocean mass and deep steric changes to sea level rise between 1993 and 2013, *J. geophys. Res.*, **119**(11), 7509–7522.

Quinn, K.J. & Ponte, R.M., 2010. Uncertainty in ocean mass trends from GRACE, *Geophys. J. Int.*, **181**(2), 762–768.

Reager, J.T., Gardner, A.S., Famiglietti, J.S., Wiese, D.N., Eicker, A. & Lo, M.-H., 2016. A decade of sea level rise slowed by climate-driven hydrology, *Science*, **351**(6274), 699–703.

Riddell, A.R., King, M.A., Watson, C.S., Sun, Y., Riva, R.E.M. & Rietbroek, R., 2017. Uncertainty in geocenter estimates in the context of ITRF2014, *J. geophys. Res.*, **122**(5), 4020–4032..
Rietbroek, R., Brunnabend, S.-E., Kusche, J., Schröter, J. & Dahle, C., 2016. Revisiting the contemporary sea-level budget on global and regional scales, *Proc. Natl. Acad. Sci. USA*, **113**, 1504–1509.

Rietbroek, R., Fritsche, M., Brunnabend, S.-E., Daras, I., Kusche, J., Schröter, J., Flechtner, F. & Dietrich, R., 2012. Global surface mass from a new combination of GRACE, modelled OBP and reprocessed GPS data, *J. Geodyn.*, **59–60**, 64–71.

Roemmich, D., 2009. Argo the challenge of continuing 10 years of progress, *Oceanography*, **22**(3), 46–55.

Sasgen, I., Martinec, Z. & Fleming, K., 2006. Wiener optimal filtering of GRACE data, *Stud. Geophys. Geod.*, **50**(4), 499–508.

Save, H., Bettadpur, S. & Tapley, B.D., 2016. High-resolution CSR grace RL05 mascons: high-resolution CSR grace RL05 mascons, *J. geophys. Res.*, **121**(10), 7547–7569.

Schrama, E.J.O., Wouters, B. & Lavallée, D.A., 2007. Signal and noise in Gravity Recovery and Climate Experiment (GRACE) observed surface mass variations, *J. geophys. Res.*, **112**(B8).

Schrama, E.J.O., Wouters, B. & Rietbroek, R., 2014. A mascon approach to assess ice sheet and glacier mass balances and their uncertainties from GRACE data, *J. geophys. Res.*, **119**(7), 6048–6066.

Shepherd, A. et al., 2012. A reconciled estimate of ice-sheet mass balance, *Science*, **338**(6111), 1183–1189.S

Stuhne, G.R. & Peltier, W.R., 2015. Reconciling the ICE-6G_C reconstruction of glacial chronology with ice sheet dynamics: The cases of Greenland and Antarctica, *J. geophys. Res.*, **120**(9), 1841–1865.

Swenson, S., Chambers, D. & Wahr, J., 2008. Estimating geocenter variations from a combination of GRACE and ocean model output, *J. geophys. Res.*, **113**(B8).

Swenson, S. & Wahr, J., 2006. Post-processing removal of correlated errors in GRACE data, *Geophys. Res. Lett.*, **33**(8).

Tapley, B.D., Bettadpur, S., Ries, J.C., Thompson, P.F. & Watkins, M.M., 2004. GRACE Measurements of mass variability in the Earth system, *Science*, **305**(5683), 503–505.

Velicogna, I., Sutterley, T.C. & van den Broeke, M.R., 2014. Regional acceleration in ice mass loss from Greenland and Antarctica using GRACE time-variable gravity data, *Geophys. Res. Lett.*, **41**(22), 8130–8137.

Velicogna, I. & Wahr, J., 2013. Time-variable gravity observations of ice sheet mass balance: Precision and limitations of the GRACE satellite data, *Geophys. Res. Lett.*, **40**(12), 3055–3063.

Wahr, J., Molenaar, M. & Bryan, F., 1998. Time variability of the Earth's gravity field: hydrological and oceanic effects and their possible detection using Grace, *J. geophys. Res.*, **103**(B12), 30 205–30 229.

Wahr, J., Nerem, R.S. & Bettadpur, S.V., 2015. The pole tide and its effect on GRACE time-variable gravity measurements: implications for estimates of surface mass variations, *J. geophys. Res.*, **120**(6).

Wahr, J., Swenson, S. & Velicogna, I., 2006. Accuracy of GRACE mass estimates, *Geophys. Res. Lett.*, **33**(6), L06401.

Watkins, M.M. & Eanes, R.J., 1997. Observations of tidally coherent diurnal and semidiurnal variations in the geocenter, *Geophys. Res. Lett.*, **24**(17), 2231–2234.

Watkins, M.M., Wiese, D.N., Yuan, D.-N., Boening, C. & Landerer, F.W., 2015. Improved methods for observing Earth's time variable mass distribution with GRACE using spherical cap mascons, *J. geophys. Res.*, **120**(4).

Whitehouse, P.L., Bentley, M.J., Milne, G.A., King, M.A. & Thomas, I.D., 2012. A new glacial isostatic adjustment model for Antarctica: calibrated and tested using observations of relative sea-level change and present-day uplift rates, *Geophys. J. Int.*, **190**(3), 1464–1482.

Wouters, B., Bonin, J.A., Chambers, D.P., Riva, R.E.M., Sasgen, I. & Wahr, J., 2014. GRACE, time-varying gravity, Earth system dynamics and climate change, *Rep. Prog. Phys.*, **77**(11).

Wu, X., 2010. Simultaneous Estimation of Global Present-Day Water Transport and Glacial Isostatic Adjustment, *Nature Geosci.*, **3**, 642–646.

Wu, X., Kusche, J. & Landerer, F.W., 2017. A new unified approach to determine geocentre motion using space geodetic and GRACE gravity data, *Geophys. J. Int.*, **209**(3), 1398–1402.

Wu, X., Ray, J. & van Dam, T., 2012. Geocenter motion and its geodetic and geophysical implications, *J. Geodyn.*, **58**, 44–61.

Yi, S., Sun, W., Heki, K. & Qian, A., 2015. An increase in the rate of global mean sea level rise since 2010: increase in sea level rise rate, *Geophys. Res. Lett.*, **42**, 3998–4006.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

**Figure S1.** Filtering parameter comparison.

**Figure S2.** Formal error in raw solutions, *a priori* error and uncertainty of GRACE raw solutions.

**Figure S3.** Evaluation of the leakage correction on a synthetic field.

**Figure S4.** Evaluation of the leakage correction on Jan 2005 for 241st case of the ensemble.

**Figure S5.** Analysis of the boundary of the coastal ocean.

**Figure S6.** Earth mask.

**Figure S7.** Effect of the Z geocentre motion trend over the Global ocean mass.

**Table S1.** GIA apparent mass change expressed in mm SLE yr⁻¹.

**Table S2.** Unitary effect of each axis of the geocentre motion on the components of the global water budget.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.