Earth and Planetary Science Letters Revisiting GRACE Antarctic ice mass trends and accelerations considering autocorrelation ✩

Previous GRACE-derived ice mass trends and accelerations have almost entirely been based on an assumption that the residuals to a regression model (including also semi-annual, annual and tidal aliasing terms) are not serially correlated. We consider ice mass change time series for Antarctica and show that signiﬁcant autocorrelation is, in fact, present. We examine power-law and autoregressive models and compare them to those that assume white (uncorrelated) noise. The data do not let us separate autoregressive and power-law models but both indicate that white noise uncertainties need to be scaled up by a factor of up to 4 for accelerations and 6 for linear rates, depending on length of observations and location. For the whole of Antarctica, East Antarctica and West Antarctica the scale factors are 1.5, 1.5 and 2.2 respectively for the trends and, for the accelerations, 1.5, 1.5 and 2.1. Substantially lower scale-factors are required for offshore time series, suggesting much of the time-correlation is related to continental mass changes. Despite the higher uncertainties, we ﬁnd signiﬁcant (2-sigma) accelerations over much of West Antarctica (overall increasing mass loss) and Dronning Maud Land (increasing mass gain) as well as a marginally signiﬁcant acceleration for the ice sheet as a whole (increasing mass loss). © 2013 The Authors. Published by Elsevier B.V. All rights reserved. temporal in time series estimate ice mass and not GRACE

In general, trends and acceleration in Antarctic time series have been computed from a simple linear regression under the assumption that the time signatures are deterministic, typically with the addition of semi-annual, annual and terms relating to tidal alias-✩ This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial-No Derivative Works License, which permits noncommercial use, distribution, and reproduction in any medium, provided the original author and source are credited.
ing. In addition, when fitting this regression (functional) model to the data it has been assumed that the residuals to the observations contain no serial correlation and can be accounted for within the stochastic model of ordinary least squares. That is, the reported uncertainties were computed assuming the observations are temporally uncorrelated, and the residuals may be modeled as white (normally distributed) noise, with a constant power spectral density. However, a white noise assumption has been shown to result in overly optimistic uncertainty estimates for other geodetic data sets (Hughes and Williams, 2010;Williams, 2003). Small changes in uncertainty may be particularly important for the Antarctic mass acceleration term; the values computed by Velicogna (2009) (−26 ± 14 Gt/yr 2 ), King et al. (2012) (−3.6 ± 5.6 Gt/yr 2 after scaling their uncertainties to 1-sigma) and Velicogna and Wahr (2013) (−12 ± 9 Gt/yr 2 ) are either within or close to the noise limit, and all are less than 2-sigma.
In regression analysis one may argue that the noise of the system is white and can be fully described by the formal errors of the data; all other variations are true changes in continental mass. The uncertainties of any estimated parameters of a functional model should simply reflect the system noise. However, the functional model may not be able to capture all of the true changes and, therefore, what is remaining may alter the estimation of the parameters and our measure of the uncertainty of those parameters. Blewitt (1998)  stochastic models. We can therefore determine a linear rate with an appropriate stochastic model or we can augment the functional model with additional terms and simplify the stochastic model and both are equivalent; meaning that any unmodeled signal/noise can be accommodated through inclusion within either the functional or stochastic models. That is, one person's noise is another's signal. Importantly, all signal/noise must be appropriately modeled, one way or the other, for the estimated parameters and their uncertainties to be robust.
Here we examine the stochastic properties of GRACE data and show that serial correlation is present in residuals to the commonly applied model. In particular, we investigate the spatial distribution of the uncertainties on GRACE-derived ice mass changes to identify areas where the ice mass trends and accelerations are statistically significant. While following most other authors in fitting such a deterministic model, we do so aware that time series non-linearity, or stochastic interannual variations (e.g. Sasgen et al., 2010), could be used to argue that there are no deterministic trends or accelerations just stochastic variations; the separation of trends into deterministic and stochastic is problematic and presents a major challenge (Fatichi et al., 2009). However, mass change rates and accelerations as simple metrics are the basis of many studies including the fourth assessment report of the Intergovernmental Panel on Climate Change (IPCC) (Solomon, 2007) and the pending fifth report. There is thus a need to understand and to quantify the observed metrics while taking into account stochastic variations. Here we use a functional model, which may include linear and quadratic terms, to describe the long term state over the measurement period particularly where the stochastic model cannot explain the changes. We limit the linear and quadratic terms to the data period only and do not assume that they have any forward or backward predictability properties. Our estimated trends simply indicate that there has been a general increase (or decrease) in mass while accelerations indicate a general increase (or decrease) in that trend over the measurement period.

