Improving the estimate of the secular variation of Greenland ice mass in the recent decades by incorporating a stochastic process

Abstract The irregular interannual variations observed in the Greenland ice sheet (GrIS) mass balance can be interpreted as stochastic. These variations often have large amplitudes, and, if not accounted for correctly in the mass change model parameterization, could have profound impacts on the estimate of the secular trend and acceleration. Here we propose a new mass trajectory model that includes both the conventional deterministic components and a stochastic component. This new model simultaneously estimates the secular rate and acceleration, seasonal components, and the stochastic component of mass changes. Simulations show that this new model improves estimates of model parameters, especially accelerations, over the conventional model without stochastic component. Using this new model, we estimate an acceleration of −1.6 ± 1.3 Gt/yr2 in mass change (minus means mass loss) for 2003-2017 using the Gravity Recovery and Climate Experiment (GRACE) data and an acceleration of −1.1 ± 1.3 Gt/yr2 using the modeled surface mass balance plus observed ice discharge. The corresponding rates are estimated to be −288.2 ± 12.7 Gt/yr and −274.9 ± 13.0 Gt/yr. The greatest discrepancies between the new and the conventional model parameter determinations are found in the acceleration estimates, −1.6 Gt/yr2 vs. −7.5 Gt/yr2 from the GRACE data. The estimated accelerations using the new method are apparently smaller than those estimated by other studies in terms of mass loss. Our quantitative analysis elucidates that the acceleration estimate using the conventional method is the lower bound (i.e., −7.5 Gt/yr2 for 2003–2017) while the acceleration estimated by the new method lies in the middle of the possible ranges. It is also found that these discrepancies between the new and the conventional methods diminish with sufficiently long (>20 yr) observation records.


