Simultaneous solution for mass trends on the West Antarctic Ice Sheet

The Antarctic Ice Sheet is the largest potential source of future sea-level rise. Mass loss has been increasing over the last 2 decades for the West Antarctic Ice Sheet (WAIS) but with significant discrepancies between estimates, especially for the Antarctic Peninsula. Most of these estimates utilise geophysical models to explicitly correct the observations for (unobserved) processes. Systematic errors in these models introduce biases in the results which are difficult to quantify. In this study, we provide a statistically rigorous error-bounded trend estimate of ice mass loss over the WAIS from 2003 to 2009 which is almost entirely data driven. Using altimetry, gravimetry, and GPS data in a hierarchical Bayesian framework, we derive spatial fields for ice mass change, surface mass balance, and glacial isostatic adjustment (GIA) without relying explicitly on forward models. The approach we use separates mass and height change contributions from different processes, reproducing spatial features found in, for example, regional climate and GIA forward models, and provides an independent estimate which can be used to validate and test the models. In addition, spatial error estimates are derived for each field. The mass loss estimates we obtain are smaller than some recent results, with a time-averaged mean rate of−76± 15 Gt yr for the WAIS and Antarctic Peninsula, including the major Antarctic islands. The GIA estimate compares well with results obtained from recent forward models (IJ05-R2) and inverse methods (AGE-1). The Bayesian framework is sufficiently flexible that it can, eventually, be used for the whole of Antarctica, be adapted for other ice sheets and utilise data from other sources such as ice cores, accumulation radar data, and other measurements that contain information about any of the processes that are solved for.


Introduction
Changes in mass balance of the Antarctic Ice Sheet have profound implications on sea level. While there is a general consensus that West Antarctica has experienced ice loss over the past 2 decades, the range of mass-balance estimates still differ significantly (compare, for example, estimates in Shepherd et al. (2012), Tables S8 and S11, which range from −84 ± 18 for the Gravity Recovery and Climate Experiment (GRACE) to −13 ± 39 Gt yr −1 for ICESat for the West Antarctic Ice Sheet (WAIS) and from −24 ± 35 to 123 ± 60 for the East Antarctic Ice Sheet). Reconciling these disparate estimates is an important problem. Previous studies have made use of satellite altimetry (Zwally et al., 2005), satellite gravimetry (Chen et al., 2006;Sasgen et al., 2013;Luthcke et al., 2013), or a combination of satellite and airborne data and climate model simulations (Rignot et al., 2011b) to provide estimates. In the latter case, the balance is found by deducting output ice flux from input snowfall in a technique sometimes referred to as the inputoutput method (IOM) or mass budget method.
Different approaches have different sources of error. A key error in the gravimetry-based estimates is a result of incomplete knowledge on glacial isostatic adjustment (GIA), which constitutes a significant proportion of the mass-change signal but leakage and GRACE errors are also important (Hor-wath and Dietrich, 2009). For satellite altimetry, uncertainties arise from incomplete knowledge of the temporal variability in precipitation Frezzotti et al., 2012) and the compaction rates of firn (Arthern et al., 2010;Ligtenberg et al., 2011), quantities which play a central rôle in determining the density of the observed volume change. For the IOM, the main sources of error stem from the surface mass balance (SMB) estimates used (usually obtained from a regional climate model) and uncertainties in ice discharge across the grounding line. Recent improvements in regional climate modelling have reduced the uncertainty in the SMB component but differences between estimates for the Antarctic Ice Sheet as a whole still exceed recent estimates of its mass imbalance. For example, a recent update of the commonly used regional climate model, RACMO, has resulted in a change in the integrated ice-sheet-wide SMB of about 105 Gt yr −1 (Van Wessem et al., 2014), which is larger than most recent estimates of the Ice Sheet imbalance. This change in SMB directly impacts the IOM estimate by the same amount. It is these hard-to-constrain biases in the forward models, such as the one just described, that has, in part, motivated our approach.
In an attempt to reduce the dependency on forward models, recent studies have combined altimetry and GRACE to obtain a data-driven estimate of GIA and ice loss simultaneously (Riva et al., 2009;Gunter et al., 2014). Here, we extend these earlier approaches in a number of ways. We provide a model-independent estimate not only of GIA but also of the SMB variations, firn compaction rates, and the mass loss/gain due to ice dynamics (henceforward simply referred to as ice dynamics). In doing so, we eliminate the dependency of the solution on solid-Earth and climate models. The trends for ice dynamics, SMB, GIA, and firn compaction are obtained independently through simultaneous inference in a hierarchical statistical framework (Zammit-Mangion et al., 2014). The climate and firn compaction forward models are used solely to provide prior information about the spatial smoothness of the SMB-related processes. Systematic biases in the models have, therefore, minimal impact on the solutions. In addition, we employ GPS bedrock uplift rates to further constrain the GIA signal. In future work the GPS data will also be used to constrain localised ice mass trends that cause an instantaneous elastic response of the lithosphere (Thomas and King, 2011). The statistical framework uses expert knowledge about smoothness properties of the different processes observed (i.e. their spatial and temporal variability) and provides statistically sound regional error estimates that take into account the uncertainties in the different observation techniques (Zammit-Mangion et al., 2014). The study reported here was performed as a proof of concept for a timeevolving version of the framework for the whole Antarctic Ice Sheet, which is currently under development. The timeevolving solution will use updated data sets and, as explained above, will also solve for the elastic signal in the GPS data. In addition, it will provide improved separation of the pro-cesses because of the additional information related to temporal smoothness that can be incorporated into the framework (discussed further in Sect. 5).

