Scatter of mass changes estimates at basin scale for Greenland and Antarctica

. During the last decade, the GRACE mission has provided valuable data for determining the mass changes of the Greenland and Antarctic ice sheets. Yet, discrepancies still exist in the published mass balance results, and comprehensive analyses on the sources of errors and discrepancies are lacking. Here, we present monthly mass changes together with trends derived from GRACE data at basin scale for both the Greenland and Antarctic ice sheets, and we as-sess the variability and errors for each of the possible sources of discrepancies, and we do this in an unprecedented systematic way, taking into account mass inference methods, data sets and background models. We ﬁnd a very good agreement between the monthly mass change results derived from two independent methods, which represents a cross validation. For the monthly solutions, we ﬁnd that most of the scatter is caused by the use of the two different data sets rather than the two different methods applied. Besides the well-known GIA trend uncertainty, we ﬁnd that the geocenter motion and the recent de-aliasing corrections signiﬁcantly impact the trends, with contributions of + 13.2 Gt yr − 1 and − 20 Gt yr − 1 , respectively, for Antarctica, which is more affected by these than Greenland. We show differences between

−20 Gt yr −1 , respectively, for Antarctica, which is more affected by these than Greenland. We show differences between the use of release RL04 and the new RL05 and confirm a lower noise content in the new release. The overall scatter of the solutions well exceeds the uncertainties propagated from the data errors and the leakage (as done in the past); hence we calculate new sound total errors for the monthly solutions and the trends.
We find that the scatter in the monthly solutions caused by applying different estimates of geocenter motion time series (degree-1 corrections) is significant -contributing with up to 40% of the total error.

Introduction
Determining reliably the mass balance (MB) of the large continental ice sheets is a major challenge. The results are of great societal importance, especially in terms of global sea level rise. The IPCC report of 2007 (Intergovernmental Panel on Climate Change, 2007) stated that the contribution from the ice sheets was insufficiently constrained, and since then a large effort has been made to improve the mass balance estimates, using different methods and data.
In 2005, Velicogna and Wahr (2005) showed for the first time the possibility of using data from the Gravity and Climate Recovery Experiment (GRACE) mission to determine the mass balance of the Greenland ice sheet. Since then many mass balance estimates of both Greenland and Antarctica have been published, both on ice sheet scale (Chen et al., 2006b;Ramillien et al., 2006;Forsberg and Reeh, 2007;Barletta et al., 2008;Velicogna, 2009) and drainage basin scale (Luthcke et al., 2006;Wouters et al., 2008;Schrama and Wouters, 2011;Sasgen et al., 2012).
The gravity changes observed by GRACE provide a measurement of the mass variations without the need to rely on volume-to-mass conversion assumptions.
However, when using GRACE data to estimate continental mass balance, many corrections are applied either for Published by Copernicus Publications on behalf of the European Geosciences Union.
hydrological purpose or for ice mass changes. The GRACE project centers remove some contributions like those due to (solid Earth and ocean) tides, the atmosphere and the ocean. Other contributions have to be removed in a second stage by the user in order to focus just on isolating the ice mass changes in selected regions. Each of these corrections is a potential source of scatter in the mass balance computations. Throughout the paper we will refer to the distribution of the results as a consequence of the different choice of methods, data sets, and model parameterization with the terms scatter and variability interchangeably. The word variability is used in GRACE literature mostly to indicate the time variation. However we use here the word in its broader meaning, i.e. indicating the general property of a mathematical quantity (e.g. a function) to vary with respect to the variation of the parameters it depends on (e.g. variables). We show, through error and uncertainties evaluations, the differences among many possible alternative solutions of the monthly mass changes, i.e. the variability of the solutions of the monthly mass changes with respect to some of the input used in the process of deriving it, and the process itself.
Glacial isostatic adjustment (GIA) has been recognized as the major cause of uncertainty in ice mass balance (Velicogna and Wahr, 2006;Barletta et al., 2008), especially in Antarctica. GIA is the solid Earth phenomenon responsible for the mantle flow from the equatorial region towards the Pleistocene deglaciated areas. This mass flow within the Earth produces variations in topography and in gravity, and the latter is detected by GRACE mainly as a positive mass trend that, in many areas, cannot easily be distinguished from surface mass accumulation.
However, for mass balance estimates, other sources of variability exist, whose importance has sometimes been underestimated. One source of uncertainty is related to geocenter motion (Chambers, 2006) and the different ways to infer it (Wu et al., 2012).
Large mass variations, like the grounded ice melting in Greenland and Antarctica, generate not only gravity changes but also Earth surface displacements. These translate into a displacement of the center of mass (CM) with respect to the geometric center of the Earth (center of figure CF): this movement is known as the geocenter motion.
Therefore, mass balance estimates should carefully take into account the geocenter motion. However, since the GRACE satellites move together with the center of mass of the Earth, in principle they cannot detect the geocenter motion, so this effect has to be recovered by other means. The choice for the correction for geocenter motion is still an open issue: there are different ways to calculate it and therefore more than one possible correction exist, as thoroughly discussed by Wu et al. (2012), and this can be a source of variability in GRACE-derived mass balance estimates.
Another source of variability in mass balance computation arises from different strategies for processing the raw data by the GRACE project centers; thus discrepancies between different solutions of GRACE data are known (Steffen et al., 2010;Barletta et al., 2012a). In this study we focus on the two official and most commonly used data sets and present the methodology we developed for analyzing different solutions, and deriving a quantitative estimate of the uncertainty due to the use of different solutions. Moreover, since the GRACE project centers have recently published a new release (RL05) of the monthly time series (though shorter and still incomplete), we have the opportunity to get an interesting overview of the difference between the RL04 and the new RL05 release for mass balance estimate.
Different approaches in extracting the mass balance from GRACE level 2 data can be another source of variability and uncertainties. Velicogna and Wahr (2005); Luthcke et al. (2006); Schrama and Wouters (2011);Horwath and Dietrich (2009) ;Sasgen et al. (2010); Barletta et al. (2008); and Sørensen and Forsberg (2010) all use different methods, and this could be part of the reason of the wide scatter in the results obtained in the literature since 2005. In principle, independent methods should produce the same results if careful calibration and cross validation have been carried out.
Our main goal is to provide as comprehensive and upto-date estimates as possible of the monthly mass changes at basin scale for the Greenland and Antarctic ice sheets, along with the secular trends of mass changes for the time period considered. For each of these products, we also provide a sound error estimate that takes into account most of the potential source of uncertainties, beside the data error.
In the following section we give a detailed description of the data sets (Sect. 2.1) and the corrections that we use in our analysis. In particular we describe the filtering of the data (Sect. 2.2), the de-aliasing correction (Sect. 2.3), the geocenter motion or degree-1 correction (Sect. 2.4) and the GIA correction (Sect. 2.5). Then we give details about the two methods, i.e. the mass inversion (Sect. 3.1) and the conversion (Sect. 3.2) methods, that we used to derive the mass changes for each basin according to the basin definition described in Sect. 3.3. For our analysis we need to define how we deal with errors and uncertainties (Sect. 3.4) and how we decided to compare and combine (Sect. 3.5) the many possible combination of solutions that we obtains in order to select our best estimate along with its total error.

