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.

1 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–2013 using 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 km3 yr−1 for the same period. A recent study, utilizing a combination of satellite radar and laser altimetry, assessed the mass balance of Antarctica for two epochs: 1992–2001 and 2003–2008 [Zwally et al., 2015]. The authors derived a mean elevation change of 1.1 and 1.3 cm yr−1, respectively, for the EAIS. They concluded that the Antarctic ice sheet had experienced a net mass gain throughout both epochs due to mass gains over the EAIS of 136 ± 28 Gt yr−1, which outweighed the losses over West Antarctica that were estimated to be −15 ± 20 Gt yr−1 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 km3. 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.

2 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, 2015b; Schoen 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. 2.1 Prior Assumptions Experiment I (E I). In this experiment we follow the assumptions made in Zammit‐Mangion et al. [ 2015a Martín‐Español et al. [ 2016a Rignot et al., 2011 Ligtenberg et al., 2011 −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). 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 Wouters et al., 2015 We carry out three experiments where two different sets of prior assumptions are considered: 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].

3 Results and Discussion 3.1 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. Figure 1 Open in figure viewer PowerPoint σ confidence interval is given by the shading. Time series of (a) SMB, (b) ice dynamics, and (c) net mass trends for EAIS over the period 2003–2013 obtained from the different experiments. Mass trends derived from E I are shown in black, those from E IIa (where SMB has been forced with RACMO2.3) in purple, and those from E IIb (forced with RACMO2.4) in light blue. Figure 1 c shows the gravimetric mass balance product in red for comparison. The 1confidence interval is given by the shading. Table 1. Average Mass Trends of the EAIS (in Gt yr−1) Over the Period 2003–2013 Assumption 1 Assumption 2 Process E I E IIa E IIb Only Altimetry Ice dynamics 17 ± 10 80 ± 6 54 ± 7 25 ± 9 SMB 40 ± 18 −31 ± 6a 3 ± 1a 14 ± 3a Mass imbalance 57 ± 20 49 ± 8 57 ± 7 39 ± 10 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 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]. Figure 2 Open in figure viewer PowerPoint Average elevation changes between 2003 and 2013 due to SMB (dh s /dt) ice dynamic (dh d /dt) and GIA (dh b /dt) processes for experiments E I, E IIa, and E IIb, where E IIa is forced by RACMO2.3 and E IIb by RACMO2.4. 3.2 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. Figure 3 Open in figure viewer PowerPoint Time series of the mass balance of East Antarctica using the assumptions taken in Experiment IIb. Red bars consider inversions with only altimetry data, while the grey bars also include GRACE data. Error bars indicate the 1σ and 2σ intervals in dark and light shadings, respectively. 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 km3 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 km3 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 km3 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]. To further explore the differences with Zwally et al. [2015], we make an additional experiment following E IIb using only altimetry data, where Zwally et al.'s [2015] ICB correction is applied to ICEsat data instead of Hofton et al. [2013]. In this case, the total mass change averages 62 ± 7 Gt yr−1. As expected, this value is larger than the one obtained using Hofton's correction, but still half the estimate given in Zwally et al. [2015] (136 ± 28 Gt yr−1 over 2003–2008). Differences could be a combination of (1) the use of Envisat as well as ICEsat over the period (2003–2009) and (2) the use of the newest RCM for Antarctica (RACMO2.4) which estimates a mean SMB anomaly of 3 ± 1 Gt yr−1 (as opposed to −11 ± 6 Gt yr−1). 3.3 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 [Ivins et al., 2013] 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.

4 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 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.

Acknowledgments This work was made possible by the NERC grant NE/I027401/1 and the Advanced ERC Grant GlobalMass, 694188. The authors would like to acknowledge the following colleagues who provided the data; without them the project would not have been possible: Melchior van Wessem and Michiel van der Brooke for providing the RCM outputs and Bert Wouters, Veit Helm, Scott Luthcke, Frédérique Rémy, Thomas Flament, Liz Petrie, and S. Ligtenberg for providing the different data sets used in this study. Remote sensing data sets used and results derived this study are available by contacting the corresponding author. The regional climate model outputs are available upon request by contacting J.M. van Wessem (Institute for Marine and Atmospheric Research).

Supporting Information Filename Description grl55780-sup-0001-supplementary.pdfPDF document, 704.3 KB Supporting Information S1 Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.