Data
We consider the Center for Space Research (CSR) Release 5 (RL05) monthly GRACE gravity fields  spanning 108 months between March 2003 and July 2012. This study is a time extension of the April 2002-May 2010 study of Velicogna (2009) and the August 2002-December 2010 study of King et al. (2012) and is comparable to the January 2003-November 2012 study of Velicogna and Wahr (2013). Perhaps more pertinently, our study and that of Velicogna and Wahr (2013) utilized the latest Release 5 (RL05) while the other studies used Release 4 (RL04). RL05 is a reanalysis of the GRACE data incorporating improvements in modeling and data quality and gives significant improvement in accuracy and spatial resolution over the previous version (Bettadpur, 2012).
The monthly GRACE gravity fields consisting of spherical harmonics to degree and order 60 were pre-processed as in the main text and supplementary material of King et al. (2012), including the approach for adding degree-1 harmonics (Swenson et al., 2008) and replacing the degree two harmonic C 20 with the value from satellite laser ranging (Cheng and Tapley, 2004). All spherical harmonics of order 8 or more were filtered  to reduce the correlated errors in the GRACE gravity field coefficients. A quadratic polynomial in a moving window of width 4 centered on the degree n was employed to filter the Stokes' coefficient of order m. It is noted that the destriping method may affect the estimated mass rates (Velicogna and Wahr, 2013).
Regional studies of mass redistribution from GRACE suffer from leakage effects from both within and outside the region of study.
To reduce leakage from hydrological and oceanic geophysical signals outside Antarctica we utilized monthly values of total water storage (TWS) from the Global Land Data Assimilation System model (GLDAS) (Rodell et al., 2004). To conserve global water mass a uniform layer was applied to the ocean mass for each month. This was derived by consideration of the total ocean mass variation (the C 00 harmonic) as supplied by the GRACE dealiasing product and the total water storage from GLDAS. It is also known that the ocean model underlying the GAD fields contains artificial trends (http://www-app2.gfz-potsdam.de/pb1/op/grace/aod_issues/ issues_aod1b_rl05.html). Analysis at GFZ has shown that the trend and low frequency variability in the ocean model are related to warming and cooling of water masses at intermediate depths. The trend is considerably less reliable than the high-frequency. Accordingly, the trend was extracted from the GAD product and restored to the GRACE gravity fields. This reduces the resultant trend within the surrounding seas, with particular impact within the Ross Sea and Ice Shelf.
The final correction is for Glacial Isostatic Adjustment (GIA). Application of a GIA model is required to recover meaningful estimates of ice-mass change over Antarctica. We use the W12a GIA model  which compares well against bedrock vertical rates from GPS observations (Thomas et al., 2011). However, despite this recent progress the GIA correction in Antarctica is still very uncertain but since GIA may be considered linear over the study period it does not affect the monthly data weights or the non-linear aspects of the time series. We subtracted the modeled GIA spherical harmonics from the GRACE fields and derive mass variation (under the assumption that all change is concentrated at the surface) by converting to surface spherical harmonics incorporating the Earth's elastic response through load Love numbers (Wahr et al., 1998).
To recover mass change over Antarctica it is necessary to reduce the noise in the higher degree and order GRACE harmonics. Two techniques were considered, namely spatial averaging and the kernel function approach. The former yields the spatial variation of ice mass change and is given here in terms of mm of equivalent water height. We produced GRACE time series for locations on a regular 500 km grid across Antarctica using the Gaussian spatial averaging approach of Wahr et al. (1998) with an averaging radius of 400 km. Since the size of the grid is at the limits of the resolution of the GRACE data, signal may leak from one grid point to another. We did not attempt to correct for this problem, nor is this effect included in our uncertainty estimates. In the second technique, areal averages of ice mass change (in Gt/yr) for the Antarctic Ice Sheet (AIS), the West Antarctic Ice Sheet (WAIS) and the East Antarctic Ice Sheet (EAIS) were derived using averaging kernels and their scale factors as given in the supplementary material of King et al. (2012). The kernels are derived by defining a function with values of unity over the required area (zero otherwise) and representing the function in terms of spherical harmonics. Since the kernel function is restricted to the same maximum degree and order 60 as the GRACE fields the representation is imperfect, particularly at the coast, where the kernel seeks to follow a step change from unity to zero. In mitigation, the area of the kernel was extended outwards from the coast of Antarctica by 150 km. Thus, values closer to unity are achieved at the coast where much of the mass change is occurring. The scale factor is applied to ensure that 1 cm of ice mass over the kernel representation is equivalent to the same over Antarctica. No further averaging is necessary as the convolution acts as a filter. We assign uncertainties to each monthly solution based on the approach of Wahr et al. (2006), although these uncertainties are only used in our analysis for variable white noise; otherwise, our derived uncertainties come from the time series themselves. We did not make the small correction related to ocean mass redistribution following ice mass changes (Sterenborg et al., 2013).

Data analysis and results
We initially began by testing pairs of functional and stochastic models using functional models that differed only in the estimation, or not, of a quadratic term, and 8 choices of stochastic model. The tested stochastic models include white noise only, time-variable white noise (using scaled individual formal errors so that the amplitude of the white noise can change on an epoch to epoch basis), autoregressive (AR) noise (order 1 only), powerlaw noise (Press, 1978) and combinations of the above. Higher order AR models were originally investigated but were dismissed as they did not achieve the criteria described below. We used Maximum Likelihood Estimation to simultaneously solve for the time series parameters (intercept, trend, quadratic when estimated, annual and semi-annual signals and the 161 day aliasing period in GRACE from S2 semidiurnal solar tide (Chen et al., 2009)) and the stochastic noise parameters (Langbein and Johnson, 1997;Williams et al., 2004). For the power-law noise we also solved for the spectral index, namely the slope of the power spectral density in log-log space, where −2 indicates random walk and 0 indicates white noise. Finally, for each stochastic/functional model pair we calculated the Bayesian Information Criterion (BIC) (Schwarz, 1978). The functional/stochastic model pair with the lowest BIC was selected as the preferred solution at each grid cell.
When testing a functional model including an acceleration term we solved for the following quadratic equation for consistency with previous work Velicogna, 2009;Velicogna and Wahr, 2013): We refer to the term c as the acceleration term. For each grid cell we adjust t 0 such that the trend, b is equal to the trend in the linear fit so that the plot of trends is independent of whether we have fitted an acceleration at that point. This does not affect the value of the acceleration term. Realistic parameter uncertainties were propagated from the optimum stochastic model for each time series. For each stochastic model and time series we also calculate a scale factor, which is the amount we have to scale the parameter uncertainty derived using a white noise assumption such that it equals the uncertainty derived from the optimum stochastic model.
The calculated linear rates and accelerations are shown in Fig. 1 based on the optimal stochastic model for each grid cell. The acceleration is shown regardless of the preferred functional model for sake of illustration and we discuss the preferred functional/stochastic model pairs for a series of points below. Also shown are the 2-sigma contour lines, that is, the contour lines representing the areas where the accelerations are significant at the 2-sigma level (black lines). Alternatively, the green boxes indicate the boundary of the grid cells where the BIC indicates that the quadratic model is a better fit to the data than the linear model, i.e. that the optimal functional model includes an acceleration. We note that the two outlines are in close agreement. Fig. 1b shows the estimated linear changes in mass. The grey areas show the regions where the trend is not significant at the 2 sigma level. Almost all of the area has a significant trend except for a couple of isolated patches in Antarctica and some larger regions in the southern oceans. The six points marked in Fig. 1 are representative of different regions and their time series are shown in Fig. 2. Point A, in the Weddell Sea, shows an insignificant rate with a preferred stochastic model of white noise. The BIC at Point B, in Graham Land within the Antarctic Peninsula, indicates that the preferred functional model is linear with a preferred stochastic model being a combination of AR and white noise, which indicates a significant time-correlation. At Point C, in the region of the Pine Island Glacier and within an area of significant negative accelerations (increasing mass loss), the preferred stochastic model is AR.
There we find an acceleration of −17.0 ± 3.3 mm/yr 2 (all quoted uncertainties are, from here-in, 2-sigma). The uncertainty using AR is 2.2 times larger than that when using a white noise estimate. Point D, along the coast of Dronning Maud Land, is in a region of significant positive acceleration (increasing mass gain). The preferred stochastic model is AR plus time-variable white noise.
Point E lies just outside a region of significant positive acceleration. Under the usual white noise assumption the acceleration in this region is 11 times greater than its estimated uncertainty and would be interpreted as highly significant. However, with the preferred stochastic model, AR plus white noise with an autocorrelation parameter of 0.96 ± 0.03, the acceleration term is no longer significant. This area was subject to snowfall-driven mass change that began in 2009 (e.g. Boening et al., 2012), as can clearly be seen in Fig. 2, and therefore a continuous acceleration over the whole period is not the most appropriate model. Point F, inland of Totten Glacier in Wilkes Land, also has no significant acceleration but very obvious time variations that manifest themselves clearly as an AR plus variable white noise model with an autocorrelation coefficient of 0.93 ± 0.08. Also plotted on the graph is the anomaly to the modeled cumulative surface mass balance (SMB) (relative to the 1979-2010 mean) for the drainage basin containing this point (Basin 13; Zwally et al., 2012) derived from RACMO2_ANT27 . The close agreement between the GRACE and SMB anomaly time series for years 2008-2010 is noteworthy and indicates that the changes in these areas are driven by surface processes and not ice dynamics.
The preferred stochastic model with respect to the preferred functional model is shown as a histogram in Fig. 3 with points partitioned based on their location on land (including the immediate 200 km offshore) or in the ocean. There is no obvious dominant stochastic model but we can make some general observations. Grid cells for which the white noise only or the time-variable white noise models are chosen, based on the BIC, are mainly in the ocean. AR(1) noise, with or without the addition of white noise or time-variable white noise, is the most dominant model but power-law noise is also close and has the largest individual share for the land points where the preferred functional model is linear (although this share is smaller than when the AR(1) models are counted together).
If we were to recommend a single model to use in studies such as this then we would choose AR(1) noise at present but note that this could change as the length of the time series increases. The scaling factor for the acceleration uncertainties for the optimal stochastic model compared with the white noise derived uncertainties is shown in Fig. 4a. We see that accounting for serial correlation in the time series leads to an increase in the acceleration uncertainty up to a factor of 4. The median scaling factor for Antarctica (including the immediate 200 km offshore) is 2 whereas in the surrounding oceans the median is 1.1. This appears to indicate that the origin of the stochastic signals is ice mass change and not noise in the GRACE system. A similar pattern (Fig. 4b) in the scaling factor for rate uncertainties is also found except the maximum increases to 6 onshore. The scaling factors are a good proxy for the level of temporal correlation in the time series; an area where the scaling factor is close to 1 indicates there is little correlation in the time series (i.e., it is white) whereas areas where the scaling factor is high indicate significant correlation. Fig. 4a are the boundaries where the BIC indicates that the preferred functional model is quadratic. We note that the areas of high temporal correlation, as implied by a high scaling factor, are mainly located along the Antarctic coastline but are not particularly correlated with the areas of significant acceleration. Indeed the largest scaling factor occurs outside all areas of significant acceleration (near Point E), and is located where Boening et al. (2012) reported significant snowfall-driven mass change. The stochastic model accounts for these variations that are not well-fit by a linear or quadratic functional model. A high scaling factor is seen near Point F for a similar reason. Mis-modeled near-coastal ocean processes may also contribute to the temporal correlations evident at locations along the coast.