GRACE L2 data releases
GRACE data are processed and provided to the scientific community by two GRACE project centers: Center for Space Research, University of Texas (CSR), and Geo-ForschungsZentrum, Potsdam (GFZ), with the Jet Propulsion Laboratory, California (JPL), providing the so-called validation solutions. Other research institutes such as the Centre National d'Etudes Spatiales (CNES) and Institut für The Cryosphere, 7, 1411-1432, 2013 www.the-cryosphere.net/7/1411/2013/ Geodäsie und Geoinformation, University of Bonn, also provide GRACE products. The different solutions are of similar quality but are not identical because of different processing strategies. The processing centers provide GRACE level-2 (L2) data (GSM product), which consist of monthly models of teh Earth's gravity field issued as fully normalized spherical harmonic coefficients, called Stokes coefficients. The analysis of the present work is based on the monthly GRACE L2 release 4 (RL04) and 5 (RL05) gravity field models provided by CSR (Bettadpur, 2007) and GFZ (Flechtner, 2007a). CSR and GFZ provide both formal and calibrated error estimates on the Stokes coefficients with their models, with the exception of CSR RL05, which at the moment do not include calibrated errors. This study is based on 113 CSR RL04, 105 GFZ RL04, 84 CSR RL05 and 72 GFZ RL05 monthly models, within the time span April 2002-February 2012.
The low-degree harmonic coefficients provided in the level-2 RL04 require some attention. The C 20 values show anomalous variability. Therefore, we replace these GRACE C 20 coefficients with estimates derived from satellite laser ranging (SLR) (Cheng and Tapley, 2004).
In the RL04 models, we also restore the C 20 , C 21 , S 21 , C 30 , and C 40 coefficients with the reference values of secular changes reported in the level-2 Processing Standards Document (Bettadpur, 2007).
GRACE data from all the processing centers are corrected at an early stage for short-term atmospheric and oceanic pressure variability. In the processing used to derive the new RL05, many improvements have been implemented compared to RL04 (Dahle et al., 2012a, b;Bettadpur et al., 2012).
Some examples of these improvements are significant upgrades to the background gravity model (GIF48), the GPS constellation being homogeneously reprocessed thereby improving the determination of GRACE satellites orbits, the reduction of the geometric bias in the K band ranging data (Horwath et al., 2011), the ocean tide model and the model for planetary ephemerides being changed and accelerometer biases estimated being improved. As a consequence, the solutions have been improved: with respect to the RL04, in the GFZ RL05 solutions the low-degree coefficients do not need to be replaced or restored as previously. However for the CSR RL05 solutions, we replace the C 20 coefficients, as specified in Technical Note 7, and we find that using original GFZ RL05 C 20 makes a significant difference in large basin (see Fig. SM12 in Supplementary Material) with respect to the use of CSR RL05 and GFZ RL04; so in order to analyze differences due to all the other degrees, we chose to replace also C 20 in GFZ RL05.
In addition, new atmosphere and ocean models have been used to correct the data. An overall visible result is a reduction of noise, especially the "stripes".

