A new approach to estimate ice dynamic rates using satellite observations in East Antarctica

Mass balance changes of the Antarctic ice sheet are of significant interest due to its sensitivity to climatic changes and its contribution to changes in global sea level. While regional climate models 10 successfully estimate mass input due to snowfall, it remains difficult to estimate the amount of mass loss due to ice dynamic processes. It’s often been assumed that changes in ice dynamic rates only need to be considered when assessing long term ice sheet mass balance; however, two decades of satellite altimetry observations reveal that the Antarctic ice sheet changes unexpectedly and much more dynamically than previously expected. Despite available estimates on ice dynamic rates obtained from 15 radar altimetry, information about ice sheet changes due to changes in the ice dynamic are still limited, especially in East Antarctica. Without understanding ice dynamic rates it is not possible to properly assess changes in ice sheet mass balance and surface elevation, or to develop ice sheet models. In this study we investigate the possibility of estimating ice sheet changes due to ice dynamic rates by removing modelled rates of surface mass balance, firn compaction and bedrock uplift from satellite 20 altimetry and gravity observations. With similar rates of ice discharge acquired from two different satellite missions we show that it is possible to obtain an approximation of the rate of change due to ice dynamics by combining altimetry and gravity observations. Thus, surface elevation changes due to surface mass balance, firn compaction and ice dynamic rates can be modelled and correlated with observed elevation changes from satellite altimetry. 25

Abstract. Mass balance changes of the Antarctic ice sheet are of significant interest due to its sensitivity to climatic changes and its contribution to changes in global sea level. While regional climate models successfully estimate mass input due to snowfall, it remains difficult to estimate the amount of mass loss due to ice dynamic processes. It has often been assumed that changes in ice dynamic rates only need to be considered when assessing long-term ice sheet mass balance; however, 2 decades of satellite altimetry observations reveal that the Antarctic ice sheet changes unexpectedly and much more dynamically than previously expected. Despite available estimates on ice dynamic rates obtained from radar altimetry, information about ice sheet changes due to changes in the ice dynamics are still limited, especially in East Antarctica. Without understanding ice dynamic rates, it is not possible to properly assess changes in ice sheet mass balance and surface elevation or to develop ice sheet models. In this study we investigate the possibility of estimating ice sheet changes due to ice dynamic rates by removing modelled rates of surface mass balance, firn compaction, and bedrock uplift from satellite altimetry and gravity observations. With similar rates of ice discharge acquired from two different satellite missions we show that it is possible to obtain an approximation of the rate of change due to ice dynamics by combining altimetry and gravity observations. Thus, surface elevation changes due to surface mass balance, firn compaction, and ice dynamic rates can be modelled and correlated with observed elevation changes from satellite altimetry.

