Mass balance of the Greenland ice sheet (2003-2008) from ICESat data The impact of interpolation, sampling and firn density

. ICESat has provided surface elevation measurements of the ice sheets since the launch in January 2003, resulting in a unique dataset for monitoring the changes of the cryosphere. Here, we present a novel method for determining the mass balance of the Greenland ice sheet, derived from ICESat altimetry data. Three different methods for deriving elevation changes from the ICESat altimetry dataset are used. This multi-method approach provides a method to assess the complexity of deriving elevation changes from this dataset. The altimetry alone can not provide an


Introduction
Different satellite based measuring techniques have been used to observe the present-day changes of the Greenland ice sheet (GrIS). Synthetic Aperture Radar (SAR) imaging reveals an acceleration of a large number of outlet glaciers in Greenland Rignot et al., 2004;Rignot and Kanagaratnam, 2006;Joughin et al., 2010). Gravity changes observed by the Gravity Recovery And Climate Experiment (GRACE) show a significant mass loss (Velicogna and Wahr, 2005;Luthcke et al., 2006;Wouters et al., 2008;Sørensen and Forsberg, 2010;Wu et al., 2010). The local elevation changes of the GrIS with significant thinning along the ice margin are revealed by laser altimetry (Slobbe et al., 2008;Howat et al., 2008;Pritchard et al., 2009).
In this study, a novel mass balance estimate of the GrIS for the period 2003-2008 is presented, derived from elevation measurements from NASA's Ice, Cloud and land Elevation Satellite (ICESat), firn compaction and surface density modelling.
Published by Copernicus Publications on behalf of the European Geosciences Union. L. S. Sørensen et al.: Greenland mass balance Different methods have been used to derive secular surface elevation change estimates dH dt of snow-or ice-covered areas from ICESat data (Fricker and Padman, 2006;Howat et al., 2008;Slobbe et al., 2008;Pritchard et al., 2009). Here we use three different methods to derive dH dt and the differences are investigated.
The total volume change of the GrIS is found by fitting a smooth surface, which covers the entire ice sheet, to the ICESat derived dH dt estimates. The conversion of the derived dH dt values to a mass change is based on various elevation change correction terms and a simple surface density model. The firn correction and the surface density models are forced by climate parameters from a regional climate model (RCM).
Other studies have linked climate models and surface mass balance models in order to estimate the mass balance of the GrIS (Li et al., 2007;van den Broeke et al., 2009;Zwally et al., 2011), but in our approach, we directly use the estimated dH dt values from ICESat to derive the total mass balance including firn dynamics, driven by the HIRHAM5 high resolution RCM (Sect. 5.2).
The first part of this paper is dedicated to the description of the ICESat data and the methods used for deriving elevation and volume changes of the GrIS (Sects. 2 to 3). The volume change estimates and their associated uncertainties are presented in Sect. 4.
In the second part of this paper, the conversion from volume to mass is described (Sects. 5 to 7). This includes the changes in the firn compaction and surface density of the GrIS. The theoretical treatment of the firn processes is presented in Sect. 5. In Sect. 6, additional elevation changes, that do not contribute to the mass balance of the ice sheet, are described and quantified. The findings from both observations and model treatment are combined to derive the total mass balance of the GrIS, which is presented in Sect. 7, along with an error analysis of the mass balance.

ICESat data
ICESat carries the Geoscience Laser Altimeter System (GLAS) instrument (Abshire et al., 2005). Technical problems with the GLAS instrument early in the mission have resulted in a significant reduction in repeated tracks and, hence, in spatial resolution. As a consequence of this and due to the inclination of the satellite, the tracks are separated by approximately 30 km in the southern part of Greenland.
The GLAS/ICESat Antarctic and Greenland Ice Sheet Altimetry Data product (GLA12) (Zwally et al., 2010) was downloaded from the National Snow and Ice Data Center. This level-2 altimetry product provides geolocated and time tagged ice sheet surface elevation estimates, with respect to the TOPEX/Poseidon reference ellipsoid. The satellite laser footprint size is 30-70 m and the distance between the footprint centres is approximately 170 m. This study is based on the 91-day repeat cycle ICESat data (release 31) from Octo- Table 1. ICESat data description. Shown is the laser campaign identifier (ID), data release number (RL), and time span of the campaigns. N and M are the number of measurements from the GrIS before and after the data culling, respectively.  Table 1.