Data
In this section we describe the data employed, which are divided into two groups. The first group contains observational data which play a direct rôle in constraining the mass trend. These include satellite altimetry, satellite gravimetry, and GPS data (Sects. 2.1-2.3). The second group comprises auxiliary data (both observational and information extracted from geophysical models), which we use to help with the signal separation (differentiating between the different processes we solve for accounting for their spatial smoothness) (Zammit-Mangion et al., 2014). These are discussed in Sect. 2.4.

Altimetry
We make use of two altimetry data sets in this study, obtained from the Ice, Cloud and land Elevation Satellite (ICESat) and the Environment Satellite (Envisat). In this study, we used ICESat elevation rates (dh/ dt) based on release 33 data from February 2003until October 2009(Zwally et al., 2011. The data include the "86S" inter-campaign bias correction presented in Hofton et al. (2013) and the centroid Gaussian correction (Borsa et al., 2013) made available by the National Snow and Ice Data Centre. Preprocessing was carried out as described in Sørensen et al. (2011). Since ICESat tracks do not precisely overlap, a regression approach was used for trend extraction, in which both spatial slope (both acrosstrack and along-track) and temporal slope (dh/ dt) were simultaneously estimated (Howat et al., 2008;Moholdt et al., 2010). A regression was only performed if the area under consideration, typically 700 m long and a few hundred metres wide, had at least 10 points from four different tracks that span at least 1 year. Regression was carried out twice, first to detect outliers (data points which lay outside the 2σ confidence interval) and second to provide a trend estimate following outlier omission. The standard error on the regression coefficient (in this case dh/ dt), SE coef , was calculated through where e = [e i ] is the vector of residuals, n is the sample size, and x = [x i ] is the input with mean x (Yan, 2009). It should be noted that this standard error is not equivalent to the measurement error but takes into account sample size as well as the variance of both input data and residuals of the regression. Only elevation changes with an associated standard error on dh/ dt of less than 0.40 m yr −1 were considered. This threshold was selected by trial and error to avoid N. Schoen et al. : Simultaneous solution for mass trends on the West Antarctic Ice Sheet 807 a noisy spatial pattern of points that are close together and opposite in sign, usually because the regression is based on a small subset of overpasses. Data above the latitude limit of 86 • S were omitted. The remaining data were gridded on a polar-stereographic projection (central latitude 71 • S, central longitude 0 • W; origin at the South Pole) at a 1 km resolution and then averaged over a 20 km grid. The a priori error used in the modelling framework was thus the standard deviation of the trends within each 20 km grid box. The Envisat mission data began in September 2002 and ended in April 2012. Compared to laser altimetry, radar altimetry is, in general, less suited for measurements over ice for several well-known reasons: the large spatial footprint, the relatively poor performance in steeper-sloping marginal areas (Thomas et al., 2008), and the variable snowpack radar penetration (Davis, 1996). Envisat data exhibit better temporal and spatial coverage over much of the WAIS, primarily because of the instrument issues associated with ICESat that resulted in a shorter repeat cycle and less frequent operation than originally planned. We use along-track dh/ dt trends which were obtained by binning all points within a 500 m radius and then fitting a 10-parameter least-squares model in order to simultaneously correct for across-track topography, changes in snowpack properties, and dh/ dt (Flament et al., 2012). The re-trended residuals were then used to obtain linear trends over the [2003][2004][2005][2006][2007][2008][2009] ICESat period for our study. As with ICESat, the data were averaged over a 20 km grid and the standard deviation of the trends were used as the error at this scale.

GRACE
The Gravity Recovery and Climate Experiment (Tapley et al., 2005) has provided temporally continuous gravity field data since 2002. Different methods have been used to provide mass-change anomalies from the Level 1 data. Most are based on the expansion of the Earth's gravity field into spherical harmonics; however, to make the data usable for ice mass-change estimates, it is generally necessary to employ further processing methods. These include the use of averaging kernels (Velicogna et al., 2006), inverse modelling (Wouters et al., 2008;Sasgen et al., 2013), and mass concentration (mascon) approaches (Luthcke et al., 2008). Spherical harmonic solutions usually depend on filtering to remove stripes caused by correlated errors Werth et al., 2009). In this paper, we used a mascon approach , although we stress that the framework is not limited to this class of solutions. The mascon approach employed here uses the GRACE K-band inter-satellite rangerate (KBRR) data which are then binned and regularised using smoothness constraints. The fourth release (RL4) of the atmosphere/ocean model correction, which utilises the European Centre for Medium-Range Weather Forecasts atmospheric data and the Ocean Model for Circulation and Tides, was used (Dobslaw and Thomas, 2007). Some concerns with this correction have been reported (Barletta et al., 2012), but a release of the mascon data using the corrected version (Dobslaw et al., 2013) was not available for this study. Contributions to degree-one coefficients were provided using the approach by Swenson et al. (2008). Our mascon approach does not call for a replacement of C20 coefficients. We assume that GRACE does not observe SMB or ice mass changes over the floating ice shelves as they are in hydrostatic equilibrium. Hence, all observed gravity changes over the ice shelves are assumed to be caused by GIA.
Although the mascons are provided at a resolution of about 110 km, their fundamental resolution is nearer to that of the original KBRR data at about 300 km . For the statistical framework, it is important to quantify the correlation among the mascons so that it is taken into account when inferring both the processes and associated posterior uncertainties. We quantify the spatial correlation by determining an averaging model such that the diffused signal is able to loosely reconstruct the mass loss obtained using only altimetry (and assuming that all height change occurs at the density of ice). The averaging strength between mascon neighbours is also estimated during the inference (Zammit-Mangion et al., 2014). The error on the mascon rates is assumed to be a factor of the regression residuals on the trends, in a similar manner to the altimeter data (Zammit-Mangion et al., 2014). The a priori errors, after these two steps, are shown  Fig. 1, which also indicates the length scale over which we estimated the GRACE mascons to be uncorrelated.

GPS
The GPS trends used in this work were taken from Thomas and King (2011). Not all of the trends were suitable for our analysis, as the length of record did not always coincide with the 2003-2009 ICESat period. We only used stations with contemporaneous data, as well as those where we could access the original time series to confirm that the trend had remained constant, within the error bounds, for our observation period. For the northern Antarctic Peninsula, we followed the approach suggested in Thomas and King (2011) and used the pre-2003 trends, ignoring the later trend estimates, which are strongly influenced by elastic signals. All other stations were corrected for elastic rebound as in Thomas and King (2011) and subsequently assumed to be measuring GIA only (the published rates were used). A more advanced approach where the estimated ice loss is fed back into a dynamic estimate of the elastic rebound, is being implemented for a spatiotemporal extension of the Bayesian framework. The GPS data used in this study are detailed in Table 1.

Additional data sets
RACMO. Elements of the Regional Atmospheric Climate Model version 2.1 (RACMO, Lenaerts et al., 2012) were used to constrain SMB properties. Spatially varying length scales describing the spatial smoothness of precipitation patterns were obtained from the 2003-2009 SMB anomalies (with respect to the 1979-2002 mean). These ranged from 80 km in the Antarctic Peninsula to 200 km east of Pine Island Glacier. The amplitude of the anomalies, which peaked at 50 mm water equivalent in the Antarctic Peninsula, was used to provide an order of magnitude annual amplitude for expected regional SMB variability (Zammit-Mangion et al., 2014). RACMO2.1 also provides a surface density map: the mean annual density of the surface layer. This was used to translate height changes corresponding to the SMB field to mass changes.
Firn correction. We used the firn correction anomalies for 2003-2009 (with respect to the 1979-2002 mean) from a firn compaction model (Ligtenberg et al., 2011). These anomalies were used to estimate, empirically, the correlation between firn compaction rate and SMB. This relationship was then subsequently used to determine jointly the SMB and firn correction processes, subject to the constraint that firn compaction is a linear function of SMB (supported by the high correlation between the respective 2003 and 2009 trends). The methodology automatically takes into account inflated uncertainties due to confounding of these two processes since they have identical length scales (Zammit-Mangion et al., 2014).
Ice Velocities. We use surface ice velocities derived from Interferometric Synthetic Aperture Radar data (Rignot et al., 2011a). In places where no observational data were available, estimated balance velocities were used (Bamber et al., 2000). This composite velocity field was employed to help in the separation of signals due to ice dynamics versus those due to SMB (Sect. 3).

Methodology
Our framework makes use of several recent improvements in statistical modelling which can be exploited for geophysical purposes. Complete details regarding the mathematical methods employed are given in Zammit-Mangion et al. (2014), and here we provide a conceptual overview of the approach. A description of the software implementation can also be found in Zammit et al. (2015). The statistical framework hinges on the use of a hierarchical model where the hierarchy consists of three layers: the observation layer (which describes the relation of the observations to the measured fields), the process layer (which contains prior beliefs of the fields using auxiliary data sets), and the parameter layer (where prior beliefs over unknown parameters are described).
The "observation model" is the probabilistic relationship between the observed values and the height change of each of the processes. For point-wise observations, such as altimetry and GPS, the observations were assumed to be measuring the height trend at a specific location. GRACE mascons, however, were assumed to represent integrated mass change over a given area. These mass changes were translated into height changes via density assumptions: upper mantle density was fixed at 3800 kg m −3 , ice density at 917 kg m −3 , and SMB at values ranging from 350 to 600 kg m −3 . Note that we used the density map from a regional climate model to specify the density of the surface layer.
In the "process model", four fields (or latent processes) are modelled: ice dynamics, SMB, GIA, and a field which combines the processes that result in height changes but no mass changes, i.e. firn compaction and elastic rebound. We model the height changes due to these as spatial Gaussian processes, i.e. we assume that they can be fully characterised by a mean function and a covariance function. For each field we assume that the mean function is 0 (we do not use numerical models to inform the overall mean) and that the covariance function, which describes how points in space covary, is highly informed by numerical models and expert knowledge as described next. The relationship between the observations, priors, and the latent process, defined by the process model, is shown schematically in Fig. 2. Those processes that are influenced by an observation are linked by a solid arrow and it is evident that the problem is underdetermined as there are less independent observations than there are latent processes. This is why the use of priors is important to improve source separation (i.e. for partitioning elevation change between the four latent processes shown in Fig. 2). It should also be noted that SMB and firn compaction have been assumed, in this implementation of the framework, to covary a priori, as discussed later.
The practical spatial range of surface processes -this describes the distance beyond which the correlation drops to under 10 % -was estimated from RACMO2.1 as described in Sect. 2.4. This analysis revealed, for example, that locations at 100 km are virtually uncorrelated in the Antarctic Peninsula but highly correlated inland from Thwaites Glacier. Similarly, GIA was found to have a large practical range (∼ 3000 km) from an analysis of the IJ05-R1 model (although version R2 is used for comparison in the results and discussion) . the individual fields. They are useful for helping to partition a height change between the different processes that can cause that change. For example, a long wavelength variation in height that spans different basins is likely associated with SMB, whereas a localised change that shows some relationship to surface velocity is likely associated with ice dynamics (Hurkmans et al., 2014). Hence, mass loss due to ice dynamics was assumed to mostly take place in areas of faster flow (Hurkmans et al., 2014). A "soft" constraint was thus placed on elevation rates due to ice dynamics such that it is small (1 mm yr −1 ) in areas of low velocities and can be large (up to 15 m yr −1 ) for velocities greater than 10 m yr −1 . A sigmoid function was used to describe this soft constraint: where v(s) denotes the horizontal velocity at location s. For illustration of how σ vel (s) is used, an altimetry elevation trend of 10 m yr −1 in Pine Island Glacier where velocities exceed 4 km yr −1 is within the 1σ vel interval and thus classified as "probable". However, a 10 m yr −1 trend in a region east of Thwaites, where velocities are 2 m yr −1 , would lie within the 2000σ vel level and thus assumed to be a virtually impossible occurrence a priori. At Kamb Ice Stream, this assumption had to be relaxed as this area shows thickening from the shutdown of the ice stream about 150 years ago (Retzlaff and Bentley, 1993). Although the velocity of the ice is low, the thickening occurs at relatively high rates. To reflect this, we fix σ vel (s) = 2 m yr −1 in this drainage basin. In Table 2, we outline the key length scale and amplitude constraints placed on the fields that are solved for in the framework. These soft constraints should be seen as ones characterising the solution in the absence of strong evidence to anything else. They can be "violated" when the data are sufficiently informative. In the Discussion we examine the sensitivity of the solution to these constraints.
Length scales and prior soft constraints are easily defined for Gaussian processes (or Gaussian fields) which nonetheless on the other hand, are also computationally challenging to use. Gaussian fields can however be re-expressed as Gaussian-Markov random fields by recognising that Gaussian fields are in fact solutions to a class of Stochastic partial differential equations (SPDEs, Lindgren et al., 2011). Numerical methods for partial differential equations, namely finite element methods, can thus be applied to the SPDEs in order to obtain a computationally efficient formulation of a complex statistical problem (Zammit-Mangion et al., 2014). Spatially varying triangulations (meshes) are used for the different processes reflecting the assumption that, for example, ice loss is more likely to occur at smaller scales near the margins of the Ice Sheet where fast, narrow ice streams are more prevalent than in the interior. We thus use a fine mesh at the margins (25 km) and a coarse mesh in the interior for this field. GIA, however, is assumed to be smooth. This allows us to use a relatively coarse mesh for this process (∼ 100 km).
We note that our methodology differs from others in that it is not an unweighted average of estimates with markedly different errors (Shepherd et al., 2012) or a sum of corrected data sources (Riva et al., 2009) but a process-based estimate. For each of the four fields (noting that elastic rebound and firn compaction covary in this implementation), we infer a probability distribution and standard deviation for every point in space. By relating pre-inference and post-inference variances, it is possible to assess the influence of different kinds of observation at each point on the resulting fields.

Results
Inferential results are available for the four processes shown in Fig. 2 in isolation. In this section we report the results for each of the processes in turn but emphasise that these are presented to demonstrate the methodology rather than provide final estimates. This is because, as stated in Sect. 2, improvements are planned both to the framework and the data sets that we use in it. In all the examples shown, green stippling indicates where the signal is greater than the marginal standard deviation.
Ice dynamics. We obtain an ice dynamics imbalance of −86.25 ± 16.12 Gt yr −1 . The results for ice dynamics (Fig. 3) are consistent with prior knowledge of disequilibria in ice flow in the West Antarctic Ice Sheet: for example, the ice build-up in the Kamb Ice Stream catchment (Retzlaff and Bentley, 1993) and the drawdowns in the Amundsen Sea Embayment (Flament et al., 2012). The strength of the approach is apparent when focusing on the Antarctic Peninsula. Due to the relatively narrow, steep terrain and northern latitude (which affects the across track spacing of the altimetry), satellite altimeter data are sparse, while GRACE data are influenced by leakage effects, making it challenging to localise the mass sources and sinks. We find that the frame- work places ice loss maxima at the outlets of several glaciers and ice streams, which are known to have accelerated (De Angelis and Skvarca, 2003). The result is a high-resolution map of ice mass loss or gain that can be linked to specific catchments. Strong ice loss can be observed on the Northern Peninsula at the Weddell Sea shore, at the former tributaries of the Larsen B ice shelf. The maximum ice loss rate is found in the area around Sjögren Glacier at −4.7 m yr −1 . Neighbouring Röhss Glacier, on James Ross Island, has been thinning considerably since the break-up of the Prince Gustav Ice Shelf (Glasser et al., 2011;Davies et al., 2011). This is also reflected in high loss rates. Hektoria and Evans, Gregory Glacier, and glaciers on the Philippi Rise also show strong ice mass loss signals, most likely as a result of the collapse of the Larsen B ice shelf (Scambos et al., 2004;Berthier et al., 2012). Other ice loss maxima are found in the region of the Wordie Ice Shelf, Marguerite Bay, and Loubet Coast, which corroborates findings from USGS/BAS and ASTER airborne stereo imagery analyses (Kunz et al., 2012). Ice loss is also observed on King George Island, which is in agreement with recent analyses of satellite SAR data (Osmanoglu et al., 2013), and on Joinville Island.
The gap in altimeter data around the pole results in spurious estimates for that region, and the black shaded area south of 86 • in Fig. 3a is not considered here. As expected, the marginal standard deviation, or error estimate (Fig. 3b), is lowest in the interior of the WAIS, where sampling density by altimetry is high, and highest on the Peninsula, where data are sparse. Also, steep coastal areas show larger errors, reflecting the dependency of altimeter errors on slope (see Bamber et al., 2005or Brenner et al., 2007. SMB and firn compaction. We obtain an SMB imbalance of 10.57 ± 4.98 Gt yr −1 . Figure 4a shows the trend of the cumulative SMB anomalies according to RACMO 2.1, calculated with respect to the 1979-2010 mean. This approx-imately corresponds to the signal we are estimating, since we are only considering trends with respect to a steady state SMB. A cursory inspection of the anomalies we obtain ( Fig. 4b) with those from RACMO2.1 suggests relatively poor agreement. It should be noted, however, that the anomalies over the 7-year interval are on the order of a few centimeters per year and only a limited area has a statistically significant trend in our inversion (stippled regions in Fig. 4b).
There is a difference in sign between the model and our inversion for the northern Antarctic Peninsula, but again the rates we obtain are below a significant threshold and the Peninsula possesses larger uncertainties than other areas for both our framework and the regional climate model. We compare our results with ice core trends from Medley et al. (2013) who conclude that, while in phase, RACMO2.1 appears to show exaggerated inter-annual variability in the Amundsen Sea Sector. The ice core trend labelled "MEDLEY" in Fig. 4 is the mean of three cores PIG2010, THWAITES2010, and DIV2010 collected in 2010 and the location is, consequently, the mean coordinates for all the cores. The trends at the single ice cores were not listed, but there appears to be qualitative agreement with our negative trend in the area. Burgener et al. (2013) also provide new ice core records for the Amundsen Sea sector (Satellite Era Accumulation Traverse, SEAT) and Fig. 4 also shows a comparison with their data. Trends were taken over the full 2003-2009 period relative to a mean for 1980-2009. The agreement is good for three out of five cores published. Following Burgener et al. (2013), we exclude SEAT 10-4 because of the high noise level in the isotope dating and surface undulations. SEAT10-5 shows a relatively strong negative trend that we do not reproduce. SEAT-01, SEAT-03, and SEAT-06 agree well with our results at the ± cm yr −1 level. We note, however, that there is substantial short-wavelength spatial variability in SMB based on the ice core data, which is below the resolution of our framewww.the-cryosphere.net/9/805/2015/ The Cryosphere, 9, 805-819, 2015  (Fretwell et al., 2013). Stippled points denote areas in which the mean signal is larger than the marginal standard deviation.
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q work. This also suggests that a single ice core measurement should be treated with caution in this type of comparison. Height changes from firn compaction and elastic rebound are estimated together in a single field. Because they take place on similar length scales, and there is no temporal evolution in our time-invariant solution presented here, they are confounded in this study. Since firn compaction occurs at relatively large rates (cm yr −1 ), we cannot make any useful inferences about elastic rebound rates. This issue will be less critical in the time-evolving version of the framework. The modelled inverse correlation between firn compaction and SMB (Sect. 2.4) is visible in the results (Figs. 4b and 5).
GIA. We obtain a GIA rate that is equivalent to a mass trend of 12.34 ± 4.32 Gt yr −1 . It is difficult to compare this directly with other published results because the domain is not the same. We can, however, examine individual basins. The GIA vertical velocities estimated by our framework are lower than some older forward model solutions (e.g. Peltier, 2004;Ivins and James, 2005). Our results, however, agree well with a recent GRACE-derived estimate, AGE-1, which also assumes that over the ice shelves, GIA is the sole process causing observed mass change (Sasgen et al., 2013). Compared with AGE-I, our maxima in vertical uplift are shifted towards the open ocean for both of the major ice shelves (Fig. 6a). Agreement with the trends at most GPS stations is reasonable; however, the imposed smoothness constraints have a larger influence. The W06A station (Table 1), which has a strong negative trend with a large error, exacerbated by a large elastic signal, stands out. Thomas and King (2011) show that its rate does not fit with any of the GIA models used in their comparison. The signal is effectively ignored in our framework due to the large spatial scale assumed for the GIA process. Figure 6b   robust GPS data exist, the errors are substantially reduced locally. Comparing our results, basin by basin, with other recent GIA estimates including two forward models (W12a, Whitehouse et al., 2012, andIJ05-R2 Ivins et al., 2013) and two data-driven solutions (Gunter et al., 2014, andAGE-1 Sasgen et al., 2013) we find that we have, in general, best agreement with AGE-1. Basin definitions are shown in Fig. 7. Both Gunter et al. and AGE-1 rely on GRACE data. W12a, while a forward model, was adjusted to better match GPS uplift rates on the Peninsula. For the Filchner Ronne Ice Shelf (basin 1), the AGE-1 estimate (2.1 mm yr −1 ) is slightly lower than ours (2.7 mm yr −1 ), while IJ05-2 is slightly higher (3.5 mm yr −1 ). W12a (7.2 mm yr −1 ) shows more than twice our rate in this area, while Gunter et al. (4.2 mm yr −1 ) lies between IJ05-R2 and W12a. At the Ross Ice Shelf (basin 18), the agreement with AGE − 1 and IJ05-R2 (both 1.9, RATES 2.0 mm yr −1 ) is very close. Gunter et al. (3.1 mm yr −1 ) and W12a (3.4 mm yr −1 ) are slightly higher. For basin 19, again the agreement with AGE-1 and IJ05-R2 is close with RATES at 2, AGE-1 at 1.7, and IJ05-R2 at 1.9 mm yr −1 . Basin 23, which connects the ASE to the Southern Peninsula, also yields a small uplift rate (0.4 mm yr −1 ). AGE-1 (0. 5mm yr −1 ) lies within the error estimate, with IJ05-R2 (1.7 mm yr −1 ) and Gunter et al. (2.0 mm yr −1 ) just outside and W12a considerably higher at 5 mm yr −1 . On the Southern Peninsula (basin 24), agreement with AGE-1 (1.2 mm yr −1 , RATES 1.3 mm yr −1 ) is very good, but W12a is close (1.8 mm yr −1 ). Gunter et al. and IJ05 both show uplift on the Southern Peninsula but at a higher rate of 2.4 and 3.1 mm yr −1 respectively. On the Northern Peninsula, again the agreement is best with AGE-1 (0.8, RATES 0.7 mm yr −1 ) followed by IJ05-R2 (0.5 mm yr −1 ). The W12a rate is higher www.the-cryosphere.net/9/805/2015/ The Cryosphere, 9, 805-819, 2015   Table 3 instead. Our basin 25 is equal to the sum of basins 25 and 26 in ; this is given here as basin 25 for the King estimate. at 1.7 mm yr −1 . Gunter et al. is the only model that shows a negative GIA trend (−0.70 mm yr −1 ) in this region. This comparison is designed to be illustrative rather than definitive as our final GIA solution from the time-evolving framework for the whole of Antarctica will contain several improvements related to both the data and the methodology. Specifically, we will improve and extend GPS time series and coverage, use a new GRACE mascon solution, solve directly for elastic deformation, and include a spatially varying length scale.

Discussion
In Fig. 8 and Table 3 we present the basin-scale combined ice and SMB mass balance in comparison with two recent studies using GRACE (King et al., 2012;Sasgen et al., 2013). The latter study spans the ICESat period and the rates were taken from the publication. The former study, however, spans the 2002-2010 period. Basin definitions are the same as those in Sasgen et al. (2013) Table 3. Overall, we obtain good agreement with Sasgen et al. (2013). Our mean time-averaged ice loss rate of −76 ± 15 Gt yr −1 deviates by less than 1 standard deviation from the value of −87 ± 10 Gt yr −1 obtained by Sasgen et al. (2013). Agreement at the basin scale is also good. For Basin 18, our error estimates are inflated because of the pole gap in the altimetry data. The largest differences occur in basins 19, 20, and 23. For 19 and 20, agreement is very good when comparing the sums of the two adjacent basins -suggesting that leakage effects might be affecting the ability of a GRACE-only solution to fully isolate the signal to each basin. For basin 23, the altimetry -both Envisat and ICESat -show a clear positive trend in this area (ICESat: +4 Gt yr −1 ), with only very localised ice loss signals on Ferrigno ice stream. This positive trend (as opposed to a negative trend from GRACE) reduces the ice loss estimate and causes the difference between the two estimates. The strong GRACE mass loss signal for the Amundsen Sea sector leads to increased leakage in the coastal basins. The  result shows basins 23 and 21 are strongly correlated at p = 0.96. When comparing the sum over the coastal basins 21, 22, and 23, the difference between the Sasgen et al. (2013) estimate (−80 Gt yr −1 ) and ours (−74 Gt yr −1 ) reduces to 6 Gt yr −1 .
We also compare our basin-scale results to ice loss rates from . Here, the observation periods are not identical, and the GIA estimates differ. Still, there is generally good agreement at the basin scale, in particular, where their GIA estimates  lie within our error ranges (basins 18, 19); the agreement is worst where their GIA uplift rate is a multiple of ours (sum of basins 1 and 24). Overall, their ice loss rate of −118 ± 9 Gt yr −1 is significantly higher than ours.
Integrated over the domain studied, our loss estimate is smaller than other recent estimates: Shepherd et al. (2012) arrive at −97 ± 20 Gt yr −1 for WAIS over the ICESat period; while Gunter et al. (2014) obtain −105 ± 22 Gt yr −1 . With regards to Shepherd et al. (2012) and other altimetry-based results, the discrepancy is partly explained by our estimate of a negative SMB anomaly in the ASE, while RACMO2.1 gives a positive trend in this region (Fig. 4a). Methodologies employing RACMO2.1 will, thus, attribute a greater loss (for a given height change) to ice dynamics. Since these losses occur at a higher density than SMB, the inferred mass loss is greater. With regards to Gunter et al. (2014), the discrepancy arises from the different GIA rates used in the ASE. One cause for this might be the different GRACE solutions used. Our GRACE data set  is equivalent to a RL04 GRACE solution and uses the same anti-aliasing products. In Gunter et al. (2014), RL05 GRACE solutions appear to yield higher overall mass loss estimates. Preliminary comparisons of new (RL05) mascon solutions with the RL04 ones appear to show, however, little impact on the trends.
The results for SMB are more challenging to interpret because the trend, over this time period, is relatively small (a few cm yr −1 ) and below 1 standard deviation for most of the domain (Fig. 4b). There is, however, some agreement with new in situ data from ice cores (Medley et al., 2013;Burgener et al., 2013). It should be remarked that in the ASE, where we also observe an ice loss maximum, the statistical framework might have difficulty in partitioning SMB and ice dynamics. The reason for this is that the density of the SMB changes tends to be higher at the coast, with higher temperatures and melt rates. Some of the large, negative trends seen in the ASE could thus be falsely attributed to SMB. This could be remedied in principle by including more information on the spatial patterns of SMB into our framework by using, for exam-ple, a more informative prior. Also, it should be noted that the uncertainties on our SMB rates, although low on a basin scale, are comparatively high on a small spatial scale. These issues will become less critical in a time-evolving solution because ice dynamics and SMB have very different temporal frequencies: the former tends to vary smoothly in time, while the latter has relatively large high-frequency variability. This important difference in temporal smoothness will elicit significant improvement in source separation.
Methods that combine altimetry and gravimetry such as Gunter et al. (2014) and also the framework presented here are sensitive to the SMB anomaly used. We illustrate this sensitivity through a simple calculation. Let the unobserved processes on a 1 m 2 unit area be as follows: SMB amounts to 0.2 m yr −1 at 350 kg m 3 density, GIA is 1 mm yr −1 at 3500 kg m 3 , and ice loss is −1.0 m yr −1 at 917 kg m 3 . This amounts to an observed height change of −0.799 m yr −1 . The observed mass change is −897.5 kg yr −1 over the unit area. We now try to explain these signals by taking into account GRACE and altimetry but erroneously assume a SMB rate that is 10 % too high at 0.22,m yr −1 (amounting to a positive mass change of 77 kg yr 1 ). The remaining mass signal that needs to be explained by ice and GIA is now −974.5 kg yr −1 . The unexplained height change is −1.019 m. We arrive at two equations, one for height and one for mass, that can be solved by finding the intersection of the two lines (see Fig. 9). Solving the equations, we arrive at an ice mass loss rate of −1.025 m yr −1 with a high, but still plausible, GIA rate of 6 mm yr −1 . Thus, in this example, a 10 % difference in SMB can result in a GIA estimate that is markedly higher (5 mm yr −1 ) than the truth. The resulting ice mass difference would be in the range of −40 Gt yr −1 when taken over the whole of West Antarctica. Naturally, this sensitivity acts both ways, so an underestimate in SMB would result in a lower GIA and less ice loss. In this context, both GRACE filtering and the treatment of the ICESat trends also play a major rôle. As the mass loss signal in West Antarctica is relatively localised, with high rates of elevation change confined to only a few percent of the area of a basin, the inclusion or exclusion of a single (informative) altimetry data point can alter the spatial distribution of height change considerably but the overall mass trend less, as this is constrained by GRACE.
It is also worth examining the sensitivity of the solution to the prior distributions that were derived from the forward models, auxiliary data sets, such as surface ice velocity, and expert knowledge. To do this, we changed the original amplitude and length-scale constraints as detailed in Table 4. The table also lists the original mass trend (using constraints detailed in Table 2) alongside the new estimates using the revised constraints. Changes in the characteristic length scale for GIA and SMB have a rather small effect on the integrated mass trend. However, the velocity threshold that is used to determine whether the signal is likely to be associated with ice dynamics appears to have a significant effect for the three www.the-cryosphere.net/9/805/2015/ The Cryosphere, 9, 805-819, 2015 basins that comprise the Antarctic Peninsula: 23, 24, 25. This is because, for the Peninsula, observed and balance velocities are missing in a number of places. Where this is the case, they were set to 5 m yr −1 . With a 50 m yr −1 soft threshold, this means that an ice dynamics signal is extremely unlikely in all locations with a missing velocity. Improving the velocity field in this area would, therefore, reduce this sensitivity. The GIA estimates from our study agree well with a recent GRACE-based estimate (Sasgen et al., 2013) and also compare well with a recent forward model . Compared to AGE-1, the spatial pattern of our uplift maximum is shifted away from the Peninsula and towards the Ronne Ice Shelf. The spatial pattern is closer to that of W12a and ICE-5G models, with a bimodal uplift maximum centred underneath the Ronne and Ross Ice Shelves (Fig. 6a). This spatial structure is likely to have resulted from the use of GPS uplift rates, which were also used in the calibration of the most recent forward models Ivins et al., 2013). The W12a model yields slightly higher estimates for most basins but shows good agreement on the southern Antarctic Peninsula. Whitehouse et al. (2012) remark that the uplift rates using the W12 de-glaciation history -which are already substantially lower than the ICE-5G (Peltier, 2004) model rates -can be viewed as an upper bound. Separating secular and present-day viscous and elastic signals from the trends in this area remains a challenging task and will be treated in greater detail in the spatiotemporal version of our framework.
For this proof-of-concept study, our focus lies mainly on ice dynamics, SMB, and GIA estimates, neglecting to a certain extent the influence of mass-invariant height changes (due to firn compaction and elastic deformation of the bedrock). At this stage, the framework solves for a single process that combines both these processes. In this timeinvariant framework, the two are confounded and cannot be separated, as they are not distinguishable by different densities or length scales. A better approach to solve for the elastic rebound of the crust would be to integrate a dynamic estimate that depends on the ice load changes. This approach is being implemented in the spatiotemporal version of the framework. Firn compaction is currently linked with SMB through a simple correlation model (Zammit-Mangion et al., 2014). This approach could be further improved by adding a temperature dependence, along the lines of a simple firn compaction model (Helsen et al., 2008). Finally, another open question concerns the extent of present-day viscoelastic rebound in the ASE.

Conclusions
Our proof-of-concept study shows that hierarchical modelling is a powerful tool in separating ice mass balance, SMB, and GIA processes when combining satellite altimetry, GPS, and gravimetry. We demonstrate that, using only smoothness criteria derived from forward models, it can provide an accurate estimate of the different processes. A time-varying version of the framework is currently being developed, which includes a number of improvements mentioned earlier. In particular, estimation of elastic rebound in the GPS time series, and more robust partitioning of ice dynamics and SMB will provide substantial improvements in source separation, error reduction, and GIA estimation. A central advantage of the framework is that new data -which need be neither regular or gridded -can be added at any point. For example, it is possible to extend the observation period forward or back in time using data from ERS-2, Cryosat2, or any other data set that contains information about one of the processes being solved for. This could include, for example, accumulation radar data or shallow ice cores for SMB variability or additional GPS sites as they become available. Preliminary tests have shown that the inference can also be performed without GRACE data.