Superimposed on
The preferred functional models, based on the BIC, using the kernel function approach for the entire AIS, WAIS and EAIS are shown in Fig. 5. For AIS, and for March 2003 to July 2012, the preferred functional model is a linear trend of −58 ± 16 Gt/yr, a marginally significant acceleration of −15 ± 13 Gt/yr 2 , and an AR stochastic model with a coefficient of 0.45 ± 0.17. The amount we need to scale the white noise only solution uncertainties to match those from the preferred stochastic model is 1.5 for both the trend and acceleration. Both EAIS and WAIS, when treated independently, show significant accelerations of +18.2 ± 10.2 Gt/yr 2 and −31.3 ± 6.9 Gt/yr 2 , respectively, with trends of +97 ± 13 and −159 ± 9 Gt/yr, respectively. The preferred stochastic model for WAIS and EAIS is AR with coefficients of 0.74 ± 0.13 and 0.45 ± 0.17, respectively, and with trend (acceleration) uncertainty scaling parameters compared to the white noise only solution of 2.2 (2.1) and (1.5) (1.5), respectively. Comparing RL04 with RL05 over the common period (Mar. 2003to Dec. 2010) showed rates to be higher in RL05 (but not significantly) by 15.5 ± 23.1 Gt/yr, 12.1 ± 18.5 Gt/yr, and 3.2 ± 4.2 Gt/yr for AIS, EAIS and WAIS, respectively. Differences in accelerations were also not significant with a change of −4.5 ± 14.8 Gt/yr 2 for EAIS, −1.4 ± 3.7 Gt/yr 2 for WAIS and −3.2 ± 18.5 Gt/yr 2 for AIS.