ICESat data pre-processing
A procedure of data culling and the application of corrections is necessary to reduce some of the systematic errors in the ICESat dataset and to remove problematic measurements (Smith et al., 2005). Saturation of the waveform can induce errors in surface elevation estimates (Fricker et al., 2005). Applying the saturation correction to the relevant measurements, which are flagged in the data files, reduces these errors (NSIDC, 2010). We have also used the standard deviation of the difference between the shape of the return signal and a Gaussian functional fit (the IceSvar parameter), to evaluate the data. Large standard deviations indicate less reliable surface elevation estimates , and measurements for which the misfit is large (IceSvar ≥ 0.04 V) are rejected from the further analysis. Multiple peaks can be caused by reflections from clouds and by topography in the illuminated footprint. All measurements that contain more than one peak in the return signal are rejected from the analysis. Besides these two criteria, we have also used data quality flags and warnings given with the data to reject problematic measurements. We find that these thresholds result in an overall reduction of crossover errors, see Supplement.
Only measurements from the GrIS and the surrounding glaciers and ice caps are considered in the elevation change analysis (Csatho et al., 2009). The total number of ICESat measurements from the ice covered areas is 10 367 807. After rejecting problematic measurements in the data culling procedure, the number is reduced by approximately 13% to 9 053 639, see

Methods for deriving surface elevation changes
An observed surface elevation difference may include a seasonal signal and a secular trend, but also components which are not related to the ice sheet mass balance. The seasonal variations are caused by variations in accumulation, flow, melt and a temperature dependent firn compaction rate. The compaction of the firn, the vertical bedrock movement caused by glacial isostatic adjustment (GIA), and presentday mass changes all cause elevation changes which are part of the observed elevation difference, but do not contribute to the ice sheet mass balance. Furthermore, a potential elevation bias between the ICESat laser campaigns must be considered, since this would also be interpreted as elevation changes.
The individual ICESat tracks are not precisely repeated but can be up to several hundred metres apart. Thus, besides the previously described signals, an observed elevation difference between tracks contains a contribution from the terrain.
The fact that the ICESat measurements are not exactly repeated, complicates the methods for deriving dH dt , since any separation between two measurements introduces a surface slope component, which can be decomposed into an alongtrack and a cross-track component. Several methods for deriving dH dt from ICESat data have previously been published (Fricker and Padman, 2006;Howat et al., 2008;Slobbe et al., 2008;Pritchard et al., 2009). We present dH dt results obtained by using three different methods (M1-M3). The methods have different strengths and weaknesses, which will be discussed in the following. M1-M3 are all set up to estimate dH dt at a 500 m along-track resolution. At track crossover locations, measurements from both tracks are used to derive the dH dt values. In all three approaches, we solve for both a secular trend, dH dt and a seasonal signal, s(t). Hence, the time dependent surface elevation,H (t), is parameterized as where the seasonal signal is given by: with amplitude D = α 2 + β 2 , period T (365 days), and phase φ. Each of the dH dt estimates from the three methods are associated with a variance from the regression procedure applied. We do not perform a full analytical error propagation through the dH dt calculation. We assume that the segment size of 500 m is small enough so that the error on the measurements (and on the DEM in M1) can be assumed constant within each segment and the variances are estimated from the regression analysis. Hence, the variances reflect both the measurement error and the goodness of the fit.

Method 1
A Digital Elevation Model (DEM) can be used to correct for the surface slope and this approach is used in the first method (M1). Unfortunately there are no independent, sufficiently accurate high resolution DEMs available which cover the entire GrIS. Following Slobbe et al. (2008), we choose the DEM generated from the first seven campaigns of ICESat data (DiMarzio et al., 2007). The grid spacing of this DEM is 1 km and the elevations are given relative to the WGS84 ellipsoid.
In order to subtract the DEM from the ICESat data, the DEM is linearly interpolated to estimate the value at each data location. The height of each ICESat measurement above the reference DEM is given by: where H ICESat is translated into elevations above the WGS84 ellipsoid, to be comparable with the DEM elevations (H DEM ). The measurements are categorized according to the ICE-Sat track (i) and 500 m along-track segment denoted j . The mean of the H M1 values of each ICESat campaign is calculated in each segment, creating time series of H M1 values along-track.
where A ij = dH dt ij , B ij is the offset between the DEM and the ICESat elevations in the segment, andt is the mean time of a campaign in a given segment. The governing equation, Eq. (4) is solved using ordinary least squares regression.
Only the long wavelength component of the terrain slope is removed, due to the relative low resolution of the DEM, compared to the spacing of the ICESat along-track measurements. The 1 km resolution is too low to capture the true topography in some areas and this will most likely be reflected in the elevation changes calculated using this method. Furthermore, since the DEM used here is based on the first seven ICESat campaigns, the reference epoch will not be the same in each segment.