De-striping and filtering GRACE data
The GRACE monthly solutions are known to be affected by both noise (instrumental and statistical) and peculiar features in the form of north-south-oriented stripes (Swenson and Wahr, 2006). The latter represent systematic errors related to orbital trajectories and data processing strategies, and are apparent as spurious spatial correlations (Kusche et al., 2009).
We apply the de-striping method presented in Kusche et al. (2009) to the GRACE L2 solutions. This is a nonisotropic smoothing procedure, based on approximate decorrelation and successive regularization of the GRACE monthly solutions. The results provided here are obtained using a weak smoothing (smoothing parameter a = 10 12 ), i.e. using the DDK3 filtering method, which is available on the website of GFZ (http://icgem.gfz-potsdam.de/ICGEM/ TimeSeries.html). In particular the DDK3-filtered data used here correspond to a resolution of 240-330 km. The effective filter has a wider radius in the E-W than N-S direction and at high latitudes, the anisotropic filtering kernel becomes close to isotropic.
Data filtering and de-striping affect the mass balance at basin scale in different ways, depending on the method used to infer the mass variations. In the present work, these effects are effectively counteracted by a proper calibration which allows us to consider the variability due to filtering as a part of the method errors.
This also holds for the data cut-off degree (i.e. the data resolution): for all the GRACE data releases employed in this study we use the harmonic coefficients up to degree-60 (max degree for CSR RL04), but in the calibration phase we also tested the differences resulting from using higher harmonic degrees when available. A detailed description of the calibration procedure is available as a supplement.

GAC correction
The gravity field generated by atmosphere and ocean is computed from 6-hourly pressure field, which is the output of the respective models. These 6-hourly simulated gravity fields are called atmosphere and ocean de-aliasing products (AOD1B), and are used locally to correct the gravity field along the track of the satellites. This product is based on a combination of the ECMWF (European Center for Medium Weather Forecast) operational atmospheric fields and the baroclinic ocean model OMCT (Ocean Model for Circulation and Tides) (Flechtner, 2007b). After the processing, the monthly average of atmosphere and ocean gravity models is generated as a collateral product of GRACE data, and it is called GAC.
The GAC are provided to be used in combination with the GSM product to restore the atmospheric and ocean contribution to gravity variations whenever necessary for the user (GRACE L2 documentation Technical Note 04 and pg. 39 of AOD1B Product Description Document). They are obtained www.the-cryosphere.net/7/1411/2013/ The Cryosphere, 7, 1411-1432, 2013 as monthly average of the AOD1B, provided in the same format as of the GSM and, according to the instructions, they have simply to be added back to the GSM; no scaling or additional processing is required. So, since GSM represents the result of GRACE processing to "data-model corrections", and GSM+GAC is assumed to represent a good approximation of GRACE processing to "data", the GAC alone represent the difference between the L2 data corrected for the AO1DB and the one not corrected, i.e. of the effect of GRACE processing on the model. If it were not the case, it would not be possible to restore the atmospheric and ocean contribution to gravity by just adding GAC to GSM. The AOD1B RL05 is based on an ocean model with increased vertical and horizontal resolution and updated parameterization. AOD1B RL05 does not have some of the coastal artifacts that were present in the AOD1B RL04 (Bettadpur et al., 2012). As soon as the new RL05 was released, it was possible to compare the GAC products of RL04 and RL05, and the difference is clearly visible (see Fig. 1). With respect to GAC-RL05, the GAC-RL04 time series shows a sudden jump in 2009 at some locations in the ocean, especially close to the coast of Antarctica (Fig. 1). The jump is artificial and is caused by an offset in the atmospheric pressure fields that are used to force the ocean model (H. Dobslaw, personal communication, 2012), and affects the monthly mass change estimate in some basins. The difference between the GAC-RL05 and the GAC-RL04 also shows pure trends in some basins instead of the jump. Notice that the difference between GSM-RL04 and GSM-RL05 is not due only to this error, but by analyzing the GAC we can appreciate its specific contribution.
We derive a correction for RL04 data which removes the artificial jump in 2009 and trends. One simple way is to use the differences between the monthly GAC-RL05 and GAC-RL04 models that we call GAC[04-05], for the common months. However, this shortens our RL04 time series to that of RL05 (January 2004-December 2010. This simple correction is used to compare the RL04 and RL05. In order to exploit the full length of the RL04, we derive a more refined correction, M GAC (t), which is derived by fitting a function (step-wise plus a linear term) to the monthly GAC[04-05], for each basin. For the M GAC (t) we also compute the standard deviation with respect to the GAC[04-05] to be used as monthly error. Furthermore, we compute the contribution of the M GAC (t) correction to the trend.

Degree-1
The geocenter motion, in harmonic coefficients, is represented by the degree-1 (C 10 , C 11 , and S 11 ) components, whose relation with the geocenter displacement in Cartesian coordinates (X, Y, Z) in the CF frame is found for example in Cretaux et al. (2002). These components are in principle not detected by GRACE, but they need to be included in the mass derived by GRACE (Chambers, 2006). Neglecting the degree-1 would lead to a systematic error in mass variations, especially for Antarctica.
The choice of the correction for the degree-1 is still an open issue, and it is subject to much debate. In the GRACE official website (NASA GRACE Tellus http://grace.jpl.nasa. gov) the degree-1 is that from Swenson et al. (2008). However, many other geocenter motion time series are available, calculated for GRACE and for other purposes (Wu et al., 2012). Those geocenter monthly time series differ in both amplitude and phase. It is also important to take into account whether a given time series includes the GIA degree-1 trend or not. For instance, Swenson et al. (2008) already have the degree-1 term corrected for GIA (based on the ICE5g-VM2 GIA model) i.e. they remove the GIA trend during their processing, while the SLR-derived degree-1 estimates (Cheng et al., 2010) are not corrected for GIA, i.e. it contains the full observed signal.
The Cryosphere, 7, 1411-1432, 2013 www.the-cryosphere.net/7/1411/2013/ The SLR-derived geocenter motion is reported to be the most precise (Wu et al., 2012). However the degree-1 derived from GRACE and ocean models (Swenson et al., 2008;Riva et al., 2012) or GRACE and GPS (Rietbroek et al., 2012a) have a smaller amplitude than the one based on SLR. Klemann and Martinec (2009) show the large variability of GIA degree-1 trends with respect to the viscosity of the lower mantle. Since one of the main problems for Antarctica is that GIA is poorly constrained there (both for the ice history and viscosity model), we are left with a very wide choice of possible degree-1 time series which also contribute to the variability of our estimated monthly mass changes as well as their trends.
To deal with this issue, we build a sensitivity kernel for the degree-1 contribution at basin scale; i.e. we estimate the mass variation for each basin in Gt for 1 mm of variation in each of the X, Y , and Z components of the geocenter motion.
This provide for each basin a linear transformation that maps the geocenter motion, the displacement of the center of mass with respect to the center of figure, into mass variations at basin scale. It is also possible to build a "reverse" kernel that maps the effect of the gravity variation due to a mass change in a specific basin into the displacement of the geocenter, but this is not relevant for the present work.
For each basin, the degree-1 correction, in terms of mass change for each component (X, Y, Z), is obtained by multiplying the sensitivity kernel (Tables 1 and 2) by the X, Y , and Z geocenter time series (in mm) or their trends (in mm yr −1 ). The total degree-1 correction is the sum of the three components.
After deriving a time series representing the degree-1 correction for each basin, we simply add it to the time series obtained from the (uncorrected) GRACE data. This simple procedure allows any user of our GRACE-derived time series to easily change the degree-1 correction using our kernel together with new and more updated geocenter motion time series (more detail in the Supplement).
In the present work we use three independent geocenter motion time series to derive the degree-1 corrections: 1. Swenson et al. (2008), for seasonal and trend component (SW).
3. SLR - Cheng et al. (2010), for the seasonal component, and for the trend we will apply a global GIA correction (Klemann and Martinec, 2009), (SLR).
These are among the most up-to-date time series, each of them representing a different method to infer the degree-1.
For the geocenter motion monthly solution we compute the monthly average with its standard deviation between the three detrended time series (SW, RR and SLR). We use this average and the degree-1 sensitivity kernel to compute our Table 1. Degree-1 sensitivity kernel for Antarctica. The first column indicates the basin number (or region), the second its area, and the other three columns indicate the variation in Gt due to 1 mm of variation in the geocenter coordinates X, Y and Z (in Gt mm −1 ). Two spherical caps of the dimension of the polar gaps (PG) and three macro regions, Antarctic Peninsula (AP), West (WA) and East Antarctica (EA), are also computed. preferred degree-1 correction in mass balance time series (for each basin). For the geocenter motion trend we use an average (with its standard deviation) between four trends as in Table 3. The first (first line of the Table 3) Table 2. Degree-1 sensitivity kernel for Greenland. The first column indicates the basin number (or region), the second its area, and the other 3 columns indicate the variation in Gt due to 1 mm of variation in the geocenter coordinates X, Y and Z (Gt mm −1 ).  (2009)

Glacial isostatic adjustment
On the NASA GRACE Tellus site (http://grace.jpl.nasa.gov/ data/pgr) a global GIA model (Paulson et al., 2007) is provided. This model is recommended as the best model, and it is associated with a uniform uncertainty of ±20 %. This model is based on the global ICE-5g ice history of Peltier (2004), and a simplified VM2 Earth model: a 4-layered approximation to VM2 viscosity profile. This model also includes the rotational feedback based on the formulation of polar wander described by Mitrovica et al. (2005). Through the sea level feedback, the center-of-mass motion is also taken into account. For Antarctica, another GIA model is the IJ05 (Ivins and James, 2005), which has often been employed as a correction in mass balance studies (Horwath and Dietrich, 2009;Chen et al., 2006a;Gunter et al., 2008Gunter et al., , 2009). Thomas et al. (2011) show that traditional GIA models (ICE5g-VM2, IJ05) predict an uplift that is too large, when compared to the uplift rates measured by GPS, and they indicate that the Riva et al. (2009) empirical model is the one which gives the best fit. Riva et al. (2009) find a GIA pattern for Antarctica by combining GRACE and ICESat data, and they find a smaller signal than the one of traditional GIA models (ICE5g-VM2, IJ05). Wu et al. (2010) also find a GIA pattern empirically by performing a global inversion with GRACE data, and the signal over Antarctica is smaller than ICE5g but higher than Riva et al. (2009). So there are now many indications that traditional GIA models overestimate GIA signal over Antarctica, and for this reason much work has been done in the last years taking into account newly available data to constrain the LGM (Last Glacial Maximum) ice models.
However, these new ice models are still not publicly available; hence we decided to derive proper GIA corrections for Antarctica, taking into consideration the most recent evi- Table 3. Trend for the geocenter motion. The first line (SW) reports the trend computed for Swenson et al. (2008). The second line (RR) reports the trend found in Rietbroek et al. (2012b). The third and fourth lines use the trend computed using the SLR time series (X = −0.131, Y = 0.352, Z = −0.637) minus the GIA geocenter motion given in Wu et al. (2012) for ICE5g/IJ05/VM2 (X = −0.12, Y = 0.24, Z = −0.48), and minus the Klemann and Martinec (2009) 0.096 0.128 0.150 dences, and we use them together with the Riva et al. (2009) empirical GIA model. One reason for the GIA overestimate can be the excess of ice melting during the LGM for the traditional ice models (ICE5g, IJ05), which clearly violates the geological evidence on the ice history (Todd et al., 2010;Ackert et al., 2011;Mackintosh et al., 2011). However due to the tradeoff between ice history and solid Earth response in GIA, another way to predict a smaller present-day signal in Antarctica is to choose a lower-viscosity profile, especially in the western part where the majority of the deglaciation took place. Barletta et al. (2008) shows that for low viscosities, the GIA corrections are smaller. In fact with respect to traditional viscosity profiles, like VM2 Peltier (2004), lower viscosities produces faster relaxation so that most of the mantle flow now is over, and the residual present-day signal is smaller. Changing the ice history or changing the Earth's response through the viscosity profile produces different GIA corrections that could be distinguished if enough constraints, both in space and time, were available. However the apparent mass changes due to a GIA correction can be effectively reproduced with a suitable choice of the viscosity profile. Since at the moment there are no new revised ice models available, we propose two alternatives, using two traditional ice models (ICE5g-VM2, IJ05) with low-viscosity profile.
We also tested the empirical GIA pattern extracted in Barletta and Bordoni (2009) and we verified that it produces mass balances similar to those produced by ICE5g-VM2. The GIA pattern obtained in that work was an upper bound of possible GIA patterns (the positive signal due to GIA was not distinguished from others due to possible present-day mass accumulation). So we decided to use ICE5g-VM2 as an upper bound of the GIA contribution.
In Greenland the GIA signal is small with respect to the total trend but, depending on the model, it is not negligible.
The Cryosphere, 7, 1411-1432, 2013 www.the-cryosphere.net/7/1411/2013/ For this region two commonly used models, ICE5g and ANU model (Fleming and Lambeck, 2004), especially for some basins, predict signals with opposite sign. We use both these two model because the difference between them ensures a reliable confidence interval for GIA prediction in Greenland. GIA corrections are based on 5-layered incompressible Earth models without lateral heterogeneities. We compute the solid Earth Green's functions (Love numbers) according to the analytical approach described in the benchmark of Spada et al. (2011) with the benchmarked code TABOO (Spada et al., 2004). The sea level equation is solved self-consistently with the pseudo-spectral approach (Mitrovica and Peltier, 1991) implemented in the optimized code TSec01 developed for Barletta andSpada (2011, 2012b) and benchmarked in Spada et al. (2012). We include the degree-1 and perform the computations in the CM reference, up to degree-128, without rotational feedback. However we remove the degree-1 to compute the apparent GIA mass changes, i.e. the GIA correction.
To summarize, we apply five different GIA corrections: 1. The revised version of Paulson et al. (2007) ICE5g-VM2 GIA model, which uses a compressible fully layered Earth model (Geruo A/Paulson). This model is used for both Antarctica and Greenland.
3. The Riva09 empirical GIA estimate for Antarctica  4. ICE5g-LV with a low viscosity (LV) of the upper mantle for Antarctica 5. IJ05-LV with a LV of the upper mantle for Antarctica By varying, within a chosen range, the lithosphere thickness and the viscosity of the upper and lower mantle we compute several corrections. The parameters are chosen within a range that gives, at basin scale, a GIA correction centered around the one obtained with the Riva09 empirical GIA . The viscosity range, in particular, is in agreement with the one inferred by Ivins for Shepherd et al. (2012). Then we compute the average and the standard deviation for each of the ensemble ice + earth models, the 2, 4 and 5 in the list above. The Earth parameter ranges for lithosphere thickness (LT), upper-(UMV) and lower-mantle viscosities (LMV) are  Due to the short timescale considered here, the GIA correction only affects the linear trend, and not the seasonal component, nor short-term changes (e.g. accelerations) of the time series. We estimate the uncertainties for all these proposed models by different approaches: the GIA uncertainties for ICE5g-VM2 are computed as in Barletta andSpada (2011, 2012b); the Riva09 model is provided with error estimates, and we propagate these to obtain error estimates of the final mass changes.
In our final preferred trend we use as GIA uncertainties the maximum on GIA uncertainty among the four models. Tables 4 and 5 report the values used in this study.
www.the-cryosphere.net/7/1411/2013/ The Cryosphere, 7, 1411-1432, 2013 Table 5. GIA trends for Greenland. The first column is the basin number, the second its area (in 10 6 km 2 ) and from the third column the values are in Gt yr −1 , and they report the values for the three models, namely ICE5g-VM2 compressible with rotational feedback, ANU, and ICE5g-VM2 incompressible without rotation. The last column reports the maximum uncertainty among the models.

Methodology
Different methods for deriving mass change estimates from GRACE data should ideally agree, when using the same input data. For gravity-derived mass changes the major difference between one method and another is the leakage treatment. In fact with low-resolution measurements (such as GRACE), discriminating the signal coming from inside the region of interest from the one coming from outside is a challenging problem. By assuming that the region of interest produces a much stronger signal than the surrounding, many methods treat the leakage by enlarging the integration areas (or the averaging kernel). Horwath and Dietrich (2009) deal with this problem in quite a systematic way and show the errors in mass changes associated with leakage.
The two methods that we use are different and independent and treat the leakage problem in different ways. Leakage is not considered explicitly as a factor affecting the scatter of the solutions. Leakage treatment is a part of the method, and it is dealt with during the calibration of the two methods, so it is implicitly taken into account in the variability, as a part of the accuracy error. The leakage between adjacent basins is an issue, but it is strictly related to the resolving power of the GRACE data, and it cannot be effectively avoided, without introducing external constraints. However, also in this case the two methods deal with this issue in different ways but produce very similar results, so we assume that also this contribution is part of the accuracy error.

Method 1: mass inversion
The direct point mass inversion method used for determining the monthly mass changes from the monthly GRACE data is based on Forsberg and Reeh (2007) and Sørensen and Forsberg (2010). Prior to the inversion, corrections and filtering are applied to the GRACE data as described in the previous section. The elastic response of the solid Earth to present-day ice mass changes involves changes in the gravity field. This must be removed from the GRACE data before deriving the surface mass changes from the observed gravity changes.
The elastic corrections based on the model of Farrel (1972), as also presented by Wahr et al. (1998), are used to make "reduced" gravity disturbances for the observation equations in the inversion. The inversion is performed on a set of N y observations y = {δg k }, k = 1 . . . N y , i.e. gravity disturbances at the altitude of the satellites, and solved for a point-like mass ensemble x = {m j }, j = 1 . . . N x located at coordinates (θ j , φ j ) which define the solution area. The linear problem y = Ax is solved using a generalized least squares inversion with Tikhonov regularization x = (A T A + λI) −1 A T y, where λ is a smoothing parameter and the observation matrix A is built upon the attraction of a point mass of the sphere to the measured gravitational attraction at the orbit level by δg k = Gm j a 2 (h + a) − a cos ψ kj r 3 kj . (1) Here G is the gravity Newton constant, a is the mean radius of the Earth, h is the height of the observation, and r kj and ψ kj are the distance and the angle, respectively, between the observation δg k and the solution point m j .
The inversion method has been refined, optimized and calibrated for this work. First we optimized the solution area using icosahedron-based grids (Tegmark, 1996) with disk elements of almost equal area. Then we improved and calibrated the solution area in order to reduce as much as possible the leakage from the signal outside the region of interest. The inversion method assume that the gravity signal is negligible outside the region of interest; otherwise such a signal would be forced to be part of the ice mass changes in the solution area. Once this assumption is verified by applying the method on synthetic data, we find that it is able to recover up to 99 % of the mass. One strategy for mitigating the effect of ocean mass changes being erroneously modeled as ice sheet changes is to force the ocean signal to be zero, especially in the far field. The signal outside the region of interest that is farther than some hundreds of km (300 to 500) from the boundary of the solution area can be forced to zero (zero mask) without compromising the signal of interest. Another strategy is to build a complementary solution area (CSA) around the primary one. The CSA is a belt around the original solution area, but separated by a gap of some hundreds of km (300 to 500), and it accounts for the signal outside the original solution area. We used a combination of the two above strategies, and we calibrated the gap for CSA and the distance for the zero mask. The parameters we chose with the calibration allow us to retrieve about 98 % of different kinds of synthetic signals. The regularization parameter (the smoothing parameter λ) used for this study was also chosen after calibration using a synthetic data set (more The Cryosphere, 7, 1411-1432, 2013 www.the-cryosphere.net/7/1411/2013/  (Tegmark, 1996), using disks of about 40 and 20 km radius for Antarctica and Greenland. detail can be found in the Supplement). Note that we do account for the degree-1 also in our calibration process because it does alter the calibration. Figure 2a shows the solution area for the inversion of the trend over Antarctica. Figure 2b shows the solution area for the inversion of the trend over and around Greenland, and the point-like mass units, used as solution area for the inversion, are also visible. The closest surrounding ice-covered areas (Ellesmere Island) are included, to reduce leakage. In fact Ellesmere Island has a strong trend and is so close to the northwest Greenland that it cannot be treated as the rest of the surrounding sea and islands around Greenland.
Once we obtained a mass grid from the inversion scheme for each month we integrate over each basin's area to derive the total mass change for Antarctica and Greenland. The mass estimates of each of the basins is a summation of the point mass changes within each basin mask definition (Sect. 3.3).

Method 2: conversion and integration
From the GRACE spherical harmonic coefficients, surface mass density in water equivalent (w.e.) is generated, as presented in Wahr et al. (1998). In order to take into account leakage, we integrate over a region more extended over the sea than the original area of interest, as done in Barletta et al. (2008). We use (for each basin) three integration areas which extend over the sea for 100, 200 and 300 km, and then we derive a weighted average of the three integrals. The largest weight is on the results obtained with the integration area which extend for 300 km from the coast, since this is the resolution of the data, and in this way we are able to reduce the leakage, and recover most of the signal.
This method is basically a simplified version of the Velicogna and Wahr (2005) and Horwath and Dietrich (2009) method, and it does not need a suitable averaging kernel nor a rescaling factor. The latter is usually used to retrieve the signal lost either for the filtering, the cut-off degree (data resolution) or the leakage.
The leakage for this method is only considered coming from the land that we are analyzing, so we do not get rid of the leakage coming from the ocean into the land. The latter is a strong a priori assumption, which is checked a posteriori by comparing the results, at basin scale, with those obtained with the inversion method, which do not rely on the same assumption (more details can be found in the Supplemental Material).

Basin definitions and resolution
For the Greenland ice sheet we use the basin definition presented in Hardy et al. (2000). The Greenland ice sheet is divided into seven major basins, which are shown in Fig. 3a. The Antarctic ice sheet is divided into 27 major basins as shown in Fig. 3b, which is the same as Zwally et al. (2012).
Both methods use these basin definitions even if the sampling (their resolution and geometry) is performed on different grids and so the contour of each region is not exactly the same, especially for the smallest basins. of about 40 and 20 km radius for Antarctica and Greenland, respectively (Fig. 2). Method 2 uses Gaussian grids of 128 × 256 cells (latitude × longitude). Notice that GRACE resolution is about 300 km and both methods use grids with higher resolution for the mass reconstruction, but the results depend strongly on data resolution and just slightly on the method resolution, i.e. the internal working resolution (and format) for each method.

Uncertainty estimates
As already mentioned, the mass estimates are associated with several sources of uncertainties. Some are coming from the data errors and others introduced through the algorithms employed, sampling grid sizes or smoothing and other parameters (in the inversion).
We derive the uncertainties which are related to the data errors provided directly with the GRACE monthly models by using a Monte-Carlo-like approach in which 100 simulations are performed. The simulations are created from Stokes coefficients drawn from normal distributions with zero mean, and the calibrated standard deviation provided with the GRACE level-2 data (Tscherning et al., 2001).
We deal with both precision and accuracy errors in our final results. The precision error accounts for the statistically distributed random error around one average value. The accuracy error accounts for how much the expected value deviates from the "true" value.
For the precision error we provide the 95 % confidence interval (2σ ) propagated from the data. For each method by using synthetic data (as in the calibration procedure) we as-sess its accuracy so that it can be added to the precision error to give the total error.
We find that the accuracy error only associated with the differences of the methods is about 2 %, which is much smaller than the difference found by the use of different data sets. So in our final estimates we neglect this small contribution, and our accuracy error is derived from the difference coming from the use of different methods and data set.

Best estimate and comparison strategies
The strategy to deal with many possible combinations of data, methods and corrections is relevant on the one hand to minimize the number of steps and so the possible sources of (human) mistakes, and on the other to make the procedure clearer and so easier to be verified and reproduced.
Both methods for computing mass balance are linear, so we estimate each of the corrections separately and combine them with the time series only at the last stage. The same holds for all of the trend corrections, especially the GIA correction which have been treated separately and added back in the last stage. This overall strategy allows us to analyze in detail the relative weight of each of the components in our final results. The combination of data and corrections in the final result M F (t) is a simple sum: where the M x (t) is the mass changes time series and the subscript x indicates the source of each contribution: the uncorrected data (data), the GIA correction (GIA), the degree-1 correction (Deg-1) and the GAC correction (GAC). , 7, 1411-1432, 2013 www.the-cryosphere.net/7/1411/2013/ For each basin, for each method and each data set we compute the mass balance (with its propagated uncertainties) for each available month. So for each basin we have four time series using GRACE data RL04 and four using the RL05 ones. We computed for each basin three time series for the degree-1 correction and two for the GAC correction (the difference GAC[04-05] and the fitted function). One of the degree-1 time series is provided on a much denser timescale than monthly, so we perform an averaging procedure. The majority of the time series do not share the same exact sequence of months, so we select only the common months for each comparison and combination in order to maximize the number of available months in each case.

The Cryosphere
Therefore the combination of the methods, models and corrections gives a set of results (mass variation time series) that are supposed to (indirectly) measure the same phenomenon and have to be considered equally reasonable. Each measure has its own formal error, but since we do not know which is the true value, the scatter of the results can be used to obtain a best estimate, and the related accuracy error, i.e. the relative distance from the true solution, by quantifying their differences. Practically, we want on the one hand to be able to compare different time series, and possibly to estimate their similarity; on the other, we want to extract from this comparison a reasonable estimate of the best value and its accuracy.
Comparing two or more time series is not trivial and the conclusion that can be drawn depends on what we are interested in. For example, if the main interest is in the trend, often all the rest of the signal present in the time series is neglected. If the focus is instead on seasonal phenomena, the straight way is to extract the seasonal signal, and compare the results. In this case, it has also to be decided how to deal with some potential scaling factor and or phase shift, for example. In a statistical approach suitable for data sets where the real form of the expected signal is not obvious, a natural choice is to focus on correlation indices and related properties. These, and others, are of course legitimate approaches, each with its own advantages.
For the present work, our aim is to find an effective way to compare scattered time series, and to extract an optimal solution (the most accurate one). We choose to follow an approach that allows us to extract different information at the same time. It is flexible enough to be used to compare two time series without many assumptions on the expected signal, but when the time series are particularly simple, like in presence of a strong dominant trend, or of an extremely good correlation, it allows for a straightforward qualitative and quantitative interpretation. It has been chosen also to address both global properties (regularity, long-term behavior) and point-like similarities by keeping them separate at the same time. Before discussing the method in formulas, let us give a simple description.
Given two or more time series to be compared and/or combined to produce a best estimate, the first way that comes to mind is to perform a simple average and to extract the standard deviation. However, different linear trends (and offsets) in time series act as a systematic bias. This spoils the statistical assumption which the standard deviation is based on, i.e. a Gaussian distribution around the mean value. For this reason we separate the discussion of global properties (e.g. trend estimate) from the monthly solutions.
Let us assume that two time series are extremely well correlated, but differ by a scaling factor. This is what often happens when comparing GRACE time series differing only in the treatment of leakage (for example via different scaling factors like in Velicogna and Wahr, 2006). Even if a visual inspection would suggest that they are very similar, a comparison of the trends, and possibly other overall properties, would show that they differ in amplitude. The similarity is instead found in the correlation.
We follow an alternative approach: first we try to match the overall (global) features of the two time series, i.e. by scaling one with respect to the other (or both with respect to their average). In this way, the rescaled time series would be directly comparable: their difference (the residue) would describe the point-by-point (or time-by-time) agreement, while the scaling factor represents the overall agreement between the two series (the closer to 1, the better is the overall agreement). If the two time series are extremely similar except for a scaling factor, the result of our method would be that the overall difference is large, while the point-like difference is negligible.
Another example shows why a direct comparison of trends in the time series could be not satisfactory, or even misleading -assuming that two time series are highly correlated (so that they seem to portray the same phenomena) and have the same trend, but with very different amplitude in the seasonal contribution. By comparing the correlation we would say that they agree well. By comparing the trends, we could say that they agree perfectly for this part of the signal. But we think they cannot be considered as being in good agreement, so in this case our approach would conclude that the point-like features are very similar, but that the two series are not in overall good agreement, and that the apparent agreement of the trends alone is misleading.
So, we decided to implement this analysis, as follows. Two time series M 1 (t) and M 2 (t) can differ for an overall factor m (regression index) and an overall offset q, and they can show also other kinds of time-dependent differences (t), e.g. M 1 (t) = mM 2 (t) + q + (t). We first identify the regression index and the offset with a least square fitting procedure and remove it in order to analyze just the monthly difference between the time series. When we analyze the difference for each basin between two time series, we compute Following the previous discussion, the overall agreement is represented by the regression index m, while the delta www.the-cryosphere.net/7/1411/2013/ The Cryosphere, 7, 1411-1432, 2013 (Eq. 3) represents the "time-by-time" agreement. Moreover we add the correlation coefficient that indicates whether the regression index is meaningful or not. When the correlation coefficient is above 75 % the (t) component is negligible and so the regression index represents the overall scaling factor between the two time series. The scaling factor m can be due to both trend and amplitude in the periodical component, and when the trend is dominant with respect to the seasonal signal (Assumption RI), the relative difference in the trends can be estimated by where T 1 and T 2 are the trends of the series M 1 (t) and M 2 (t), respectively. When the above assumption (RI) is not verified, Eq. (4) represent the relative ratio between amplitude in seasonal signals, alone or combined with a trend or other longterm signals. If the regression index m is close to one, i.e. δ T 21 /T 2 and δ T 12 /T 1 are close to zero, then the two time series do not need to be scaled relative to each other.
For averaging the time series we follow an analogous strategy: we first find the average between N time series, M i (t) with i = 1. . .N , as M A (t) = i M i (t)/N, then we find, for each of the series M i (t), the regression index with respect to the average M A (t), i.e. the m i and the q i parameters, so that M A (t) = m i M i (t) + q i + i (t) for each i. Then by rescaling each solution with these parameters we find for each series the monthly residual standard deviation with respect to the average (assuming normally distributed values): where δ A (t) is considered as the accuracy error when averaging between time series obtained with different methods and different data sets. The monthly precision error (the 95 % confidence interval) is 1.96 × σ A (t) computed from the monthly square sum of the σ i (t) of the time series. So the total error on M A (t) is Extending the definition of Eq. (4) for δ T 12 , we define the quantity δ T A (assuming normally distributed values) as where T A is the trend of the average. Under the above assumption (Assumption RI) we can use δ T A as an estimate for the trend accuracy. For those time series (i.e. basins) that have a low correlation index, the comparison with the average series M A (t) will have a better correlation index, thus making the related regression index more meaningful. Moreover when the above assumption (Assumption RI) does not hold, Eq. (7) represents the difference that the trend will have after applying the scaling factor m. Therefore, as explained above, by using an overall scaling factor to derive accuracy error we are accounting for the overall differences between time series.
The trend on the time series is computed with a weighted least square fitting procedure using a function composed by a linear trend, annual and semi-annual term, and a constant. As for the monthly weight we use the monthly error of the time series. The precision error (the 95 % confidence interval) 1.96 × σ T A for the trend is computed using the variance of the fit.
Note about q: all our GRACE-derived mass changes are computed with respect to the same reference gravity model, so they should have small or no q offset from one another. We find that a small offset difference between data obtained with different release and methods exists.

Results and discussion
All the monthly time series shown in the following pictures are not corrected for the trend contribution of degree-1 and GIA, but these are accounted for when analyzing the trend.

Method vs. data
We perform a cross validation of the two methods using both CSR and GFZ GRACE data, 113 CSR RL04, 105 GFZ RL04, 114 CSR RL05 and 107 GFZ RL05 monthly solutions (see note in SM for the RL05 time span), from April 2002 to February 2012. We compared all the 27 + 7 basin time series, and in addition to the discussion of the results for all of them, we show one example for each region, Antarctica and Greenland.
A general good agreement between the two methods for each of the data sets (CSR and GFZ) for each basin is clear by a simple visual inspection (pictures for each basins are in the Supplement Material). This agreement is particularly clear in the Amundsen sector (Fig. 4a, Antarctica basin 21) and in southwest Greenland (Fig. 4b, Greenland basin 6), which we choose as examples. By visual inspection the time series for the inversion ("i", solid lines) and the conversion method ("c", dash-dotted lines) overlap over almost the entire time interval, for both data sets, CSR (blue) and GFZ (red) lines. For each basin, we quantify the agreement in the time series, regardless of the trend as explained in Sect. 3.5, by using the point-to-point difference δ 12 (Eq. 3), computed for each basin with respect to the average monthly error (Fig. 4c). We find that the differences between the two methods (lightpurple bars) are much smaller than the differences (lightgreen bars) between the use of CSR and GFZ, for almost all basins: basins 19 and 25 show the same amount of difference The Cryosphere, 7, 1411-1432, 2013 www.the-cryosphere.net/7/1411/2013/   Comparison between the two methods and two data sets (CSR and GFZ): in the legend the inversion method (Method 1) is indicated by "i" and solid lines, and the conversion method (Method 2) is indicated by "c" and by dash-dotted lines. The use of CSR is indicated by light-blue and blue lines, and the GFZ by light-red and red lines. Each of the small dispersion graphics shows the time series obtained with inversion versus conversion methods (purple square), and the use of CSR versus GFZ data set (green dots). For each of these couple of time series the value of their regression index m (as in Eq. 3) is indicated in square parentheses. In panel (c), the vertical axis indicates the basin (number, the region (EA, WA, AP) and area in 10 6 km 2 ) in descending order from the largest. For each basin, differences are plotted for the two methods (purple) and for the two data sets (green). The light colors represent (the quadratic sum of) the monthly difference with respect to (the average of) the monthly errors. The normal colors represent 1 minus regression index m, i.e. δ T 21 /T 2 as in Eq. (4), and the grey bar is the error on the same trend as Figs. 8d and 9d with respect to the trend. The light-blue and yellow bars represent 1 minus the correlation coefficient for the use of different method and data set, respectively. When the correlation is below 75% (i.e. the related bars are larger than 0.25) the associated regression index is not meaningful. between the two methods and the use of different data sets. Basin 25 is the smallest, so either resolution or leakage can give noisy signal. Basin 19 has a very small trend and it is very close to basin 20, which has one of the strongest trend, so resolution or leakage can give again noisy signal.
The different combinations of methods and data sets also result in different overall factor among the mass time series. This difference is quantified by the linear relation, the regression index m and offset q parameters (as explained in Sect. 3.5), between two time series (shown in the small dispersion graph inside each of the two plots of Fig. 4a, b). The two data sets (green dots) in the Amundsen sector produce a difference in regression index of 11 % and the two meth-ods (purple squares) only 1 %. For southwest Greenland the differences in regression index are closer, i.e. 6 % vs. 8 %.
We analyze for each basin the relative difference δ T 21 /T 2 of Eq. (4), (i.e. 1 − m, scale on the bottom), (Fig. 4c) caused by the two methods (purple bars) and the two data sets (green bars). We also take into account the average correlation coefficient obtained by comparing the two method (light-blue bars) and the two data sets (yellow bars), and in Fig. 4c we find it convenient to show its distance from one (scale on top). We find that the comparison of the two methods gives for all basins except basin 19 an extremely high correlation (light-blue bars close to zero), and therefore the associated analysis based on regression index is meaningful. The correlation among time series obtained with the two data sets Each of the small dispersion graphics shows the time series obtained with RL04 versus RL05 with the use of CSR (blue) and GFZ (red). Note that GFZ RL05 has several months missing in 2004 and that is the reason for the anomaly in the plot. In panel (c), for each basin differences between RL04 and RL05 are plotted for the use of CSR (blue) and GFZ (red). The light colors represent (the quadratic sum of) the monthly difference with respect to (the average of) the monthly errors. The normal colors represent 1 minus regression index m, i.e. δ T 21 /T 2 as in Eq. 4, and the grey bar is the error on the same trend as Fig. 8d and Fig. 9d with respect to the trend. The very light-blue and dark-pink bar represent 1 minus the correlation coefficient for the use of the two releases of CSR and GFZ, respectively. When the correlation is below 75% (i.e. the related bars are larger than 0.25), the associated regression index is not meaningful. NOTE: the difference in the time span of the RL05 time series is due to an update in data; see note in the Supplement.
is instead in several cases below 75 % (yellow bars larger than 0.25, scale on top), and in all those cases the associated relative difference δ T 21 /T 2 is meaningless. For small basins (area, indicated in the vertical axes of Fig. 4c, smaller than 0.1 million km 2 ) the two methods result in differences in regression index between 20 % and 30 %, larger than the error on the trends (the grey bars) and than the difference produced by the two data sets. In particular the conversion method produces trends with absolute values larger than the one obtained with inversion method. We also performed analogous analysis directly on the differences in trends (Fig. SM10) and we can confirm the same kind of differences for the use of different methods. Using different data sets, regression index and direct trend analysis show the same kind of differences only for basins with high correlation coefficient, namely (from larger to smallest) 3, 6,7,4,18,21,22,20,23,26,25 for Antarctica, and all but 2 for Greenland. All of those basins are experiencing strong trends and for all of them the use of CSR gives absolute values of trends slightly larger than GFZ.

RL04 versus RL05
The differences between the use of RL04 and RL05 data are also related to the difference GAC[04-05] in their de-aliasing product (green line in Fig. 1)  accounts for the possible variability during the full GRACE period.
In order to compare the two distributions as discussed in Sect. 3.5, we correct our RL04 solution using the GAC[04-05] and we compare it with the RL05 just for the inversion method (Fig. 5). Also for this comparison, for each basin, we quantify the agreement in the time series by using the point-to-point difference δ 12 (Eq. 3), with respect to the average monthly error ( Fig. 5c the light-blue and light-red bars). The results obtained from the two releases are quite different. We computed analogous differences using detrended time series (Fig. SM11) and we confirm that we obtain very similar distributions of average monthly differences even though in some cases they are larger. The average monthly differences are larger than 40% of the monthly total error and this is compatible with our accuracy error. Comparison between RL04 and RL05 is also performed with the correlation coefficients, and where those coefficient are very high (values close to zero, scale on top, for very-light-blue and dark-pink bar in Fig. 5c) also the regression coefficients give us information about differences (blue and red bars), i.e. in basin 3,6,7,4,18,21,22,20,23,26,25 for Antarctica and all except basin 2 for Greenland. In those basins' differences, the δ T 21 /T 2 represents also differences in trends (Fig. SM11) which are comparable to or smaller than the total error on the trend except for basin 18, where GFZ-RL05 shows differences with respect to the GFZ-RL04 larger than the error on trends. In the monthly differences, in the correlation coefficients and in the trend differences, GFZ show in several basins larger differences than CSR. The explanation could be found in our previous analysis on the differences between CSR and GFZ release 04: in the new release 05, once we replace the same C 20 in both solution they give much more similar results.

Errors in each monthly solution and degree-1 contribution
Yet another important contribution to the mass change variability comes from the degree-1 correction (see Sect. 2.4). The sensitivity kernel for degree-1 (Tables 1 and 2) can also be interpreted as the relative weight of each component X, Y and Z of the geocenter motion on each basin. From Table 1, it is clear that basin 17 in Antarctica is the most affected by degree-1, but it is also the largest basin, and it is quite straightforward to see that the impact of the degree-1 is proportional to the area of the basin. The three degree-1 time series obtained with this sensitivity kernel (Table 1 and 2) and the three geocenter motion detrended time series have quite similar phase but different amplitude (two examples in Fig. 6). The mass changes time series are the sum of four contributions (Eq. 2, Fig. 7). When estimating the errors on the monthly mass changes, we neglect the GIA because it only gives a contribution to the trend. For taking into account the effect of using different data sets, we use the monthly average (Sect. 3.5) and its error (Eq. 6), being both precision (blue band) and accuracy (light-blue band) errors. We then add the monthly average of the degree-1 and its standard deviation (green band), and as a last step, we add the GAC correction M GAC (t) and its standard deviation (yellow band). We find that the degree-1 variability, represented by the monthly standard deviation with respect to the monthly average (green band), has a considerable impact on the total error as can be appreciated in Fig. 7c, where it is presented as average per basin with respect to the average of the total error. The contribute of degree-1 uncertainty to the total error ranges from 20% for the smallest basin up to 40% for the largest. The monthly variability of the GAC correction is quite small but not negligible (   of the accuracy error (light-blue bars) represents the variability in the use of different methods and different data sets, and its impact on the total error is comparable with the degree-1 variability. The three sources of monthly variability together (the accuracy, the degree-1 and the GAC) are for almost all the basins larger than the precision error propagated from the data alone; i.e. the blue bars are almost all below 50 %.

Uncertainties and contributions to the trends
Basins 20, 21 and 22 (in front of the Amundsen Sea) have the largest trend in Antarctica (Fig. 8c). The Amundsen sector (basin 21 and 22) has also the lowest relative errors (violet bars in Fig. 8a). In the Amundsen sector also all the other contributions different from the GRACE data (the blue bar, Fig. 8a) have a very small impact, namely the degree-1 (skyblue bars, Fig. 8a), the GIA (green bars) and also the GAC (orange bars) corrections. Due to its character, the GIA correction has been computed separately (Sect. 2.5). We compute the four GIA cor-rections described in Sect. 2.5 (Table 4), but we show its relative contribution only for the average of the model Riva09 and IJ05-LV (green bars, Fig. 8a). Each of the four GIA corrections contributes to the trend with quite different proportions, and this becomes clear when observing the GIA variability (Max-Min) (green bars, Fig. 8b). So the average of all the corrections is not very meaningful nor representative, and in some cases, when the maximum and minimum value among the four GIA corrections have opposite sign, the average is zero. The Riva09 and IJ05-LV in the typical GIA basin (1, 2, 18) give similar corrections, and in some large basins in East Antarctica (17, 2, 3, 10, for example) they both show much smaller correction than both the ICE5g (VM2 compressible and LV). Using the average of the model Riva09 and IJ05-LV is just one reasonable choice among other reasonable choices.
The GIA variability (Fig. 8b) is largest in most of the largest basins. However the most typical GIA pattern affects mainly basins 1, 2, 18 and 19, which are among those  with largest GIA trend. In several basins, the GIA contribution (green bars, Fig. 8a) is significant compared to the data trend (blue bars), and in basins 17, 15, 2, 19, 11 and 24 it is dominant. So the GIA contribution has a relative importance also in basins not within the typical GIA regions but which have small trends in the data (blue bar, Fig. 8c), e.g. as in basins 17, 15, 11 and 24. The GIA contribution for the whole of Antarctica that must be subtracted to the mass trend is +71 Gt yr −1 . The degree-1 contribution is not dominant in any of the basins (sky-blue bars, Fig. 8a), and it only exceeds 10 % of the sum of contributions in those basins associated with a low mass trend (blue bar, Fig. 8c), e.g. in basins 17, 10, 11 and 16. The degree-1 contribution for the whole of Antarctica that must be added to the mass trend is +13.2 Gt yr −1 .
The GAC correction (orange bars, Fig. 8a) is more important than the degree-1 correction for almost all basins, and in some cases it has the same impact as the data trend (blue bars). The GAC correction for the whole of Antarctica that must be added to the mass trend is −20 Gt yr −1 .
The accuracy errors on the trends (which account for difference between data set and methods, and so also account for leakage errors) have about the same impact as the error on data trend computation (blue and light-blue bars, Fig. 8d). It is interesting to notice that where the two largest data trends occur (blue bars, Fig. 8c), in basin 21 and 22 (Amundsen sector), we also find the largest (absolute) accuracy errors. However basin 20, the third largest data trend, has one of the smallest accuracy error. As noticed above, the relative errors for these basins are the smallest.
The degree-1 uncertainties on the trend, which are primarily caused by the GIA degree-1 uncertainties (sky-blue bars, Fig. 8d), are comparable with the GIA uncertainties (green bars) for many basins. And for both these contributions, the uncertainties are larger for larger basins.
The uncertainties on trend computation for the GAC correction are very small for all the basins (orange bars, Fig. 8d) For Greenland the trends computed only on data are much more important than the correction contributions (Fig. 9a), and the errors are below 20 % for all basin but one, the largest. The largest basin (number 2), in northwest Greenland, has a small trend and small GIA correction which have similar contributions in the total trend. Also the GAC correction has the largest value in northwest Greenland. The GAC correction for the whole of Greenland is −6.3 Gt yr −1 .
The largest GIA variability is located in basin 6, in southwest Greenland (green bars Fig. 9b), where also the largest total trend (blue bars, Fig. 9c) is found, and yet the trend on data in basin 6 dominates with almost 85 % of the total trend.  Fig. 9. Trend summary for Greenland with same meaning as in Fig. 8. The GIA contribution that we used in this case (i5g-CP) that must be subtracted to the mass trend is −5.3 Gt yr −1 . The degree-1 correction has small impact on the total trend (sky-blue bars, Fig. 9a) but contributes significantly the total error on trend (sky-blue bars, Fig. 9d). The degree-1 contribution for the whole of Greenland is −3.5 Gt yr −1 .

Accelerations in mass changes
As the last result, we present the trends on four different periods (Table. 6) and related accelerations for our derived time series using the same configuration as in Figs. 8 and 9. The total mass balance for the whole GRACE period (2003-2011) is found to be −234 Gt yr −1 for Greenland and −83 Gt yr −1 for Antarctica where most of the mass loss is going on in West Antarctica with −111 Gt yr −1 (errors are in Table 6). A more rapid mass loss clearly takes place in Greenland and West Antarctica (Table 6)  This increase or acceleration is even clearer at basin scale ( Figs. 10 and 11). For Antarctica, however, the acceleration in mass loss in most of the basins in the western part (Fig. 10) is counteracted by an increase of accumulation in the eastern part. In western Greenland (Fig. 11) the mass loss increased in the last 5 yr, while in the eastern part the mass loss has decreased. Accelerations are computed as differences between trends obtained over two time periods with respect to our robust error estimate. We analyze two kinds of acceleration, one (blue bars in Figs. 10 and 11) related to the last time frame of GRACE data (2007)(2008)(2009)(2010)(2011) against the full one (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011). And the other acceleration (green bars) is related to the last time frame of GRACE data (2007)(2008)(2009)(2010)(2011) against the first one (2003)(2004)(2005)(2006)(2007). When both accelerations are larger than the error on trends, we have reliable information about this quantity.

Conclusions
In light of the consistent and systematic error analysis that we have performed, the results presented about time series and mass acceleration in this study are statistically meaningful, and yet trends on mass balances especially in Antarctica can still be systematically biased by GIA correction.
For the first time the various sources of variability in mass change estimates have been altogether systematically assessed, and we quantify their associated uncertainties for the mass trends as well as the monthly mass change solutions. We find monthly differences between the results obtained with use of RL05 with respect to the use of RL04 comparable with our derived total monthly error. The correlation coefficients between the two release being below 75 % for several basins in Antarctica confirm a difference in noise content in the new release 05.
We cross-validate our two independent methods, and the clear agreement between the two confirms that the low resolution of the input GRACE data allows us to use very simple leakage treatment like the one employed in the conversion method. A surprisingly large part of the variability in the monthly solution arises from the use of different data sets rather than different methods. However, the uncertain--4.00 -3.00 -2.00 -1.00 0.00 1.00 2.00 3.00 ties on degree-1 also largely contribute to the variability in the monthly solutions.
Moreover the performance of the inversion method is also tested in Shepherd et al. (2012), which, together with our validation, makes us confident that we have addressed most of the methodological uncertainties.
The choice of degree-1 is still an open issue; hence we build a degree-1 sensitivity kernel (at basin scale) which represents a solid and lasting tool to perform straightforward computation of the degree-1 correction on mass balance using any geocenter motion time series. For our preferred estimate, we then choose to use an average of three available geocenter motion time series based on different methods.
The trend in the geocenter motion correction is found to be +13.2 Gt yr −1 for Antarctica and −3.5 Gt yr −1 for Greenland. By considering that our preferred GIA correction is about +70 ± 11 Gt yr −1 for Antarctica and +2 ± 7 Gt yr −1 for Greenland, the degree-1 correction represents a significant contribution to both the trend and uncertainties. We find that the scatter in the monthly solutions caused by applying different estimates of geocenter motion time series (degree-1 corrections) is significant -contributing with up to 40 % of the total error.
We show the impact of the correction for the de-aliasing GRACE product (GAC), especially on the trends in Antarctica mass balance, and we generate an alternative correction that can be applied to the whole RL04 time series, while we wait for the RL05 to be completely released and tested. The GAC correction allows us to quantify an www.the-cryosphere.net/7/1411/2013/ The Cryosphere, 7, 1411-1432, 2013 error affecting previous mass balance estimates, which corresponds to −20 Gt yr −1 for Antarctica and −6.3 Gt yr −1 for Greenland. The outcome of this systematic analysis is a set of our preferred monthly solutions and their associated error estimate, which is a combination of precision error (propagated from the data) and accuracy error due to the method and the different data set. Furthermore, we provide our preferred degree-1, GIA and GAC corrections for both the monthly solution and the trends. We provide trends for each basin with their associated total error estimate composed by precision and accuracy error added to uncertainties due to degree-1 and GIA choice, and GAC correction. The computed total error is more than double the simple precision error on the trend, and in large basins the degree-1 uncertainty is as important as the GIA one.
Since trends often depend on the choice of the time interval, we compute trends over the whole period 2003-2011 and sub-periods 2003-2006 and 2007-2011. We find a clear increase in ice loss in the sub-interval 2007-2011 only for West Antarctica and Greenland.