Deterministic Antarctic accelerations
In the above we have identified regions where the acceleration is significant in the sense that there is an acceleration-like signal that exceeds that due to stochastic variation as quantified by the error bar. The error bar itself describes the range of possible accelerations due to the stochastic model given in the covariance matrix. However, it is possible that apparent deterministic accelerations are just stochastic variations (Fatichi et al., 2009) and that fitting acceleration terms removes long period trends from the signal, biasing our estimation of the stochastic properties towards whiteness. To investigate the spatial distribution of deterministic accelerations let us consider the following. If we assume there is no deterministic acceleration we could look at the stochastic noise model due to the linear model only and use those to predict the range of accelerations by propagating the covariance matrix. The regions where we can ascribe deterministic accelerations corre-spond to the regions where the predicted "deterministic" acceleration is still too large for the stochastic variation to encompass. In Fig. 6 the black contours are those from Fig. 1a with the red contours the new 2-sigma areas using the linear-only stochastic parameters. The area with significant accelerations has shrunk but the regions of large acceleration persist and we conclude they are deterministic.
In one further test we assumed that the time series were purely stochastic in nature and estimated the best fitting stochastic model. Power-law noise was the optimal model for over 85% of points, increasing to 92% for those points on land. The average spectral index (see Fig. 7) for the land-only points is −1.5, with a range of between −1.0 and −2.3 indicating that the time series are non-stationary in behavior. For the points in the ocean the spectral index is closer to 0 (white) with an average of −0.8 and range of −0.1 and −1.8. In terms of the BIC, only one quarter of points have a lower value for the purely stochastic model compared to the linear or quadratic models, and these points are not situated where the accelerations and trends are largest. In this case, because we are testing whether a stochastic model can reproduce both a trend and an acceleration of a certain magnitude, as opposed to either the trend or acceleration separately, we examine the combined two-dimensional confidence ellipse to test for significance. Fig. 8 shows the confidence we can assign to each grid cell that the estimated trend and acceleration cannot be generated by chance from the pure stochastic model. There are large areas over the continent where the significance is still greater than 80%. Point C for instance has a confidence of 94% that the estimated acceleration and trend combined would not be generated by chance from the purely stochastic model. So while we cannot rule out the purely stochastic results we are equally confident that we can describe the results with deterministic trends and accelerations together with the realistic uncertainties described here. If we were to use the pure stochastic model to assign uncertainties then the scaling from the white noise assumption is up to 10-15 times larger for the acceleration, and up to 20 times for the trend in the region of Pine Island Glacier. While confidence levels above 80% are shown for large parts of the interior of East Antarctica in Fig. 8, this is a result of very low noise in this region (likely reflecting the low accumulation rates) rather than large magnitude trends or acceleration (see also Fig. 1).