Introduction
Assessing and understanding ice mass balance of the Antarctic ice sheet (AIS) is challenging due to the remoteness and extensive ice cover of the continent, resulting in a sparse network of field observations to provide information about the climate, mass balance, and bedrock uplift rates. In order for an ice sheet to be in balance, the amount of ice lost due to the processes of meltwater runoff and solid ice discharge over the grounding line needs to be balanced by accumulated snowfall. If one exceeds the other, the ice sheet either gains or loses mass, resulting in a change in ice sheet mass balance (Cuffey and Paterson, 2010). The surface processes of snowfall, snowmelt and subsequent runoff, sublimation, evaporation, and snowdrift add, remove, or distribute snow and define the surface mass balance (SMB) (e.g. Van Wessem et al., 2014). Changes in SMB occur primarily in the firn layer that covers the AIS, the intermediate product between snow and ice (Ligtenberg et al., 2011). Temperature variations, overburden pressure, deformation, and repositioning of snow grains cause snow to densify until it reaches the density of glacier ice ( ∼ 917 kg m −3 ) (Herron and Langway, 1980). This results in a change in the ice sheet surface elevation without changing the mass of the ice sheet.
When thoroughly evaluated with field observations and downscaled using statistical interpolation methods, regional climate models can be used to simulate fields of SMB components, temperature, and near-surface wind speed. Ice loss rates can be obtained by combining individual estimates of accumulation, ablation, and dynamic ice loss, with the difference between mass input and mass output providing the mass balance of the ice sheet. While SMB can be taken from re-gional climate models, estimates on ice discharge are limited and difficult to obtain. The amount of ice discharge can be estimated by obtaining the product of ice velocity and ice thickness across the grounding line. Satellite radar interferometry is used to retrieve information about ice velocity rates. The ice thickness is estimated from airborne radar or, in the absence of direct observations, using surface elevation observations under the assumption that the ice is floating once it has crossed the grounding line (Rignot and Thomas, 2002;Rignot et al., 2008;Allison et al., 2009). Commonly, changes due to ice dynamics are either estimated using satellite altimetry observations (Shepherd et al., 2012;Sasgen et al., 2013) or assumed to be insignificant when studying shortterm changes (e.g. Ligtenberg et al., 2011). However, unexpected changes in ice sheet dynamics have been observed in the past decades, with some glaciers found to accelerate, while others decelerated (Rémy and Frezzotti, 2006). In general, ice dynamics are not well known and information about ice dynamic variations is limited (Rignot, 2006;Rignot et al., 2008). This becomes an issue when assessing ice mass balance and surface elevation changes, or establishing ice sheet models.
Although satellite observations help provide information about temporal and spatial changes in ice mass and ice volume, large uncertainties remain when interpreting the signals and assigning the origin of change. Ice mass balance can be measured directly from gravity observations but needs to be separated into the possible changes caused by SMB; ice dynamics; and glacial isostatic adjustment (GIA), which is the response of the lithosphere to changes in surface loading. Changes in ice sheet thickness can be obtained from altimetry observations but need to be separated into the change caused by SMB, ice dynamics, GIA, and/or firn compaction. Observed elevation changes can subsequently be converted to changes in mass by employing firn densities.
In this study we obtain an estimate of ice sheet dynamic elevation changes by combining modelled SMB rates using the Regional Atmospheric Climate MOdel (RACMO2); Gravity Recovery And Climate Experiment (GRACE); and laser altimetry observations from the Ice, Cloud, and land Elevation Satellite (ICESat). We found that the attained estimates of ice dynamic changes obtained from GRACE and ICESat are of similar magnitude. In conjunction with our estimates on our rate of change due to ice dynamics we model the rate of change of the ice surface and compare our results with direct observations taken from ICESat measurements. A study site in East Antarctica has been chosen due to the increase in mass that has been observed there by GRACE and altimetry, suggesting a thickening of the ice sheet.

Study area
The chosen study area combines Enderby Land, Kemp Land, and Mac.Robertson Land, as well as parts of Dronning Maud Ice velocity rates are plotted, sourced from the NASA MEaSUREs Program (Rignot et al., 2011;Mouginot et al., 2012), to identify glaciers and regions with dynamic ice loss.
Land and Princess Elizabeth Land (hereafter referred to as Enderby Land for simplicity). The study area is assumed to be a stable region (e.g. Rignot et al., 2008), with the ice sheet predominantly located on bedrock above sea level, making it less vulnerable to changes in ocean temperatures. The major outlet glaciers of this region are the Lambert and Mellor glaciers feeding the Amery Ice Shelf in the east, together with the smaller (∼ 3000 km 2 ) Fisher, Scylla, and American Highland glaciers. Only smaller glaciers are found along the remaining coastal region of Enderby Land, including the Shirase, Rayner, Thyer, and Robert glaciers (Fig. 1). Previous research based on the mass budget method found the ice sheet to be largely in balance across this area, possibly even slightly thickening (Rignot, 2006;Rignot et al., 2008Rignot et al., , 2013. A general positive mass trend across this region has also been recorded by gravity and altimetry observations (e.g. Shepherd et al., 2012;Sasgen et al., 2013).

Data sets and implemented models
We use observational measurements of mass variations from GRACE and surface elevation changes observed by laser altimetry using ICESat. The regional climate model RACMO2/ANT ) is used to model the trend in SMB and to force the firn compaction model. Two versions of the RACMO2 model are used here, RACMO2.1 and RACMO2.3. The SMB used throughout this paper is the sum of snowfall, evaporation/sublimation, snowdrift, and runoff. The SMB components are provided in units of kg m −2 t −1 , where t is the temporal resolution of the model.