Introduction
The Greenland ice sheet (GrIS) loses its ice mass through runoff of meltwater and ice discharge processes from the ice sheet to the ocean. Both processes have increased over the past two decades and have resulted in accelerated ice mass loss (Hanna et al., 2013;Velicogna et al., 2014;Khan et al., 2015;Seo et al., 2015;King et al., 2018;Aschwanden et al., 2019;Bevis et al., 2019). Ice loss from the GrIS has become one of the largest contributors to sea level rise (Cazenave and Remy, 2011;Shepherd et al., 2012;Jacob et al., 2012;Chen et al., 2013;Andersen et al., 2015). This acceleration in GrIS ice mass loss was also accompanied by strong interannual * Corresponding author.
E-mail addresses: sggzb@whu.edu.cn (B. Zhang), liulin@cuhk.edu.hk (L. Liu), ybyao@whu.edu.cn (Y. Yao), tonie.vandam@uni.lu (T. van Dam), abbas@space.dtu.dk (S. Abbas Khan). and seasonal variations (Velicogna et al., 2005;Chen et al., 2011;King et al., 2018;Zhang et al., 2019). The seasonal variations are much more regular than the interannual variations so that they are relatively straightforward to quantify and separate from the other components in the time series. However, the acceleration and the interannual variations that are associated with climate variability (Wouters et al., 2013;Seo et al., 2015;Bevis et al., 2019;Zhang et al., 2019) are not as predictable as the seasonal variations. In addition, the interannual variation can behave like an acceleration over certain time periods, which makes it more difficult to quantitatively separate the acceleration from the interannual variations. As a result, the acceleration will be incorrectly determined if the potential aliasing between the interannual variation and the acceleration is not given adequate consideration.
Greenland mass change is usually characterized by a mass trajectory model consisting of a constant, a rate, an acceleration, and some periodic terms (e.g., Svendsen et al., 2013;Wouters et al., https://doi.org/10.1016/j.epsl.2020.116518 0012-821X/© 2020 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 2013; Bevis et al., 2019). Some studies used a simpler model consisting of only a constant, a rate and a quadratic term to describe Greenland mass changes (e.g., Velicogna, 2009;Seo et al., 2015). Many other studies have also used one of the above two models to estimate the rate and acceleration of Greenland mass loss (e.g., Rignot et al., 2011;Schrama and Wouters, 2011;Velicogna et al., 2014). However, significant differences exist between published estimates of these rates and accelerations estimated from even the same Gravity Recovery and Climate Experiment (GRACE) data. Despite the differences in the time spans that were analyzed, a significant part of the discrepancies can be attributed to an imperfect understanding of the nature of Greenland mass change.
This is particularly evident in acceleration estimates of the GrIS mass loss. Based on GRACE data, Velicogna (2009) Rignot et al. (2011) also reported an acceleration of −21.9 ± 1 Gt/yr 2 for 1992-2009 based on the mass budget approach (modeled surface mass balance minus observed ice discharge). Wouters et al. (2013) reported an acceleration of −25 ± 9 Gt/yr 2 for 2003-2012. Bevis et al. (2019) reported that Greenland ice mass had an acceleration of −27.7 ± 4.4 Gt/yr 2 from 2003 to mid-2013 but decelerated after that. The deceleration after mid-2013 was also pointed out by Ruan et al. (2019) and . After 2013, it is unclear whether the strong acceleration in ice loss will persist into the future. Mouginot et al. (2019) reconstructed the mass balance of the GrIS from 1972 to 2018 and found that the acceleration in ice loss switched from positive in 2000-2010 to negative in 2010-2018 due to a series of cold summers. All these results suggest that the acceleration is highly variable and that we should exercise caution when extrapolating long-term trends from short-term records.
We share the opinion of Wouters et al. (2013) that ice sheet mass balance contains both a secular component and a stochastic component. The latter is related to stochastic behaviors of climate variability and can be large in amplitude (several hundred gigatonnes) (Lenaerts et al., 2019). Williams et al. (2014) found that correlated noise is present in GRACE-derived ice mass time series in Antarctica and pointed out this stochastic process is related to ice mass changes. Such strong stochastic variability in mass change can temporarily subdue or amplify any underlying secular signals, especially when the span of the observations is not sufficiently long to eliminate this influence. However, most of the previous studies did not give adequate consideration to the stochastic process in ice mass change, which may result in inaccurate modeling of past mass change and inaccurate extrapolation of future mass change.
To improve our quantification of Greenland ice mass changes and to better understand the mass change process, this study proposes a new model to describe Greenland ice mass changes (including deterministic components and stochastic components). Using this new model, we improve the estimates of the secular mass changes from 2003 to 2017. We show that such improvement is important if one simply extrapolates the observation for future ice change. We also quantitatively assess the impacts of stochastic variations on estimating secular components. We also calculate the Greenland mass changes from Surface Mass Balance (SMB) and ice discharge (ID) data. The SMB model used is the regional climate model RACMO2.3p2 . The RACMO2.3p2 provides 1 × 1 km monthly incremental SMB covering the whole GrIS and peripheral ice caps and 5.5 × 5.5 km monthly cumulative precipitation, evapotranspiration, and runoff covering the surrounding tundra region. We obtain the RACMO2.3p2 SMB for the tundra covered regions by using precipitation minus ablation from evapotranspiration and runoff. We use the monthly ice discharge data (version b) over Greenland from Bamber et al. (2018aBamber et al. ( , 2018b. This new data set includes results from 1958 to 2016 by supplementing and extending Rignot et al. (2008) with a new data set that captures surface velocities at sub-annual temporal intervals for 195 outlet glaciers across Greenland. Ice thickness data are generally obtained from ice penetrating radar data measurements. Ice discharge is the product of surface velocity and ice thickness along an outlet glacier flux gate. The error in ice discharge is conservatively estimated at 6% for years when ice thickness was observed (1992-2016) (Bamber et al., 2012). The sum of the SMB and ID gives the total mass change over Greenland.

A new model for describing temporal changes of Greenland ice mass
The conventional model used to describe Greenland ice mass changes usually consists of a constant, a rate, an acceleration, and four seasonal (annual and semi-annual) terms. However, this conventional model does not take into account the interannual variations in Greenland ice mass, a contribution that has been proven to play an important role in Greenland mass change processes and impacts the estimation of secular rates and accelerations (Wouters et al., 2013;Zhang et al., 2019). Therefore, the conventional model is limited in its ability to describe the Greenland mass changes.
Since the interannual variations are not regular and behave more like a stochastic process, we therefore cannot use deterministic functions to describe them. Here, we improve the conventional model by introducing a stochastic term to describe the interannual variations. Such a model is called augmented stochastic model in Blewitt (1998) and equivalent to the conventional model plus a stochastic term and can be mathematically described as where m is the mass; t is the time in years; a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 are the coefficients for the deterministic components that are invariable with t; ξ s describes the irregular interannual variations and changes with t; and ξ w is the noise that obeys a normal distribution when the mass change is correctly modeled. Using this new model, we can not only estimate the trend, seasonal and interannual variations simultaneously but also assess the impacts of interannual variations on the estimates of other parameters. For simplicity, we refer to this model as the new model hereafter.
We use a First-Order Gauss-Markov (FOGM) process (equivalent to autoregressive AR(1) noise) to describe the interannual variations (ξ s ). The spectrum of an FOGM process is controlled by a variance (σ 2 f ogm ) for the overall power and by a correlation time, τ , for the transition from flat to the decreasing power. As τ approaches zero, an FOGM process degrades into a white noise pro- cess. When τ approaches infinity, an FOGM process becomes a random-walk. Therefore, an FOGM process can describe stochastic processes ranging from white noise to random-walk noise. Practically, the FOGM process can be easily implemented in a Kalman Filter (KF) (Ji and Herring, 2013). In the KF model, a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 are constant, the state transition for ξ s is denoted by ϕ f ogm k and the filter process noise variance is denoted by q f ogm k .
They are related to the FOGM parameters τ and σ 2 f ogm by the Yule-Walker equations: where t k = t k+1 − t k , k orders the observations. We follow Ji and Herring (2013) to use a forward-backward KF to solve Eq. (1). Before running this KF model, we use a 2-D grid search technique to search for the FOGM parameters (Ji and Herring, 2013). In this KF approach, two criteria are used: one is the least squares criterion and the other is the maximum likelihood criterion. Therefore, this approach belongs to the maximum likelihood estimation method. Additional technical details about the determination of the optimized FOGM parameters and the KF model can be found in Ji and Herring (2013).

Determination of the FOGM parameters
When using the 2-D grid search technique (Ji and Herring, 2013) to search for the FOGM parameters, we have difficulty finding the global minimum of the misfit Root Mean Square (RMS). We attribute this to the short time span of the GRACE observations and the aliasing between the quadratic trend and the stochastic process. Fig. 1 shows the misfit RMS when different FOGM variances and correlation time are used. We identify the misfit RMSs that result from normally distributed fitting residuals to indicate which FOGM parameters can be used to separate stochastic process from white noise process. It is observed that when the residuals obey normal distribution, the misfit RMS decreases with larger FOGM variance and increases with longer correlation time. It shows that the fitting residuals start to obey a normal distribution when the correlation time is greater than 180 days. Since the misfit RMS increases with increasing correlation time, we choose 180 days as the correlation time for the GRACE-estimated Greenland mass changes. When the correlation time is fixed at 180 days, the residuals start to obey a normal distribution when the variance is greater than 1600 Gt 2 . Although the misfit RMS decreases with increasing variance, it changes very slowly once the variance is greater than 1600 Gt 2 . We use the new method to estimate the rate and the acceleration using a fixed correlation time (180 days) and different variances. The estimated rates and accelerations are shown in Fig. 2. Fig. 2a shows that the rate increases slowly with increasing variance and the increment change is less than 2 Gt/yr when the variance increases from 100 Gt 2 to 10000 Gt 2 . This finding indicates that the rate estimate is insensitive to the variance parameter. Fig. 2b shows that the estimated acceleration increases with larger variance until the variance is greater than 4000 Gt 2 . After this threshold, the acceleration remains essentially constant. Fig. 3 shows the rate and acceleration estimates derived by using different correlation time and fixing the variance at 4000 Gt 2 . It shows that when the correlation time increases from zero to 1000 days, the rate estimate decreases from −284.5 Gt/yr to −340.6 Gt/yr while the acceleration estimate increases from −7.5 Gt/yr 2 to 11.7 Gt/yr 2 . This indicates that the rate and acceleration estimates are sensitive to the correlation time. Therefore, we should carefully determine the correlation time. Based on the above findings, we choose τ = 180 days and σ 2 f ogm = 4000 Gt 2 as the FOGM parameters.

Simulation tests to validate the new method
We conduct three simulations to evaluate the appropriateness and effectiveness of the new method, as detailed in Sections S1-S3 of the Supplementary Materials (SM). We create synthetic observations based on Eq. (1). The reference deterministic components are obtained by fitting Eq. (1) without stochastic components to the GRACE-estimated GrIS mass changes from 2003 to 2017, and are shown in Table S1 of the SM. The stochastic component includes an FOGM process and white noise. The synthetic observations are calculated by adding different FOGM processes and white noise to the reference deterministic component. We test the performance of the new model in comparison with the conventional model in three cases by each time only varying: (1) the FOGM correlation time, (2) the FOGM variances, or (3) the length of the observations.
In Test 1 of the SM, we evaluate the Mean Absolute Errors (MAEs) of the estimated rate, acceleration, and seasonal amplitudes when the FOGM correlation time varies but the FOGM variance and the white noise variances are fixed. Results show that all the parameters estimated by the new method have smaller MAEs than those estimated by the conventional method and the improvement signifies with increasing correlation time (up to 900 days). These results suggest that the new method provides better estimates of rate, acceleration, and seasonal amplitudes than the conventional method. The most significant MAE improvement appears in the acceleration determination, which, on average, can be as high as 14.1% of the acceleration itself. This result suggests that the stochastic variations are more likely to impact the estimation of acceleration than the other parameters. This is explained by the fact that the stochastic variations in the Greenland ice mass change may have similar behaviors with the quadratic term (see the example presented in Fig. S2 of SM). Fig. S1 of the SM also shows that the MAEs of the rate and acceleration estimates from both methods increase with increasing correlation time up until a correlation time of 900 days. This result suggests that larger correlation times have a greater impact on the estimates of rate and acceleration than smaller correlation times. Longer than 900 days, the MAEs of estimated rate and acceleration stay relatively stable. However, the MAEs of annual and semiannual amplitudes decrease with increasing correlation time, which suggests that larger correlation times have less impact on the estimates of both amplitudes. It should be noted that the impacts associated with the correlation time also depend on the length of the observations.
In Test 2 of the SM, we evaluate the MAEs of the rate, acceleration, and seasonal amplitudes when the FOGM variance varies but the FOGM correlation time and the white noise variance are fixed. Fig. S3 of the SM shows that the MAEs derived from the new method are smaller than those derived from the conventional model and the improvement becomes greater with increasing FOGM variance. This result indicates again that the new method outperforms the conventional model when FOGM process noise is present in the time series. We also confirmed the expectation that MAEs derived from both methods increase with increasing FOGM variance.
In Test 3 of the SM, we investigate the impact of data length on the parameter estimates. We use time spans ranging from 3 to 50 yr to estimate the rate, acceleration, and seasonal amplitudes. Fig. S4 of the SM shows again that the rates, accelerations, and seasonal amplitudes estimated by the new method have smaller MAEs than those estimated by the conventional method. Furthermore, the MAE differences become smaller as the data time span increases. This result indicates that the new method yields better estimates of rate, acceleration, and seasonal amplitudes than the conventional method, but this improvement decreases with increasing data time span. Fig. S4 of the SM also shows that the MAEs of all the parameters estimated by both methods decrease with increasing data length. It needs ∼20 yr for the MAE of rate to drop below 5 Gt/yr and ∼15 yr for the MAE of acceleration to drop below 5 Gr/yr 2 .
However, the MAEs of the annual and semiannual amplitudes decrease quickly in the first 20 yr but decreases very slowly or even keep constant after that.
In summary, the new method provides improved estimates of the rate, acceleration, annual and semiannual amplitudes than the conventional method, especially when the data length is short (e.g., <20 yr). The errors in the above estimates become smaller and even negligible when the data length gets longer, which is especially true for the rate and acceleration estimates.

Results from GRACE Mascon Data and SMB plus ID data
When using the FOGM parameters determined in Section 2.3 to constrain the stochastic process and use the GRACE data to fit and solve Eq. (1), we successfully separate the trend (linear and quadratic terms), the seasonal variations, and the interannual variations. The Greenland mass change time series and the best-fit mass trajectory model are shown in Fig. 4. The new model fits well with the GRACE measurements having a misfit RMS of only 14.5 Gt. The rate and acceleration are estimated to be −288.2 ± 12.7 Gt/yr and −1.6 ± 1.3 Gt/yr 2 . When we use the conventional model to run the KF (in this case the KF is equivalent to the least square method), the rate and acceleration estimates become −284.5 ± 1.4 Gt/yr and −7.5 ± 0.2 Gt/yr 2 . The difference between estimates from the two models is small for the rate but is significant for the acceleration. This suggests that the stochastic component has a significant impact on the estimate of acceleration, which is consistent with our simulations. Fig. 5 shows the separated seasonal component, interannual component, and the fitting residuals. The peak-to-peak amplitude of the seasonal variations is 335.7 ± 11.4 Gt, which is greater  than the long-term mass loss rate. The interannual variations are characterized by a gradual increase from 2003 to 2009, a sharp decrease from 2009 to 2013, then followed by an increase from 2013 to 2017. The peak-to-peak amplitude of the interannual variations is as high as 568.8 ± 40.1 Gt, which is much greater than the long-term rate and the seasonal amplitude. Fig. 6 shows the Greenland mass change time series estimated from the SMB plus ID data and its best-fit mass trajectory model. It shows that the new model fits well with the modeled mass change data. The rate of the mass change is estimated to be −274.9 ± 13.0 Gt/yr and the acceleration is estimated to be −1.1 ± 1.3 Gt/yr 2 . The mass loss rate from SMB plus ID data is 13.3 Gt/yr lower than the GRACE-estimated rate and the mass loss acceleration is 0.5 Gt/yr 2 smaller than the GRACE-estimated acceleration. These differences are within the uncertainty level. Fig. S5 of the SM shows the mass change time series derived from SMB, ID, and their sum, and Fig. S6 of the SM shows the detrended ones. They together show that the interannual and seasonal variations of mass changes were mainly contributed from SMB with little contribution from ID. This may indicate that the stochastic climate variability has more impacts on the SMB than ID over the GrIS. Fig. S7 of the SM shows the separated seasonal component, interannual component, and the fitting residuals derived from the SMB plus ID. The FOGM component shows the same pattern as that from the GRACE result as shown in Fig. 5.

Impacts of stochastic variations on estimates of deterministic parameters
As highlighted in Section 4 that the interannual variations have a much greater amplitude than the secular rate in 2003-2017, it will certainly bring impacts on the estimates of secular terms, especially when the time span of observations is not long enough to eliminate this influence. When investigating the long-term trend using data only ranging from 2003 to 2013, we may easily estimate a large mass loss acceleration (−32.8 ± 1.9 Gt/yr 2 using the conventional model) because the interannual variations from 2003 to 2013 resemble quadratic variations (Fig. 5b). But if we extend the time span beyond 2013, the estimated acceleration would decrease since the interannual variation pattern changed and started to decelerate after 2013. Table 1 compares the rate and acceleration estimates from the new and the conventional models using the GRACE-estimated and the SMB plus ID mass changes. Although the estimates based on the GRACE and SMB plus ID data are different, these differences are relevant to the input data and irrelevant to methods. In this section, we focus on the differences caused by the estimating method itself. The differences between estimates from different methods in rate or acceleration parameter are no more than 4 Gt/yr or 6 Gt/yr 2 ( Table 1). The differences are small for the rate given their magnitudes but are significant for the acceleration since the magnitude of the acceleration itself is only a few Gt/yr 2 . Taking the results based on the GRACE data for example: the conventional method estimates the acceleration to be −7.5 ± 0.2 Gt/yr 2 , which is actually in the lower range of acceleration estimates; whereas the new method estimates the acceleration to be −1.6 ± 1.3 Gt/yr 2 . The difference in the acceleration estimate can be as large as 5 Gt/yr 2 . If simply extrapolating the same quadratic trend into the future, such difference of 5 Gt/yr 2 in acceleration would introduce a bias of ∼16,000 Gt in the GrIS mass balance and a resultant bias of ∼44 mm in global sea level rise between 2020 and 2100.
When referring the uncertainties in Table 1, we note that the uncertainties of rate and acceleration of the new method are much larger than that of the conventional method. This is reasonable since the conventional method assumes white noise in the observations and thus estimates the least uncertainties. Williams et al. (2014) pointed out that the uncertainties of rate and acceleration estimates of Antarctica ice mass change should be scaled up when we use correlated noise model instead of white noise. This is also applicable to Greenland ice mass changes. According to the theory of Williams et al. (2014), the significant difference in rate and acceleration uncertainties between the new and the conventional method implies strong stochastic mass variations, which is consistent with our previous findings in Fig. 5. The studies mentioned in the introduction did not separately consider the acceleration and the stochastic process and take the stochastic process as a quadratic trend, therefore they actually estimate the lower bounds of the acceleration, i.e., they overestimated the mass loss acceleration. Using this new method, we simultaneously estimate the acceleration and the stochastic variations that result in much smaller mass loss acceleration. As the length of GRACE time series is not sufficiently long to eliminate the impacts from stochastic process, using the new method can improve the estimates of the deterministic parameters, especially the acceleration estimate, which is important for extrapolating future mass changes.

Discussion
To fully reveal the impacts of stochastic processes on the estimates of rate and acceleration of the Greenland ice mass change, here we give the lower and upper bounds of these estimates. Since there is a strong aliasing between the acceleration and the stochastic component, we assumed two extreme situations: (1) the acceleration does not exist, i.e., the quadratic trend is only part of a stochastic process; (2) the acceleration does exist and the stochastic processes are neglected.
Under case (1), we remove the quadratic term from Eq. (1) and use the FOGM process to represent the possible stochastic process and quadratic trend. We first reestimate the FOGM parameters and then use these FOGM parameters to estimate the rate and acceleration using the same method as in Section 2, whereby we obtain the upper boundary of the acceleration estimate (5.3 ± 4.0 Gt/yr 2 from the GRACE data and 0.7 ± 2.2 Gt/yr 2 from the SMB plus ID data). The reestimated FOGM parameters and the corresponding rate and acceleration estimates are listed in Table 2. In case (2), we essentially use the conventional mass change model. We have estimated the rate and the acceleration, which represents the lower bound of the acceleration estimate (−7.5 ± 0.2 Gt/yr 2 from the GRACE data and −3.5 ± 0.2 Gt/yr 2 from the SMB plus ID data).
Based on the results from the two extreme cases, we find that the possible range of acceleration of Greenland mass change from 2003 to 2017 is estimated to be −7.5 to 5.3 Gt/yr 2 from the GRACE data and −3.5 to 0.7 Gt/yr 2 from the SMB plus ID data. The corresponding rate is estimated to be −284.5 to −301 Gt/yr from the GRACE data and −273.0 to −277.6 Gt/yr from the SMB plus ID data. As expected, estimates of the rate and the acceleration by the new method in Table 1 fall in the middle of the extreme values.
As revealed by this exercise, the conventional model, which does not consider the stochastic processes in Greenland mass loss, provides an estimate of the lower boundary of acceleration and the upper boundary of rate. In this study, we assume the interannual variation to be a stochastic component of mass change. This stochastic mass change is due to climate variability  and if not well accounted for, will add uncertainty to extrapolation of present-day changes to estimates of future mass loss (Wouters et al., 2013).

Conclusions
We propose a new mass trajectory model to describe Greenland mass change by introducing a stochastic process to describe the interannual variations. Simulations verify that this new model better estimates the rate, the acceleration, and the seasonal components than the conventional model, and can also reliably extract the interannual variations.
Our results show that the interannual variability is an important component in Greenland mass change and can greatly impact estimates of the long-term rate and acceleration. The conventional model neglects this component and thus represents the lower bounds of acceleration and the upper bounds of rate. Using the new model, we estimate the acceleration to be −1.6 ± 1.3 Gt/yr 2 from the GRACE data and −1.1 ± 1.3 Gt/yr 2 from the SMB plus ID data for 2003-2017. The corresponding rate is estimated to be −288.2 ± 12.7 Gt/yr and −274.9 ± 13.0 Gt/yr, respectively. Under extreme assumptions, the acceleration is estimated to be −7.5 to 5.3 Gt/yr 2 from the GRACE data and −3.5 to 0.7 Gt/yr 2 from the model data. The corresponding rate range is −284.5 to −301 Gt/yr and −273.0 to −277.6 Gt/yr, respectively.
This new model improves the estimates of Greenland mass change and changes our understanding of Greenland mass change patterns. The presence of the stochastic mass change is due to climate variability and has nonnegligible impacts on long-term rate and acceleration estimations. We should pay more attention to the stochastic mass change and its linkage to climate variability and exercise more caution when extrapolating future ice mass change and sea level rise from short-term observations.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.