Discussion
Discussion of GRACE errors to date has focused on the dominant systematic errors related to correcting for mass change due to glacial isostatic adjustment (GIA) (e.g., Velicogna and Wahr, 2006b); or due to different raw data analysis strategies (e.g., Sasgen et al., 2007). Other important discussions have surrounded methods for obtaining reliable uncertainties for a given GRACE monthly solution (e.g., Horwath and Dietrich, 2006;Wahr et al., 2006) or reducing spatially-correlated errors (e.g., Sasgen et al., 2006;Swenson and Wahr, 2006). We are not aware of detailed discussion of temporal correlations in GRACE mass time series apart from a very brief comment by Horwath and Dietrich (2009) who state that a factor √ 2 scaling of white noise uncertainties is required. Similarly, Velicogna and Wahr (2013) acknowledge that the actual ice mass loss is not perfectly represented by the regression and include a contribution to the stochastic uncertainty to account for autocorrelation. Wahr et al. (2006) also briefly discuss temporal correlations in the raw GRACE K-band data. Recently Wouters et al. (2013) examined temporal correlations in GRACE time series using long-term mass balance time series to calculate the effect on robustness of trends and accelerations. Our calculated uncertainty for the acceleration of the AIS of ±13 Gt/yr 2 is equal to their estimate of ±13 Gt/yr 2 and this reinforces our conclusion that the origin of the stochastic signals is mainly ice mass change and not noise in the GRACE system. Our study considers both temporal correlations in the GRACE mass time series as well as detection of deterministic or stochastic trends in the time series. Our results show that temporal correlations in GRACE cannot be ignored (either as part of the stochastic model or the functional model) and suggest that previously published GRACE mass rate and acceleration uncertainties are un- derestimated. In addition, by studying time series at individual locations we concluded that the trends are deterministic. Of particular significance are areas of Antarctica where the acceleration is estimable and deterministic. By combining the preferred stochastic model for each location with the linear model we identified regions where the acceleration can be shown to be deterministic and hence indicate areas of increasing ice mass change over the measurement period. These areas of acceleration, although smaller than the regions identified when using the quadratic functional model, are still significant.
While our findings pertain directly to Antarctic time series, it is true that considering the presence of time-correlated noise can only inflate parameter uncertainties generated under a white noise assumption. It is highly likely that similar mis-modeled signals exist in other locations. Consequently, published GRACE parameter (e.g., trend, periodic and acceleration) uncertainties for other regions must also be regarded as a lower bound.
Realistic uncertainties are particularly important for studies which seek to optimally combine GRACE time series or rates with other data sets (e.g., Hill et al., 2010;Riva et al., 2009;Wu et al., 2010). In these cases, the weighting assigned to GRACE is critical and we conclude that GRACE data have been incorrectly weighted in previous solutions of this type. However, analysis of the other data sets may equally underestimate uncertainties due to temporal correlations (although not always, see for example Ferguson et al., 2004) and, if the correlation comes from the ice mass signal, as is likely with altimetry data for example, then these data sets would suffer similar underestimation of trend and acceleration uncertainties.
Other studies have used GRACE rates to test or constrain models of GIA (e.g., Steffen et al., 2010). In these cases, statistical tests based on GRACE uncertainties may have been used to reject models which would not have been rejected with more rigorous uncertainties.
In terms of our estimates of ice mass change, our finding of statistically significant accelerated mass loss being limited entirely to the region draining into Pine Island Bay is in agreement with satellite altimeter studies that have identified accelerations in this region (Flament and Rémy, 2012;Wingham et al., 2009). , using a white noise scaling factor of 2, found that statistically significant ice loss acceleration was limited to just one drainage basin in this region; Fig. 4 suggests this scale factor may have been too small for this region and too large for much of the rest of the continent.
Our finding that AIS mass changes support a marginally significant acceleration is at odds with the findings of Velicogna (2009) and King et al. (2012) who, based on shorter RL04 time series, both reported accelerations of mass loss that were insignificant at 2-sigma (−26 ± 28 Gt/yr 2 (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009), errors converted to 2-sigma, and −4.4 ± 16.0 Gt/yr 2 (2002-2010), respectively). However, the range of estimates and their uncertainties are such that we can also say that they are not significantly different from each other either. Whilst we extend their findings to show that an acceleration term is, at this stage, representative of the AIS time series we believe that another five years is probably required before we are truly confident of this. King et al. (2012) also report separate acceleration terms for EAIS (+4.1 ± 11.7 Gt/yr 2 ) and WAIS (−8.4 ± 8.4 Gt/yr 2 ) and these are substantially different to those reported here, presumably not due to differences between GRACE RL04 and RL05, but different analysis approaches and time series spans and the general brevity of the time series. Finally, using RL05 Velicogna and Wahr (2013) report changes of −83 ± 49 and −147 ± 80 Gt/yr for two different GIA models (both different to that which we use here), larger magnitude rates than we report here (−58 ± 16 Gt/yr) over a similar time period, with an acceleration of −12 ± 18 Gt/yr 2 (converted to 2 sigma) that is in agreement with our estimate and dominated by the southeast pacific sector of West Antarctica and the Antarctic Peninsula.

