Constraining the mass balance of East Antarctica

We investigate the mass balance of East Antarctica for the period 2003–2013 using a Bayesian statistical framework. We combine satellite altimetry, gravimetry, and GPS with prior assumptions characterizing the underlying geophysical processes. We run three experiments based on two different assumptions to study possible solutions to the mass balance. We solve for trends in surface mass balance, ice dynamics, and glacial isostatic adjustment. The first assumption assigns low probability to ice dynamic mass loss in regions of slow flow, giving a mean dynamic trend of 17 ± 10 Gt yr−1 and a total mass imbalance of 57 ± 20 Gt yr−1. The second assumption considers a long‐term dynamic thickening hypothesis and an a priori solution for surface mass balance from a regional climate model. The latter results in estimates 3 to 5 times larger for the ice dynamic trends but similar total mass imbalance. In both cases, gains in East Antarctica are smaller than losses in West Antarctica.


Introduction
Obtaining reliable estimates of the current mass balance of the ice sheets is needed for predicting their future evolution and contribution to sea level rise. While the Greenland Ice Sheet is known to be experiencing an increasing mass loss since the 1990s [Sasgen et al., 2012], uncertainties in the mass balance of Antarctica remain large [Hanna et al., 2013]. The East Antarctic Ice Sheet (EAIS) has the greatest potential to contribute to sea level rise in the future. It is, however, the most challenging ice mass to study because of factors such as the relative paucity of in situ observations and the poor signal-to-noise ratio of satellite data related to mass balance.
Several studies have investigated the mass balance of the EAIS using different techniques. Using the Center for Space Research Release 5 monthly Gravity Recovery and Climate Experiment (GRACE) solution, Williams et al. [2014] found a trend of 97 ± 13 Gt yr −1 between March 2003 and July 2012 after correcting for glacio-isostatic adjustment (GIA) using the W12a model [Whitehouse et al., 2012]. For the period August 2002 to July 2016, the gravimetric mass balance (GMB) product from TU Dresden [Groh and Horwath, 2016] over the EAIS gave a trend of 63.9 ± 30.8 Gt yr −1 .
Employing radar altimetry data from the European Remote Sensing satellites ERS-1 and ERS-2, Davis et al. [2005] found a mass gain of 45 ± 19 Gt yr −1 from 1992 to 2003 using a near-surface snow density of 350 kg m −3 . Zwally et al. [2005] found a slightly lower mass gain of 16 ± 11 Gt yr −1 between 1992 and 2002. The ice mass balance intercomparison exercise (IMBIE) [Shepherd et al., 2012] reported an overall mass balance of 58 ± 31 Gt yr −1 for 2003-2008, by taking the arithmetic average of estimates from different techniques.
A more recent study using a Bayesian hierarchical framework combining altimetry, gravimetry, and GPS measurements estimated a mass balance of 56 ± 18 Gt yr −1 for the period 2003-2013[Martín-Español et al., 2016a. For the most recent years 2010-2013using Cryosat-2 data, McMillan et al. [2014 estimated a mean elevation change of 0.1 ± 0.2 cm yr −1 , which they inferred as a mass change of −3 ± 36 Gt yr −1 . V. Helm (personal communication, 2016 after a review of Helm et al. [2014]) derived a change in volume of −15 ± 60 km 3 yr −1 for the same period. and −25 ± 15 Gt yr −1 for the two periods, respectively. The large mass gain over the EAIS was primarily a consequence of the assumption that the elevation increase was due to a long-term dynamic thickening in response to an increase in snowfall since the start of the Holocene. Thus, most of these studies agree that there has been a positive elevation trend over the EAIS but disagree on (i) its magnitude and (ii) its origin, and, as a consequence, (iii) the density that needs to be assigned to convert the observed height change to mass change. The challenge over the EAIS is that the elevation rate is small and close to (or possibly below) the combined errors of the observations. A 1 mm elevation change, averaged over the EAIS, equates to a volume change of about 10 km 3 . If this is due to ice dynamics, then the inferred mass change is 9.2 Gt. If it is due to snowfall variations, it is closer to 3.5 Gt. If it is due to changes in firn compaction rate, the mass change is zero. In reality, the total mass change is some combination of all three of these factors plus a poorly constrained component arising from GIA.
The purpose of this study is to investigate the origin of elevation changes and mass trends over the EAIS by considering three experiments where we apply an inversion from satellite data using different assumptions. In the first experiment, we assume that dynamically driven elevation changes are correlated with high surface velocities [Hurkmans et al., 2012] and are, therefore, not expected to occur in regions of slow ice flow. In the second and third experiments, we use a hypothesis based on Zwally et al. [2015], which allows for dynamic thickening anywhere across the EAIS. For the second and third experiments, as in Zwally et al. [2015], we need to prescribe the surface mass balance (SMB) anomalies. We do this using two versions of a regional climate model (RCM), the Regional Atmospheric Climate Model (RACMO) version 2.3 [van Wessem et al., 2014] and RACMO2.4 (unpublished). We find differences in SMB anomalies implied by the two RCM versions with the more recent version being more consistent with the combined satellite data.

Data and Methods
We use a Bayesian hierarchical framework to resolve spatiotemporal mass trends for the EAIS, described and tested in previous studies [e.g., Zammit-Mangion et al., 2015a, 2015bSchoen et al., 2015;Martín-Español et al., 2016a]. We simultaneously solve for surface processes (SMB and firn densification), ice dynamics, and GIA in a probabilistic inversion of remote sensing observations from different techniques: altimetry, gravimetry, and GPS. Since we have more unknowns than observations, this is an underdetermined problem. Therefore, we need to apply smoothness constraints and prior beliefs (in space and in time) on the unknown processes by appropriately configuring the covariance matrices. The spatiotemporal wavelength parameters are obtained from physical numerical models and other prior information to help apportion the mass change to the different processes in what is termed "source separation." By configuring the covariance matrices, we can separate the signals if they have different spectral characteristics. Alternatively, we can apply restrictive estimates (e.g., a regional climate model (RCM) output for SMB anomalies) to predetermine the mean of a process. In this case, a biased prior estimate in one process will translate into biased posterior estimates in all the other processes.
Three altimetry data sets are used to obtain time series of elevation change for the period 2003-2013. Along-track elevation measurements from ICESat (data release 33) provide high-resolution altimetry observations for the period February 2003 to October 2009. These data were corrected for range determination from Transmit-Pulse Reference-Point Selection (Centroid versus Gaussian) [Borsa et al., 2014] and for the intercampaign bias (ICB) [Hofton et al., 2013]. Envisat radar altimetry data were available for the period January 2003 to November 2010, and elevation rates were obtained following Flament and Rèmy [2012]; Cryosat-2 provided elevation rates for the period July 2010 to December 2013. The rates were derived from level 1B data processed as described in Helm et al. [2014]. We produced spatially resolved annual elevation rates using a 3 years moving window from ICESat and Envisat and a single trend with acceleration coefficients from Cryosat-2.
Data from the GRACE mission were used to estimate annual mass trends for the period January 2003 to December 2013. We used the latest release of the NASA Goddard Space Flight Center global mass concentration (mascon) solution, RL2v15, following Luthcke et al. [2013]. Subannual variability was removed from these data using the ensemble empirical mode decomposition adaptive filtering approach described in Loomis and Luthcke [2014] and Luthcke et al. [2013] prior to trend estimation. Annual mass trends were obtained for each mascon by fitting linear trends to the filtered data in each year.
A network of 80 elastic-corrected GPS stations was used to constrain elevation changes due to GIA. Processing was carried out as described in Martín-Español et al. [2016a]. The elastic effects from surface loading changes were estimated using the Regional Elastic Rebound Calculator [Melini et al., 2015] and removed from the GPS uplift rates using the loading inputs derived in Martín-Español et al. [2016a].
The RACMO versions 2.3 [van Wessem et al., 2014] and 2.4 (unpublished) are used to constrain SMB in the second and third experiments. SMB trends modeled from RACMO are also used to extract the spatiotemporal length scales describing the spatial smoothness of SMB anomalies for 2003-2013 (with respect to the 1979-2008 mean) as described in Zammit-Mangion et al. [2015b].
Prior information is used to help separate the different processes. For example, we assume that changes in ice dynamics are temporally smooth when compared to SMB anomalies (i.e., the "weather" signal oscillates from 1 year to another and also at a subannual time scale while changes in ice dynamics have not been seen to have similar variability anywhere in Antarctica and glaciological theory would preclude such behavior) in order to obtain a decomposition of the observed signal. This prior information appears in the spatiotemporal covariance matrices of the different processes. We adopt this approach since it is likely that the geophysical models exhibit first-order systematic biases that are hard to quantify (as illustrated, for example, by the differences between RACMO 2.3 and 2.4, see supporting information). There is, however, broader confidence in their second-order properties such as spatiotemporal length scales.

Prior Assumptions
We carry out three experiments where two different sets of prior assumptions are considered:

Experiment I (E I). In this experiment we follow the assumptions made in Zammit-Mangion et al. [2015a] and
Martín -Español et al. [2016a]: (1) changes in height due to surface processes (primarily snowfall) have high temporal variability at large (ice sheet) scales [Rignot et al., 2011] and exhibit spatial coherence. Both the spatial and temporal length scales were extracted from RACMO. Firn densification processes are assumed to be highly correlated to surface processes locally. We capture this interaction by finding the empirical covariance matrix between the surface and firn processes from RACMO and a firn densification model [Ligtenberg et al., 2011].
(2) Changes in height due to ice dynamics are assumed to be linear over the relatively short time scales considered and to have short spatial length scales (i.e., individual ice streams behave independently). In addition, a soft constraint is used, so that dynamically driven height changes are less likely (but not impossible) in areas of low velocities and more probable for velocities greater than 10 m yr −1 .
(3) GIA processes are assumed to be temporally invariant and have long spatial wavelengths (from 500 km in West Antarctic Ice Sheet to 1500 km in EAIS). 2. Experiments IIa and IIb. These experiments attempt to replicate the hypothesis that elevation changes in the interior of the EAIS are due to an increase in accumulation rate at the beginning of the Holocene, which has not yet fully reached equilibrium with ice motion, following Zwally et al. [2015]. To implement this hypothesis, we apply the following modifications to the assumptions used in E I: (1) given that the interior of EAIS is characterized by low horizontal velocities at the surface, we remove the soft (probabilistic) constraint imposed in E I for the dynamically driven elevation changes, (2) elevation changes due to SMB processes are not coestimated through the statistical inversion, but predefined using a RCM output. We compare two different versions of RACMO: RACMO2.3 in E IIa and RACMO2.4 in E IIb. A conservative uncertainty of 20% is added to the SMB anomalies [Wouters et al., 2015].
In all experiments, height changes are converted into mass changes via the following density assumptions: upper mantle density at 3800 kg m −3 over West Antartica and 4500 km m −3 over EAIS, ice density at 917 kg m −3 , and SMB values ranging from 350 to 600 kg m −3 , using the density map of Ligtenberg et al. [2011].

Estimates of Mass Balance of East Antarctica: Results From Experiments
The SMB and ice dynamics trends over the EAIS for each experiment are shown in Figures 1a and 1b. We also estimate a time-invariant solution for GIA for each case (see section 3.3). Note that as the net mass change is largely constrained by GRACE data, the actual impact of the assumptions considered on each experiment reflects how the total mass signal is apportioned between the individual processes. Mean mass balance values for each component are given in Table 1.
The assumptions used in E I lead to small positive dynamic trends with a mean rate of 17 ± 10 Gt yr −1 , despite the constraint imposed for the ice dynamics. By removing the constraint over low-velocity regions in E IIa and E IIb, the dynamic signal over the EAIS increases substantially. However, the version of RACMO used to  Figure 1c shows the gravimetric mass balance product in red for comparison. The 1 confidence interval is given by the shading.
constrain SMB considerably affects the results: RACMO2.3 simulates more negative SMB anomalies for the period 2003-2013 (with a mean rate of −31 Gt yr −1 over the EAIS) than the newer version RACMO2.4 (mean rate of 3 Gt yr −1 ). As a consequence, the ice dynamics trend in E IIa is considerably stronger, with an average rate of 80 ± 6 Gt yr −1 compared with E IIb (54 ± 6 Gt yr −1 ). However, irrespective of which version of the RCM is used, we obtain smaller ice dynamic trends than those derived by Zwally et al. [2015] (147 Gt yr −1 over the periods 1992-2001 and 2003-2008 following similar assumptions). Differences in net mass balance between Experiments E IIa and E IIb (Figure 1c) can be explained by the large differences in SMB estimates. It demonstrates the importance of the RCM SMB reconstructions when they are used to constrain the SMB processes over vast regions such as the EAIS.
We compare the net mass balance derived from the different experiments with independent time series from gravimetric mass balance (GMB) over the EAIS (Figure 1). The GMB product over the EAIS agrees well with the results from E I. E IIa shows the largest discrepancies with the GMB, particularly in the years 2006, 2007, and 2010. Version 2.4 of RACMO adjusts the previous model's SMB anomalies between 2005 and 2007, bringing them into closer agreement with results from E I. However, in 2010 it is not clear whether the discrepancy is due to a positive bias in the RCM or due to a sudden ice shelf breakup, triggering negative ice dynamics trends, which, due to their short temporal wavelength, have been incorrectly attributed to SMB in E I. Figure 2 shows the spatial distribution of mean elevation changes between 2003 and 2013 due to SMB, ice dynamics, and GIA, over Antarctica, as obtained from each experiment. Spatial patterns of elevation changes due to ice dynamics (dh d ∕dt) are similar in the two experiments, but the regions showing dynamic thickening are slightly larger, both in magnitude and in spatial extent, in E II than in E I. On the other hand, we find large differences in the pattern of SMB-driven elevation changes (dh s ∕dt) between E I and E II. Recall that in E I, SMB anomalies are coestimated in the inversion process while in E II, these are predefined from RCM outputs. The largest differences are found over Victoria-Wilkes Land region (basins 11-13, using definitions from Sasgen et al. [2013] shown in Figure S1) ( Figure S2). Annual SMB estimates from E I and from model outputs RACMO2.3 and RACMO2.4 for the Victoria-Wilkes Land region are shown in Figure S3. In that region RACMO2.3 estimates a strong negative trend of −49 Gt yr −1 , which does not match in situ measurements (as noted in Wang et al. [2016]) nor satellite observations ( Figure S4). However, in RACMO2.4 the trend over the  same region is reduced to −27 Gt yr −1 . From the Bayesian inversion, with no predefined SMB, we obtain slightly positive trends (15 ± 13 Gt yr −1 ) for the entire region. This is consistent with the increase in accumulation rates reported by Frezzotti et al. [2013] and the values modeled by the Modern-Era Retrospective Analysis for Research and Applications reanalysis product [Bosilovich et al., 2008;Rienecker et al., 2011].

Experiment II Using Only Altimetry
To allow a more direct comparison with Zwally et al. [2015], we repeat E IIb excluding GRACE data. In this case, the resulting mass trends are only dependent on the input altimetry data (ICESat, Envisat, and Cryosat-2) and the length scales and assumptions which are used as prior information (GIA is fixed to that obtained in E IIb, see section 3.3). A comparison between these mass trends and those estimated including GRACE is shown in Figure 3. We find that by using only altimetry data, the mean mass trend (over 2003-2013) becomes smaller (39 ± 10 Gt yr −1 ) than when GRACE data are included in the statistical inversion (57 ± 7 Gt yr −1 ), therefore increasing discrepancies with the results found in Zwally et al. [2015] (note that their study period was 2003-2008). However, without GRACE, apportioning between SMB and ice dynamics processes is challenging unless a RCM output is used. Furthermore, altimetry data might be subject to systematic errors such as snowpack penetration and the ICB for ICESat.
The IMBIE study, covering the period October 2003 to December 2008 [Shepherd et al., 2012], includes four different ICESat trends for the EAIS, averaging 78 ± 19 km 3 yr −1 (or 109 ± 57 Gt yr −1 based on SMB and firn compaction models), while estimates derived from Envisat averaged 22 ± 39 Gt yr −1 . The difference may be due to penetration effects of the radar signal and how it is corrected for in the preprocessing stage.
Our altimetry-only trend is less positive than ICESat trends estimated by Zwally et al. [2015] and IMBIE. One reason could be the use of radar altimetry (Envisat and Cryosat-2) data in the inversion, which observe lower trends. Another reason could be the different ICB corrections applied to ICESat data in the different studies. We use the Hofton et al. [2013] ICB correction which corrects for a trend of 1.04 ± 0.48 cm yr −1 (over the entire ICESat period). After applying this correction, we obtain a net ICESat-volume change of 40 km 3 yr −1 over the EAIS. This correction differs from that used in the IMBIE study, 0.65 cm yr −1 , and even more from that used by Zwally et al. [2015], −1.43 cm yr −1 (whose large discrepancy largely explains the positive bias of their trends when compared to the others). Note that these two studies consider the period 2003-2008. The equivalent of our ICB correction for that period would be 0.87 cm yr −1 [Hofton et al., 2013]. If compared to that used in the IMBIE study, Hofton's ICB correction would explain a difference of 22 km 3 yr −1 , enough to reconcile our ICESat-derived volume change with that from IMBIE. Evidently, the ICB correction has a strong impact in the EAIS, and any error in this correction results in a biased estimate of the EAIS mass balance using ICESat alone [Scambos and Shumman, 2016].

Estimated Glacial Isostatic Adjustment
For each experiment we coestimate a GIA solution along with the ice sheet mass balance ( Figure 2). GIA estimates from E II show a strong subsidence in the interior of the EAIS. This is a consequence of the different mass spatial distribution at the surface due to the strong ice dynamic trends. Interestingly, it compares well with the forward GIA solution W12 [Whitehouse et al., 2012], which shows a subsidence pattern over the same region (although with a slightly larger magnitude). The total GIA-induced mass change over the EAIS is 23 ± 6 Gt yr −1 for E I, 12 ± 4 Gt yr −1 in E IIa, and 5 ± 4 Gt yr −1 in E IIb. To support the assumption of the large dynamic thickening trend while reconciling GRACE data, Zwally et al. [2015] proposed an average GIA for the EAIS of −1.2 mm yr −1 . This was obtained by removing 1.6 mm yr −1 from the forward GIA model IJ05R2  as a result of considering an additional ice loading during the Holocene. This would be equivalent to a mass change of −54 Gt yr −1 , which is much more negative than any result derived from available GIA models. The average GIA over the EAIS, as obtained from several forward and inverse models, is equivalent to 25 ± 18 Gt yr −1 when using a rock density of 4500 kg m −3 [Martín-Español et al., 2016b] ( Figure S5). Furthermore, it should be noted that models such as W12 and IJ05R2 already account for Holocene accumulation increases.

Conclusions
We apportion mass trends to SMB and ice dynamics for the EAIS, based on two different assumptions, different remote sensing data, and two RCMs. In the first experiment, the model apportions about a third of the mass trend to ice dynamics, +17 Gt/yr, and two thirds, +40 Gt yr −1 to SMB, resulting in a total mass trend for the 10.1002/2017GL072937 EAIS of 57 ± 20 Gt yr −1 . In the second and third experiments, we allow dynamic-thickening trends in regions of low velocities, replicating dynamic thickening hypothesis proposed in a recent study [Zwally et al., 2015]. We obtain a ∼5 times larger dynamic trend (80 ± 6 Gt yr −1 ) when SMB is constrained using RACMO2.3 than EI, but this model is now known to produce negatively biased estimates over large parts of the EAIS. When forcing SMB with RACMO2.4, we obtain a ∼3 times larger dynamic trend (54 ± 7 Gt yr −1 ) than E I.
Our experiment using only altimetry data (ICESat, Envisat, and Cryosat-2) agrees well with the estimates produced by combining different techniques; therefore, the discrepancy with Zwally et al. [2015] estimates cannot be attributed to the use of GRACE. Neither can it be reconciled when applying the ICB correction from Zwally et al. [2015] over ICESat data, if we combine those with Envisat and force SMB with RACMO2.4.
The coestimated GIA-induced mass changes range from 5 to 23 Gt yr −1 for the different experiments. These are plausible values within the range of other published GIA solutions. The GIA-induced mass rate needed to reconcile the mass change estimates from Zwally et al. [2015] with those from gravimetry data is much smaller than our smallest estimate.
We are unable to reproduce the large magnitude trend obtained in Zwally et al. [2015] using the experiments presented in this paper. We conclude that irrespective of any assumption made about the density of surface elevation changes, mass gains in the EAIS do not exceed mass losses for the West Antarctic Ice Sheet over the studied period. Furthermore, the differences in SMB estimates found between the different versions of RACMO demonstrate the importance of improving the accuracy of the SMB models if they are to be used to constrain SMB processes over vast regions such as the EAIS.