Method 2
In the second method (M2), data from two ICESat campaigns are used to create a reference surface, which is assumed to represent the topography in each segment. The reference surface is represented by a centroid point ( dH dx , dH dy and it is found by a least squares fit of these surface parameters to the measurements from two campaigns. The choice of the two campaigns, which are used to generate the reference surface, is based on two criteria. The first criterion is that the two campaigns are separated by one year in time. This ensures that both the seasonal signal and the actual change in elevation between the two campaigns are minimized. The second criterion is that the ICESat tracks, used to generate the reference surface, are the ones that span the largest area. These criteria help to ensure that the reference surface is representative of the surface slope. Hence, this reference surface is considered the reference for all other ICESat measurements in a given along-track segment, similar to the use of a DEM in M1: The height of the reference surface at a point (x,y) is given by: The approach of solving for dH dt is similar to Eq. (4). In spite of the criteria used to select the ICESat campaigns from which the reference surface is generated, method M2 is sensitive to seasonal variations and actual elevation changes between the two campaigns chosen. The dH dt estimates will, therefore, be biased. If the surface elevation has changed significantly, due to a change in mass balance in the time between the two tracks used, this method will underestimate the dH dt values, since some of the elevation change signal is removed.

Method 3
The third method (M3) is similar to the one presented in Howat et al. (2008) and Smith et al. (2009). In each alongtrack segment, the surface elevation H M3 is assumed to vary linearly with position (x,y), time (t) and a sine and cosine term, describing the seasonal signal: where A ij = dH dt ij , dH dx is the along-track slope, dH dy is the cross-track slope, and B ij is an estimate of the topography underlying the elevation changes. (x 0 ,y 0 ) is the centroid point of the area spanned by all the measurements in the track segment. In each segment, a least squares linear regression is performed to estimate all these parameters.
This method is sensitive to the track constellation in a segment. If the change in time dH dt is strongly correlated with the change in position (e.g. dH dx ), this method will not be able to separate the two components.

Elevation change results
The elevation changes obtained by the three methods show that there is good agreement between the patterns of elevation changes (see Fig. 1a-c). A distinct thinning of the ice sheet is generally found along the southeast and west coast, while a smaller but consistent thickening is found in the interior part of the ice sheet, which is in agreement with other altimetry studies Thomas et al., The Cryosphere, 5, 173-186, 2011 www.the-cryosphere.net/5/173/2011/ Table 2. The total mass balance of the GrIS estimated from three different methods for deriving dH dt , and different assumptions in the firn compaction model. The contributions to the total mass balance from above and below the ELA are specified, along with the mass balance above an altitude of 2000 m. Note that the mass balance below the ELA is unaffected by firn model processes and is, therefore, the same for with and without firn compaction correction firn assumptions.

Deriving volume changes
In order to estimate the total annual volume change, a smooth surface that covers the entire ice sheet is fitted through the dH dt estimates. The uncertainty of the total volume change is quantified using a bootstrap method.

Interpolation of volume changes
The dH dt estimates are interpolated onto a 5 × 5 km grid, using ordinary kriging. For all three method results, an exponential variogram model with a range of 50 km is used. The variogram model is based on all data and, for simplicity, it is assumed to be isotropic, see Supplement. The range and the choice of model are based on the experimental variogram. Due to the large number of dH dt estimates, only a local subset of points is used in the kriging procedure. Cross-validation analysis is applied to determine the sufficient number of the closest points to be used in the interpolation. In order to pass on the variances from the regression analysis, these are added to the variogram model (Pebesma, 1996). The R package gstat is used for the kriging procedure (Pebesma, 2004).
The estimated volume changes are summarised in Table 2. The estimates are of little significance without knowing their associated uncertainties. It is often difficult analytically to keep track of the error when different calculations have been performed on data and, therefore, a bootstrap method (Davison and Hinkley, 2006) is used here to quantify the uncertainty.

Bootstrapping
In the bootstrap method, data is repeatedly re-sampled to create numerous artificial datasets (Davison and Hinkley, 2006). The original dataset consists of m tracks of dH dt estimates. For each method, 1000 new bootstrapped datasets are created by randomly drawing m tracks with replacements from the tracks in the original dataset. Hence, the bootstrapped dataset will most likely contain multiple copies of some of the original tracks. Each bootstrapped dataset contains the same number of tracks as the original dataset. From each of these bootstrapped datasets a new estimate of the volume change is calculated in the same way as for the original dataset. Finally, the standard deviation of these 1000 bootstrapped volume estimates is used as an estimate of the standard deviation of the original volume change estimate (in a frequentist sense).

Volume change results
The 1000 bootstrap re-samples make up the distributions of the volume changes. For all methods, these distributions are approximately Gaussian and centred around the volume change estimate based on the original dataset (see Fig. 2). Hence, the 95% confidence interval of the volume change will be ±2σ , where σ is the standard deviation. The error estimates of the volume changes are summarized in Table 2.
There is a relatively large spread in the resulting volume changes. We believe that the volume estimate found from M2 of −187 ± 21 km 3 yr −1 is probably an under-estimation. It is likely that the reference surface, which is created in M2, contains an actual elevation change and this will result in biased dH dt values. M1 and M3 are similar, with volume change estimates of −231 ± 24 km 3 yr −1 and −239 ± 26 km 3 yr −1 , respectively.

Volume to mass conversion
In order to convert the derived elevation changes for the GrIS to mass change, the involved physical processes have to be known. Generally, the change in surface elevation can be written as whereḃ is the surface mass balance, ρ is the density of the snow or ice and w c is the vertical velocity of the surface due to change in firn compaction, in the following referred to as the firn compaction velocity. w ice is the vertical velocity of the ice matrix,ḃ m is the basal mass balance, w br is the vertical velocity of the underlying bedrock associated with glacioisostatic adjustment, u s is the horizontal ice velocity of the surface, S and u b is the horizontal velocity of the ice at the bed B (Paterson, 2002;Zwally and Li, 2002;Helsen et al., 2008;Zwally et al., 2011). A Cartesian coordinate system with a vertical axis pointing upwards is used, and we define accumulation positive and ablation negative.
As seen from Eq. (8), firn compaction and surface densities must be taken into account in order to convert the ICESat volume change to mass change. The firn responds to changes in surface temperature and precipitation and this response will not contribute to the mass balance. The firn response, the intercampaign bias and the glacio-isostatic adjustment are corrected for, before converting elevation change into mass change estimates. Based on Eq. (8), we then write the mass change as where H ICESat corrected is the elevation change from ICESat corrected for non-ice mass related processes. A distinction is made between the elevation changes caused by snow accumulation, surface melt or dynamical changes. Therefore, the densityρ value is chosen based on the physical process involved in the mass change (Thomas et al., 2006;Zwally et al., 2011). In the ablation zone, defined here as the area below the equilibrium line altitude (ELA), all elevation changes are assumed to be caused by either surface melt or dynamical changes. In the accumulation zone above the ELA, an elevation increase is assumed to be caused by an addition of snow/firn, while an elevation decrease is assumed to be caused by ice dynamics. Therefore, we define the densityρ to bẽ where ρ s is the surface density of firn including ice lenses, written as Here, r is the amount of refrozen melt water inside an annual firn layer, ρ i is the ice density (917 kg m −3 ) and ρ 0 = 625+18.7T +0.293T 2 is the temperature dependent density of new firn before formation of ice lenses (Reeh et al., 2005), T is in • C. The assumptions that define the applied density in the volume to mass conversion is adding to the uncertainty of the mass estimate, and this will be discussed in detail in Sect. 5.5. Comparing with other studies (e.g. Thomas et al., 2006) we also perform mass change calculations, using an alternative densityρ which replacesρ in Eq. (9), and which is defined aŝ Ifρ is applied, the elevation changes above the ELA, caused by dynamic mass losses, are not accounted for.

Firn compaction modelling
In order to estimate the effect of firn compaction on short time scales, a time-dependent densification model is needed. Following Reeh (2008), the time-dependent contribution to the elevation change from changes in firn compaction, is the sum of firn layer anomalies with respect to a steady state reference. The steady state reference is defined as the The Cryosphere, 5, 173-186, 2011 www.the-cryosphere.net/5/173/2011/ youngest layer in the firn column, which is unaffected by the inter-annual variability in the surface temperature and surface mass balance. The firn compaction velocity is then defined as where t 0 is the time of deposition, t 2 is the time of the addition of a new surface layer, λ(t 0 ,t) is the annual layer thickness at a time t = t 0 +t i after deposition and λ ref is the steady state reference. λ(t 0 ,t) depends on the local mass balance and is given by where τ is a time constant which, for the present study, is one month and δ is the Kronecker delta (Reeh et al., 2005). The firn density ρ f (t 0 ,t) can be derived from the Zwally and Li (2002) parametrization of the Herron and Langway (1980) densification model where ρ c is the critical firn density of 550 kg m −3 defined by Herron and Langway (1980), t c is the time it takes for the firn to reach the critical density and c is the densification constant describing the linear change in air volume in the firn, caused by the overlaying pressure (Reeh, 2008). Following Arthern et al. (2010), the densification constant c is given by a Nabarro-Herring type creep: where g is the gravity, E c and E g are the activation energies (60 kJmol −1 and 42.4 kJmol −1 , respectively), and T av is the average temperature. T is the seasonal temperature at depth z derived by surface temperature fluctuations, described by the general heat equation, where C is the specific heat capacity, K is the thermal conductivity and u,v,w are the velocities at the spatial coordinates x,y,z (Paterson, 2002). Equation (17) is solved following Schwander et al. (1997).

HIRHAM5 -forcing of the firn compaction model
The monthly mean surface temperature, runoff, snowfall and precipitation variables, that are required for the firn compaction model, are produced by the HIRHAM5 RCM (Christensen et al., 2006). The HIRHAM5 RCM is a hydrostatic RCM developed at the Danish Meteorological Institute (DMI) and is based on the HIRLAM7 dynamics (Eerola, 2006) and ECHAM5 physics (Roeckner et al., 2003). The HIRHAM5 RCM used here is an upgraded version of the HIRHAM model that has been used in several studies of accumulation and climate over Greenland. This model has been validated both with ice core data and automatic weather station data and is shown to perform well over Greenland (Dethloff et al., 2002;Kiilsholm et al., 2003;Box and Rinke, 2003;Stendel et al., 2008;Lucas-Picher, 2011). The lateral boundary condition for HIRHAM5 are taken from the European Centre for Medium Range Weather Forecast (ECMWF) ERA-Interim re-analysis (Simmons, 2007) at T255 (∼0.7 • or ∼77 km), which is a comprehensive re-analysis of the state of the atmosphere, using measurements from satellites, weather balloons and ground stations. A continuous simulation with HIRHAM5 at 0.05 • (∼5.55 km) resolution on a rotated grid is realised. The sea-surface temperature and sea-ice distribution, taken from ERA-Interim, are interpolated to the HIRHAM5 grid and prescribed to the model. The wind components, atmospheric temperature, specific humidity and surface pressure from ERA-Interim are transmitted to HIRHAM5 every six hours for each atmospheric model level of the HIRHAM5 RCM. At the lateral boundaries of the model domain, a relaxation scheme according to Davies (1976) is applied with a buffer zone of ten grid cells. The high 5.5 km horizontal resolution data are appropriate to determine the precipitation distribution over the sharp edge of the ice sheet, prominent in the ablation zone.
A comparison of the publicly available 1.5 • × 1.5 • ERA-Interim dataset and the output from the HIRHAM5 of 0.05 • , that has been forced with the same dataset, is shown in Fig. 3. It is clear that the high resolution HIRHAM5 RCM output captures the complex coastal topography of Greenland, which the low resolution forcing field cannot. The high resolution precipitation pattern impacts on the area above the ELA, where the firn compaction correction is applied and, therefore, the benefit of the high resolution forcing field is clear (see Fig. 3).

Interpolated metric grid
In order to derive the mass change of the GrIS, the area of each grid cell must be known. To ensure the equal area of each grid box, the high resolution data from the HIRHAM5 RCM is interpolated onto the equal distance 5 × 5 km grid by a nearest neighbour interpolation. The snowfall of 2008 in two different map projections is shown in Fig. 3. The interpolation onto an equal distance grid preserves the pattern of snowfall, but introduces a latitude dependent noise which is, however, only significant over the northernmost part of Greenland (e.g., at Station Nord). Despite this noise, the interpolation of the HIRHAM5 climate output on an equaldistance grid provides a good representation of the fields and it is used here to force the surface density and firn compaction models.
www.the-cryosphere.net/5/173/2011/ The Cryosphere, 5, 173-186, 2011 [m] Fig. 3. The 2008 snowfall on a scale at 0 to 2 m of water equivalent (from blue to red). (Left) Precipitation field from the ERA-Interim re-analysis, linearly interpolated from the 1.5 • × 1.5 • resolution onto an equal distance 5 km × 5 km grid. The Greenland coastline is marked in yellow, the ice boundary in red and the green diamond marks the location of Station Nord. (Middle) The precipitation field from the HIRHAM5 RCM on its original map projection, with a grid spacing of 0.05 • × 0.05 • . This projection gives a metric resolution of ∼5.5 km × 5.5 km. (Right) Nearest neighbour interpolation of the precipitation field from the HIRHAM5 RCM onto an equal distance 5 km × 5 km grid. The highly dynamic behaviour of the precipitation from the HIRHAM5 model is preserved in the transformation of the map projections.

Refreezing of melt water and formation of ice lenses
On the GrIS, 60% of the run-off given by the HIRHAM5 RCM is assumed to refreeze in the snowpack (Reeh, 1991). The accumulation is calculated as the sum of snowfall and the refrozen run-off. To simplify the following derivation of a time dependent densification model, the refrozen run-off is assumed to refreeze inside the firn layer, from which it originates, and the water is not allowed to penetrate deeper into the firn column. This assumption is a simplification. Observations from the Arctic snowpack show that melt water often penetrates through the snowpack until it reaches a hard layer, where it flows along until it refreezes or finds a crack to propagate downwards into the deeper firn (Benson, 1962;Bøggild, 2000). In order to model this behaviour (Jansson et al., 2003), the percolation depth has to be accounted for and knowledge of grain growth in water-saturated firn is required. Development of such models is outside the scope of the present study of firn compaction, where the overburden pressure is believed to be the driving force, despite the fact that melt water percolation may redistribute the load on a layer.

Results of firn compaction and density modelling
The monthly firn layer thickness is computed from Eq. (14), using the output from the HIRHAM5 RCM as forcing. To derive the firn compaction velocity (Eq. 13), a steady state reference (λ ref ) must be defined. The time span of the climate record is too short to define a robust steady state reference for the firn compaction model. Moreover, the interannual variation in temperature and precipitation will bias a chosen reference to the climate pattern which is dominating in the time span of the reference period. To avoid defining a steady state reference layer thickness, the thickness of the top firn layers is compared over the period from 2003 to 2008. The maximum number of layers, that can be evaluated at the beginning of 2003 is 169. Hence, the thickness of the top 169 layers is compared from month to month during the period 2003 to 2008 at each grid point above the ELA. The change in thickness is shown in Fig. 4a, along with the error in the linear fit in Fig. 4d. The change in the thickness of the 169 layers is a combination of changes in accumulation/surface melt and changes in the firn compaction. The change in accumulation, given in ice equivalent, for the top 169 layer thickness, is shown in Fig. 4b. By subtracting the change in the thickness of the 169 layers, in ice equivalent, from the 169 layer firn thickness, the change in air volume of the top firn is found. The rate of change in this air volume in the firn, is equivalent to the firn compaction velocity defined in Eq. (13). The approach of evaluating the relative change in air volume in each grid point above the ELA, avoids the definition of a steady state reference for the firn compaction. The resulting firn compaction velocity is the linear trend in air volume of the top 169 layers for the period from 2003 to 2008 (Fig. 4c). Figure 4c shows how the firn compaction velocity is mainly increasing in the central area of the GrIS, whereas it is decreasing in the coastal areas. This pattern shows the importance of taking the firn compaction into account, when converting the ICESat derived volume change to a change in the total mass balance of the GrIS. Depending on the assumed density of the volume changes, the firn correction corresponds to a mass loss of 18 or 36 Gt yr −1 . This corresponds to a reduction of up to 13% in the mass loss estimates when compared to the estimate from the ICESat measurements, without applying any firn compaction correction.
It is difficult to quantify the error in the firn compaction model. Further studies have to be carried out in which the modelled firn densities are compared with in situ measurements, in order to determine the error in the firn compaction velocity. The error estimate of the firn compaction correction is found here from the error in the linear fit of the interannual variability of the firn column. The 95% confidence interval is shown in the lower panel of Fig. 4. The error associated with the firn compaction velocity is most pronounced in the coastal areas near large outlet glaciers, where the forcing field from HIRHAM5 shows the largest variability. The error in the fitted firn compaction velocities will result in an error in the estimate of the total mass loss of the GrIS. The error shown in Fig. 4f is summed over each of the 5 × 5 km grid cells above the ELA, to estimate the resulting volume error. This volume is then converted into mass, using errors between 7-15 Gt yr −1 , depending on which ice or firn density is assumed. With the possible underestimation of the firn compaction, caused by only deriving the change in compaction of the top 169 layers, the higher error estimate is probably the more realistic.
The surface densityρ, as defined by Eq. (10), takes both changes into account due to ice dynamics and surface mass balance, based on a simplified assumption of the processes causing the elevation changes. Any elevation increase, above the ELA, is assumed to be caused by surface mass imbalance, while an elevation decrease is assumed to be caused by ice dynamics. Using this assumption, two processes are neglected: basal melting and the possibility that an elevation increase is caused by ice dynamics. A maximum basal melt rate of 15-20 mm yr −1 has been re-ported by Fahnestock et al. (2001) and Buchardt and Dahl-Jensen (2007). Assuming a more realistic basal melt rate of 1 mm yr −1 everywhere above the ELA, the error of neglecting basal melt corresponds to 0.9 Gt yr −1 . Acceleration of outlet glaciers is known to cause thinning of the ice sheet further inland, and conversely, a slow down of outlet glaciers could result in a build up of ice further inland. The assumptions neglect the latter possibility above the ELA. This is reasonable, because the GrIS is generally experiencing retreat and acceleration along the margins , and any changes would be most significant below the ELA. If all elevation increases above the ELA are assumed to be caused by ice instead of firn, the corresponding error is 37 Gt yr −1 . An elevation increase caused by pure ice is, however, unlikely. The Parallel Ice Sheet Model (Bueler and www.the-cryosphere.net/5/1/2011/ The Cryosphere, 5, 1-14, 2011 the surface density model, and the result is firn compaction induced errors between 7-15 Gt yr −1 , depending on which ice or firn density is assumed. With the possible underestimation of the firn compaction, caused by only deriving the change in compaction of the top 169 layers, the higher error estimate is probably the more realistic. The surface densityρ, as defined by Eq. (10), takes both changes into account due to ice dynamics and surface mass balance, based on a simplified assumption of the processes causing the elevation changes. Any elevation increase, above the ELA, is assumed to be caused by surface mass imbalance, while an elevation decrease is assumed to be caused by ice dynamics. Using this assumption, two processes are neglected: basal melting and the possibility that an elevation increase is caused by ice dynamics. A maximum basal melt rate of 15-20 mm yr −1 has been re-ported by Fahnestock et al. (2001) and Buchardt and Dahl-Jensen (2007). Assuming a more realistic basal melt rate of 1 mm yr −1 everywhere above the ELA, the error of neglecting basal melt corresponds to 0.9 Gt yr −1 . Acceleration of outlet glaciers is known to cause thinning of the ice sheet further inland, and conversely, a slow down of outlet glaciers could result in a build up of ice further inland. The assumptions neglect the latter possibility above the ELA. This is reasonable, because the GrIS is generally experiencing retreat and acceleration along the margins , and any changes would be most significant below the ELA. If all elevation increases above the ELA are assumed to be caused by ice instead of firn, the corresponding error is 37 Gt yr −1 . An elevation increase caused by pure ice is, however, unlikely. The Parallel Ice Sheet Model (Bueler and Brown, 2009;Aschwanden and Khroulev, 2009)  used to identify regions where thickening caused by ice dynamics can occur, which reduces the error of neglecting build up of ice inland to 14 Gt yr −1 .

Additional elevation change corrections
The elevation changes observed by ICESat include signals from processes which do not contribute to the mass balance of the GrIS. The most significant contribution is the firn compaction, but it is also necessary to correct for GIA, elastic uplift caused by the present-day mass changes and the ICESat intercampaign elevation biases.

Vertical bedrock movement
Elevation changes which are not related to ice volume changes are also detected by ICESat, and the estimated dH dt values must be corrected for these changes in order to determine the mass balance of the ice sheet. A bedrock movement (w br ), caused by GIA and elastic uplift from presentday mass changes, is part of the elevation changes observed by ICESat.
The GIA contribution, according to Peltier (2004), is used. It is based on the ice history model ICE-5G and the VM2 Earth model (http://pmip2.lsce.ipsl.fr/design/ice5g/). The rate of vertical motion caused by GIA is subtracted from the ICESat dH dt estimates. This correction contributes to the mass balance of the GrIS with an amount of approximately 1 Gt yr −1 .
The present-day ice sheet mass changes cause an elastic response of the bedrock (e.g., Khan et al., 2010). These vertical displacements are estimated by solving the Sea Level Equation, the fundamental equation that governs the sea level changes associated with glacial isostatic adjustment (Farrell and Clark, 1976). Since the time scale of the mass changes considered here is extremely short compared to the Maxwell relaxation time of the mantle (Spada et al., 2010), any viscoelastic effect is neglected and the ice thickness variations deduced by ICESat are spatially convolved with purely elastic loading "h" Love numbers. Sea level variations associated with melting are computed first, taking into account the elastic response of the Earth and the gravitational interaction between the ice sheets, the oceans and the mantle. Then, vertical displacements are retrieved by the surface load history over the entire surface of the Earth, associated with ice thickness variations and sea level changes. The result in Fig. 5 is obtained from a suitably modified version of the code SELEN 2.9 (Spada and Stocchi, 2007), which solves the Sea Level Equation iteratively, essentially following a variant of the pseudo-spectral method introduced by Mitrovica and Peltier (1991). A maximum harmonic degree l max = 128 is used here. Vertical displacement is computed in the reference frame with the origin in the centre of mass of the system (Earth + Load), and includes the harmonic component

ICESat intercampaign bias correction
It has been documented that there are elevation biases between the different ICESat laser campaigns. Following the method described in Gunter et al. (2009), the trend in the ICESat intercampaign bias is estimated by O. B. Andersen and T. Bondo (personal communication, 2010). The GLA15 release 31 ocean altimetry elevations are compared to a mean sea surface topography model (DNSC08). The trend is found to be 1.29 ± 0.4 cm yr −1 , when corrected for an assumed actual sea level rise of 0.3 cm yr −1 (Leuliette et al., 2004). This trend in intercampaign biases contributes with approximately 14 ± 0.4 Gt yr −1 to the mass balance.

Mass balance of the GrIS
Determining the mass change of the GrIS is a complex problem and the result depends on the type of observation and on the level of complexity of the volume to mass conversion. This may explain differences in the estimates of the total The Cryosphere, 5, 173-186, 2011 www.the-cryosphere.net/5/173/2011/ mass balance of the GrIS, which appear in the literature. To summarise the results of this study, the total mass balance estimates of the GrIS, are listed in Table 2. We have chosen to derive the mass change both with and without applying the firn compaction correction, to highlight the importance of this correction. The second key assumption of the mass loss derivation is the densityρ, from which the volume change is related to mass. The assumption, that an elevation decrease above the ELA is caused by a loss of glacial ice somewhere in the ablation area due to ice dynamics, increases the estimated mass loss of the GrIS. The total mass balances estimates (Table 2) are derived both withρ andρ. The total error of the mass balance estimate, δ dM dt , is a function of the many different contributing components to the mass balance calculation. The error of the combined volume to mass conversion is given by the following contributors. (1) The error of the ICESat derived volume change given by the bootstrap method. The methodology of volume error determination, has also been applied on the mass change estimate, resulting in an estimated error, δ ICESat , of 18-23 Gt yr −1 . (2) The error of change in firn compaction, δ Firn . (3) The error of bedrock movement, δ br . (4) The error of the intercampaign bias correction, δ IntCamp . (5) The error from neglecting basal melt in areas of corrected eleva-tion increase, δ Melt . (6) Error of the neglecting ice build up above the ELA, δ ELABuildUp . If the errors are assumed to be independent, the sum of errors is written Our estimates of the total mass balance of the GrIS are in the range −191 ± 23 to −240 ± 28 Gt yr −1 , based on the comprehensive error analysis of the ICESat derive volume changes and the theoretical treatment of the surface density and firn compaction modelling. The spatial distribution of the mass balance is seen in Fig. 6. The total mass loss, based on M3, is equivalent to a global sea level rise of 0.66 ± 0.08 mm yr −1 .
The mass loss of the major outlet glaciers is evident in Fig. 6 and the interior part of the GrIS shows little changes over the period. West of the South Greenland ice divide, the ice sheet is gaining mass, which may be caused by an increase in precipitation (cf. Fig. 4c). The most prominent area of mass increase is found in the upper area of the Storstrømmen (Bøggild et al., 1994) outlet glacier in Northeast Greenland. The ice sheet drainage basin ending in Storstrømmen is believed to originate in the central part of the GrIS near the summit area (Rignot and Kanagaratnam, 2006). Therefore, changes in Storstrømmen glacier may be caused by effects inland, or the dynamical response of the GrIS due to changes in climate. However, this has to be verified by additional studies of this area.

Discussion and conclusions
Using three different methods to derive elevation changes of the GrIS from ICESat data during the period October 2003-March 2008 reveals a consistent picture of massive ice thinning along the margin of the GrIS and a smaller elevation increase in the interior parts. The thinning is most evident along the southeast and the west coasts. An interpolation and bootstrap approach is applied, to derive a total annual volume change of snow/ice together with the corresponding uncertainties for all three methods. Volume changes of −231 ± 24 km 3 yr −1 , −187 ± 21 km 3 yr −1 and −239 ± 26 km 3 yr −1 are computed, depending on the method used. The difference in volumes obtained, confirms that mass balance studies from satellite altimetry are sensitive to the approach chosen for deriving elevation changes. From this analysis it is concluded that the volume estimate based on M2 probably is an underestimation of the volume change. Method M1 and M3 results in similar volume change estimates. Based on the uncertainty estimate alone, it is not possible to conclude which method performs the best.
In order to correct the observed elevation changes for processes not contributing to the ice sheet mass balance, we have estimated the change in firn compaction, the vertical bedrock movement caused by GIA and elastic uplift and the trend in the ICESat intercampaign elevation bias.
www.the-cryosphere.net/5/173/2011/ The Cryosphere, 5, 173-186, 2011 The largest elevation change correction, corresponding to 36 ± 7 Gt yr −1 is the firn compaction correction. The trend in the ICESat intercampaign bias is found to be −1.29 ± 0.4 cm yr −1 which corresponds to a mass gain of approximately 14 ± 0.4 Gt yr −1 . The elastic uplift of the bedrock, caused by the present-day mass changes are found to contribute with −4 to −2 Gt yr −1 to the total mass balance and the GIA correction is 1 Gt yr −1 .
The firn compaction model, besides its application shown here, can also be used to validate the RCM forcing by comparing the modelled density of the firn with in situ observations from the GrIS. However, a model comparison study for the GrIS is not within the scope of the presented work, but might be elaborated on in the future.
Modelled surface densities are used to convert the volume change into mass balance. Based on the different methods, for deriving elevation changes, we obtain mass balance estimates of the GrIS for 2003-2008 of −233 ± 27 Gt yr −1 , −191 ± 23 Gt yr −1 and −240 ± 28 Gt yr −1 , respectively.
These mass balance estimates, are in good agreement with results obtained by others. Based on GRACE data, Velicogna (2009)