Conclusions
We examined GRACE mass change time series for Antarctica and found serial correlation in regression residuals. We found that an autoregressive model fits the residuals sufficiently in most places but we could not rule out the power-law model. When applying an autoregressive model we found that GRACE trend and acceleration uncertainty estimates for Antarctica have previously been underestimated by a factor of between 1 and 4 (6 for trends) when using monthly uncertainties computed according to . We note, however, that other errors exist which we have not considered, including those related to raw data analysis strategy and areal averaging effects (see Horwath and Dietrich, 2009;Velicogna and Wahr, 2013). Systematic bounds in GIA model uncertainty must also be considered Velicogna and Wahr, 2013;Whitehouse et al., 2012).
While we have focused on Antarctic mass changes here, these results will apply to other regions, as serial correlations are likely to be present everywhere. Since the adopted noise models are apparently driven by geophysics rather than GRACE analysis noise, the preferred noise model may change as GRACE time series extend further in length, meaning that it is especially important that all analysts consider the appropriate functional and stochastic model for their own time series.
Over March 2003 to July 2012, East and West Antarctica ice mass change was +97 ± 13 and −159 ± 9 Gt/yr, respectively, with accelerations +18 ± 10 and −31 ± 7 Gt/yr 2 , respectively (2-sigma uncertainties) not considering GIA model error bounds. Mass change for the entire Antarctic Ice Sheet is also best modeled with a linear plus acceleration functional model and a stochastic model that considers temporal correlations, giving an ice mass trend of −58 ± 16 Gt/yr and an acceleration of −15 ± 13 Gt/yr 2 .
Over this time period contribution to global-mean sea-level change was, therefore, +0.16 ± 0.04 mm/yr with an acceleration of 0.04 ± 0.03 mm/yr 2 .