GRACE
We use the monthly gravity field solutions CNES/GRGS RL03-v3, provided by the Groupe de Researches de Géodésie Spatiale (GRGS). The RL03 solutions have a spatial resolution of degree and order 80 (Lemoine et al., 2013) and have been chosen due to the stabilisation process that is applied to reduce noise in form of north-south striping. This is achieved by regularising the inversion for spherical harmonic coefficients (Bruinsma et al., 2010).
Temporary changes in the Earth's gravity field can be related to changes in surface mass due to the distribution of mass, as well as the elastic and viscoelastic (GIA) response of the lithosphere, the instantaneous and long-term signal to changes in surface load (Wahr et al., 1998). We obtain mass anomalies by applying the equations that relate mass changes to gravity changes (Wahr et al., 1998) to obtain the change in mass due to SMB, and due to the viscoelastic deformation, or GIA: where R is the Earth's radius; P nm are the fully normalised Legendre functions; n and m are degree and order of the spherical harmonic coefficients, respectively; θ and λ are colatitude and longitude, respectively; and C nm and S nm are the spherical harmonic coefficients, at time t, of the GRACE anomaly fields. k n and h n are the elastic Love loading numbers (e.g. Pagiatakis, 1990) and the ratio of viscoelastic Love loading numbers (Purcell et al., 2011), depending on the degree. Purcell et al. (2011) showed that this empirical approximation permitted the accurate computation of viscoelastic uplift that was independent of any particular GIA model, provided that there has been no change in load for the past 5000 years.

ICESat
Various methods are used to estimate surface elevation changes from ICESat observations, using either along-track measurements or measurements directly taken from the crossover location (e.g. Slobbe et al., 2008;Gunter et al., 2009;Pritchard et al., 2009;Sørensen et al., 2011;Ewert et al., 2012). Due to perturbations in the orbit, deviations of the repeated ground track occur, and it is necessary to determine the surface topography to correct for cross-track variations in surface elevation due to surface slope rather than changes in ice mass.
Here we use the estimated rate of change of ice sheet elevation obtained from a newly developed technique that combines both crossover and along-track observations (Hoffmann, 2016). The method allows estimation of the local surface slope using a digital elevation model that has been derived from gridded estimates of ice height at ICESat crossover points. Over a crossover grid that geographically spans all campaign crossovers of a location, a static grid was created on which heights were interpolated at the epochs of all campaigns. The estimate of the elevation change over time is made by computing a weighted least-squares regression of the height time series of each grid node and then computing a weighted mean value for all grid nodes to derive the rate of change at the crossover. This allows not only changes in height rates to be assessed at one location over time but also a digital elevation model (DEM) to be evaluated for each crossover region directly from the data. The DEM is then used to estimate the cross-track slope at the crossovers (Hoffman, 2016).
The slope estimates at the crossovers are then interpolated along-track to remove the cross-track slope from the along-track measurements. Although the elevation change estimates from along-track measurements are naturally less precise than the rate estimates at crossovers, combining both methods significantly increases the accuracy of the crosstrack slope correction applied to the along-track data (Hoffman, 2016).

RACMO2/ANT
The RACMO2/ANT regional climate model, used to obtain SMB estimates, adopts the dynamical processes from the High Resolution Limited Area Model (HIRLAM) and the physical atmospheric processes from the European Centre for Medium-range Weather Forecasts (ECMWF) (Reijmer et al., 2005) and is forced by ERA-Interim reanalysis data at the lateral boundaries (e.g. Ligtenberg et al., 2011;. The latest version, RACMO2.3 (Van Wessem et al., 2014), extends available model data from 1979-2012 (RACMO2.1) to 1979-2015 (RACMO2.3) and improves the temporal resolution from 6-hourly (RACMO2.1) to 3-hourly (RACMO2.3) (S. R. M. Ligtenberg, personal communication, 2016). The horizontal resolution is 27 km, and the vertical resolution 40 levels. Individual SMB components are provided, including snowfall, evaporation/sublimation, and snowmelt, as well as snowdrift in RACMO2.3. Over Antarctica RACMO2/ANT is coupled with a multilayer snow model, which estimates meltwater percolation, refreezing, and runoff, as well as surface albedo and snowdrift (Van Wessem et al., 2014). The update in the physical parame-ters of RACMO2.3 results in a general increase in precipitation over the grounded East Antarctic Ice Sheet -which is in good agreement with in situ observations, ice-balance velocities, and GRACE measurements -and shows a general improvement of the SMB (Van Wessem et al., 2014).

Firn compaction
We developed a firn compaction model based on the firn densification model of Ligtenberg et al. (2011), using nearsurface climate provided by RACMO2.1. It is a onedimensional, time-dependent model that estimates density and temperature individually for each layer and at each time step in a vertical firn column. The firn densification model of Ligtenberg et al. (2011) adds new snowfall instantly to the current top layer until the layer thickness exceeds ∼ 15 cm (S. R. M. Ligtenberg, personal communication, 2016), at which time it is divided into two layers. The properties of each layer are passed on to both layers. If a layer becomes too thin, due to compaction or surface melt, the layer is merged with the next layer and assigned the average properties of both layers. Our model has been simplified to improve the computational time. Rather than adding new snowfall instantly to the top layer, we compute the monthly sum of SMB and use the monthly averaged surface temperature to estimate the densification rate, density, and new temperature to obtain the vertical velocity of the surface due to monthly firn compaction.
The model starts with a new firn layer created by the total SMB of 1 month and is built up by adding a new layer each month using monthly SMB values and mean surface temperatures. The surface snow density of each top layer is estimated using the proposed parameterisation of Kaspers et al. (2004), together with a proposed slope correction to improve the fit in Antarctica by Helsen et al. (2008): where T is the average annual temperature (in K), A the average annual accumulation (in mm water equivalent (w.e.) yr −1 ), and W the average annual wind speed 10 m above the surface (in m s −1 ). The densification rate is obtained using a dry-snow densification expression proposed by Arthern et al. (2010): where C is the grain-growth constant (m s 2 kg −1 ), independently calculated for densities below (C = 0.07) and above (C = 0.03) the critical density of 550 kg m −3 ; A is the accumulation rate (mm w.e. yr −1 ); g the gravitational acceleration; and ρ and ρ i are the local density and the ice density (kg m −3 ), respectively. The exponential term includes the activation energy constants (kJ mol −1 ) for creep and for grain growth, E c and E g , respectively; the gas constant R (J mol −1 K −1 ); and the local temperature T and annual average temperature T av (K). The process of liquid water percolation and refreezing is incorporated as a function of snow porosity P s and density, as proposed by Coléou and Lesaffre (1998) (Ligtenberg et al., 2011;Kuipers Munneke et al., 2015): with the snow porosity where ρ is the density of the layer and ρ i the density of glacier ice. The heat transport throughout the firn column is solved explicitly using the one-dimensional heat-transfer equation (Cuffey and Paterson, 2010) with κ being the thermal diffusivity and z the depth. Initially the heat-transfer equation consists of a term for heat conduction, advection, and internal heating. However, initial heating is small within the firn layer and therefore neglected, and the contribution of heat advection is taken into account by the downward motion of the ice flow (Cuffey and Paterson, 2010;Ligtenberg et al., 2011). Finally, once the densification rate is estimated, the vertical velocity of the surface due to firn compaction, V fc , can be assessed by integrating over the displacement of the compacted firn layers over the length of the firn column (Helsen et al., 2008): where z is depth, ρ density, and dρ(z)/ dt the densification rate. Ligtenberg et al. (2011) found that Eq. (4) over-predicts the rate of densification for most regions in Antarctica, with the effect of the annual average accumulation being too large on the densification rate. They reintroduced an accumulation constant that previously had been proposed by Herron and Langway (1980) as α in A α (below 550 kg m −3 ) and β in A β (above 550 kg m −3 ), initially chosen between 0.5 and 1.1 but later assumed to be α, β = 1 (Zwally and Li, 2002;Helsen et al., 2008). Ligtenberg et al. (2011) applied a modelled-toobserved ratio to correct for the accumulation dependence. We also found that Eq. (4) over-predicts the rate of densification, depending on the rate of the average annual accumulation.
However, due to our use of monthly layers, the ratio proposed by Ligtenberg et al. (2011) is no longer valid and we Therefore, we continued to use our firn compaction model using RACMO2.1 near-surface climate data.
In Fig. 2a we show the average annual rate of firn compaction across the study site, and in Fig. 2b the differences between our model and the model of Ligtenberg et al. (2011). Along the ice sheet margins and the Amery Ice Sheet our model overestimates their firn compaction rates by 5-10 cm yr −1 , while it underestimates rates by 7-12 cm yr −1 in most other areas further inland, with up to 15 cm yr −1 at two individual location near 28 • E and between 68 and 70 • E. These differences are within our estimated uncertainty, based on the uncertainties provided for the modelled SMB from RACMO2.

Method to estimate the rate of change due to ice dynamics
A change in surface elevation, dH /dt, as measured by satellite altimetry is caused by a combination of processes that affect ice sheet thickness as well as the effect of GIA. The temporal change in surface height can be described as with the individual components representing elevation changes related to SMB (dH SMB /dt), firn compaction (dH fc /dt), ice dynamics (dH ice /dt), and the elastic and viscoelastic response of the lithosphere combined under the term of GIA (dH GIA /dt). While the process of firn compaction plays an important role in surface elevation changes, it does not affect the overall mass balance of the ice sheet. Therefore, the general change in ice mass as detected by GRACE can be expressed as with the individual components representing a change in mass due to SMB (dM SMB /dt), ice dynamics (dM ice /dt), and GIA (dM GIA /dt).
With the components that assemble dM SMB /dt being represented by regional climate models simulating near-surface climate in Antarctica, and dM GIA /dt modelled by available GIA models, dM ice /dt remains the only unknown in Eq. (10). Therefore, an estimate of dM ice /dt can be obtained by removing dM SMB /dt and dM GIA /dt from the GRACE observations: Similarly, the same approach can be used to obtain dH ice /dt from altimetry: The solutions to Eqs. (10) and (11)  , associated with changes in ice dynamics. We assume that changes within the firn layer have been taken into account by removing the rate of change due to SMB and firn compaction from the observations, and that the remaining signal is solely due www.the-cryosphere.net/11/1235/2017/ The Cryosphere, 11, 1235-1245, 2017 to changes within the glacier ice. Therefore, we can convert to (from) the rate of change in mass and surface elevation by dividing (multiplying) by the density of glacier ice. Thus, observations from each satellite mission can provide an independent estimate of the ice dynamics. We first correct both observational measurements, GRACE and ICESat, for GIA using three available GIA models: the W12a model of Whitehouse et al. (2012), the ICE-6G_C (VM5a) model of Peltier et al. (2015), and the recomputed version ICE6G_ANU of Purcell et al. (2016). Changes due to SMB are modelled using RACMO2.3/ANT, and the total trend due to SMB, for the period 2003-2009, is obtained using the monthly SMB (kg m −2 mth −1 ). The change in dH SMB /dt is acquired by dividing dM SMB /dt by the density of surface snow (Eq. 3), and the rate of change due to firn compaction, dH fc /dt, is taken into account by using our modelled firn compaction rates. Each month, the total SMB is computed and a monthly average firn compaction rate is removed from the SMB, before calculating the overall trend dH SMB /dt over 2003 indicates that there must be an error, which can be attributed either to errors in the data processing techniques or the inability of the models to realistically simulate surface changes due to SMB, firn compaction, and/or GIA.

Results and discussion
The chosen region is part of a vast area in East Antarctica that shows an increase in mass, suggesting that the ice sheet is growing in this region. The signal the GRACE satellites detect includes changes in mass due to accumulation, ice discharge, and GIA. In Fig. 3 we show the observed change in mass measured by GRACE. Figure 3a shows the map of the GRACE mass change signal, and Fig. 3b shows a time series for a coastal location near 67 • S, 54 • E for the entire operational period. In order to obtain the signal that is solely due The time series shows a change in gravity at a chosen location in Enderby Land (67 • S, 54 • E) over the total observational period. The green line illustrates the change, assuming the gravitational change is caused by a surface mass load, and is expressed in water equivalent (w.e.) (Eq. 1); the purple line illustrates a change due to viscoelastic deformation (GIA) (Eq. 2).
to ice mass changes the contribution of GIA needs to be removed. In Fig. 4 we show the GRACE signal corrected for GIA uplift rates using the ICE-6G_C (VM5) model by Peltier et al. (2015), W12a model by Whitehouse et al. (2012), and the recomputed version ICE6G_ANU of Purcell et al. (2016). Using ICE-6G_C (VM5) (Fig. 4a) significantly reduces the observed positive anomaly in Enderby Land, while applying W12a (Fig. 4b) and ICE6G_ANU (Fig. 4c) results in a smaller reduction of the mass anomaly, yielding a similar corrected GRACE signal. Due to the similarity between the W12a and ICE6G_ANU model the W12a model was chosen to correct the satellite observations for GIA, although the effect on the rate of change due to ice dynamics is insignificant between the models due to very small uplift rates across our study region. With the contribution of GIA removed, the signal should only comprise contributions from snowfall and ice discharge. The GIA-corrected GRACE observations suggest  (Fig. 4b).
The modelled trend in SMB and surface elevation due to SMB and firn compaction can now be removed from the GRACE and ICESat observations (Eqs. 11 and 12), to ob- by dividing (multiplying) by the density of glacier ice. We converted the rate of change of surface elevation due to the ice dynamic signal obtained from ICESat into spherical harmonics to be comparable with dH ice GRACE dt . By doing this, we represent the ice height information with the same spa-tial resolution as the mass change information and impose the same potential leakage on to the altimetry observations. The estimated rate of change due to ice dynamics is shown in Fig. 5, comparing estimates obtained using two different SMB models: RACMO2.1 and RACMO2.3.
We obtained similar rates of change due to ice dynamics by removing the modelled SMB estimates from both RACMO2 models and GIA uplift rates from GRACE and ICESat observations. When using SMB estimates from RACMO2.3, the ice dynamic estimates are significant smaller and primarily present between 30 and 60 • E with estimated rates between −0.08 and −0.13 m yr −1 obtained across the region. Using SMB estimates from RACMO2.1 yields a change due to ice dynamics of −0.08 m yr −1 and above along the entire ice sheet margin of our study region, stretching across to 75 • E. Generally, when using RACMO2.3 the SMB estimates show a smaller difference between the obtained ice dynamic estimates obtained from GRACE and ICESat, improving results across the study area. However, regions remain that exhibit differences in the obtained ice dynamic signal of up to ±0.05 m yr −1 (Fig. 5c and f). Significant changes emerge between the rate of change due to ice dynamics obtained using the former and latter RACMO2 versions, with a root mean square error, averaged over the study region, of 0.019 and 0.021 m yr −1 for RACMO2.3 and RACMO2.1, respectively.
In both dH ice ICESat dt rates a positive trend is estimated across the centre of the region. This is the result of a slightly positive elevation trend that has been recorded by ICESat observations in region D (Fig. 6b).
Finally, the total change in surface elevation is modelled, based on dH SMB /dt, dH fc /dt, dH GIA /dt, and   the general positive trend across the region is modelled, together with the positive signal near 70 • E, as well as a slight negative trend across the margin. However, the strong negative trend at the Mellor Glacier is lacking, though the region does show a slight negative trend. Although the modelled trend in surface elevation suggests similar behaviour to the altimetry observations, the signal generally appears damped compared to the ICESat observations. This is likely caused by the loss of spatial resolution through the use of degree 80 spherical harmonics (the resolution of the GRACE gravity fields) to remove the ice dynamic signal. Uncertainties are estimated for the satellite observations and models individually, and error propagation is used to obtain the uncertainty of the modelled ice dynamic estimates and modelled surface elevation changes. The uncertainty estimated for the modelled surface elevation trend varies between near zero and ∼ 6 cm yr −1 across the interior and along large parts of the ice sheet margins, and up to 12 cm yr −1 for the two locations with high SMB rates. The uncertainty of the monthly GRACE solutions are derived following the method of Wahr et al. (2006) and are ∼ 8 mm w.e. yr −1 (Fig. 7a), reducing towards the polar regions due to denser ground track coverage . The uncertainties of the ICESat observations are below 0.05 m yr −1 in the interior, where a dense network of ground-tracks exists, and between 0.15 and 0.3 m yr −1 along the ice sheet margins due to greater distances between the ground tracks and steeper slopes along the margins (Hoffmann, 2016) (Fig. 7b).
For both RACMO2 models the overall uncertainty is given as 8 % for the grounded ice sheet Van Wessem et al., 2014), resulting in an estimated uncertainty of less than 1 cm yr −1 in the interior and up to 6 cm yr −1 across the high-SMB locations proposed in Enderby Land. The firn compaction model contains several error sources. In general, the complex physics of firn densification are still not fully understood, and the density of snow and firn is not well known, introducing large uncertainties into the computations  (Sutterley et al., 2014). Error sources include the parameterisations to estimate surface snow density (Eq. 3) and the densification rate (Eq. 4), together with uncertainties within the forcing climate model RACMO2. As the firn compaction model is tuned to fit observations, it is difficult to obtain realistic uncertainty estimates. However, following the idea of Helsen et al. (2008), we obtain our error estimate for the firn compaction model by assessing the propagation of the major error sources that affect firn compaction rates. This was done by applying a bias to the accumulation (8 %) and temperature (10 K; Reijmer et al., 2005;Maris et al., 2012), as well as to the surface snow density (± 20 kg m −3 ; Helsen et al., 2008). The propagation of the errors is calculated to obtain the total uncertainty of the firn compaction model (Fig. 7c). Across most of the study site the uncertainty is estimated to be around ± 2-3 cm yr −1 . However, at the two locations with the high SMB rates the uncertainty is significantly larger and is estimated to be up to 8 cm yr −1 . Uncertainties for GIA models are not provided, as the models are tuned to fit observations and the best-fitting ice sheet history and earth rheology values (e.g. Velicogna and Wahr, 2006). However, uncertainties within our study region are small due to small uplift rates and differences between the models of < 2 mm yr −1 . Therefore, the error in the modelled GIA signals in our study region is considered negligible.
To estimate the uncertainty of the modelled ice dynamics and modelled surface elevation change, the propagation of er-rors of the particular error source is obtained (Fig. 7d and e). Depending on the incorporated satellite mission the uncertainty for the modelled rate of change due to ice dynamics is up to 6 cm yr −1 (GRACE, Fig. 7d) and up to 30 cm yr −1 (ICESat, Fig. 7e), due to the larger error of the ICESat observations. The uncertainty of the modelled elevation change is 0-12 cm yr −1 (Fig. 7f), with the greatest error source being the firn compaction model.

Conclusions
The rate of change due to ice dynamics can be estimated independently from GRACE and satellite altimetry observations through the removal of GIA signals; SMB; and, in the case of altimetry, firn compaction signals. Both approaches depend upon a separate SMB model, albeit in different ways since SMB causes a mass change in GRACE observations but a height change in altimetry observations. Therefore, any errors in the modelled SMB lead to differences in the ice dynamic estimates derived from GRACE versus altimetry. Thus, this approach provides a new and independent means of assessing the accuracy of SMB models. We showed that the differences between the old and new RACMO2 versions yield significantly different ice dynamic estimates, with RACMO2.3 producing smaller differences between the GRACE-and ICESat-derived estimates.