Mass balance of the Greenland ice sheet - a study of ICESat data, surface density and firn compaction modelling

ICESat has provided surface elevation measurements of the ice sheets since the launch in January 2003, resulting in a unique data set 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. 5 Four di ﬀ erent methods for deriving the elevation changes from the ICESat altimetry data set are used. This multi method approach gives an understanding of the complexity associated with deriving elevation changes from the ICESat altimetry data set. The altimetry can not stand alone in estimating the mass balance of the Greenland ice sheet. We ﬁnd ﬁrn dynamics and surface densities to be important factors in de-10 riving the mass loss from remote sensing altimetry. The volume change derived from ICESat data is corrected for ﬁrn compaction, vertical bedrock movement and an in-tercampaign elevation bias in the ICESat data. Subsequently, the corrected volume change is converted into mass change by surface density modelling. The ﬁrn compaction and density models are driven by a dynamically downscaled simulation of the 15 HIRHAM5 regional climate model using ERA-Interim reanalysis lateral boundary conditions. We ﬁnd an annual mass loss of the Greenland ice sheet of 210 ± 21 Gt yr − 1 in the period from October 2003 to March 2008. This result is in good agreement with other studies of the Greenland ice sheet mass balance, based on di ﬀ erent remote sensing 20 techniques.


Introduction
Different satellite based measuring techniques have been used to observe the presentday changes of the Greenland ice sheet (GrIS).Synthetic Aperture Radar (SAR) imaging reveals an acceleration of a large number of outlet glaciers in Greenland (Abdalati et al., 2001;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).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).
We provide a novel mass balance estimate of the GrIS for the period 2003-2008, derived from elevation measurements from NASA's Ice, Cloud, and land Elevation Satellite (ICESat), firn compaction and surface density modelling.
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 four 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 mass changes is based on firn compaction and surface density modelling, 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), 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.
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 volume to mass conversion 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 dynamics involved in elevation changes without Figures contributing to the mass balance of the GrIS is presented in Sect. 5.The findings from both observations and model treatment are combined to derive the total mass balance of the GrIS in Sect.7.

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.Successive tracks are separated by approximately 30 km in the southern part of Greenland, because the GLAS instrument has been operating only 2-3 months per year.
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 centers is approximately 170 m.This study is based on the 91-day repeat cycle ICESat data (release 31) from October 2003 to March 2008.The time span and release number of the laser campaigns in the data set are listed in Table 1.

ICESat data pre-processing
A procedure of data culling and application of corrections is necessary to reduce some of the systematic errors in the ICESat data set, 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 difference between the shape of the return signal and a Gaussian fit (the IceSvar parameter), to evaluate data.Large differences indicate Introduction

Conclusions References
Tables Figures

Back Close
Full 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.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 a satisfactory size of crossover error.
Only measurements from the GrIS and the surrounding glaciers and ice caps are considered in the elevation change analysis.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.The details are listed in in Table 1.

Methods for deriving surface elevation changes
The individual ICESat tracks are not precisely repeated but can be up to several hundred meters apart.Thus the observed elevation difference between tracks contains contributions from terrain, seasonal variations and secular trends.
The fact that the ICESat tracks are not exactly repeated, complicates the methods for deriving dH dt along-track, due to the presence of a cross-track slope, caused by the topography.The cross-track slope must be determined and subtracted in order to derive the actual elevation change.Several methods for doing this 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 four different methods (M1-M4).The methods have different strengths and weaknesses, which become apparent when comparing the results.M1-M3 are along-track analysis and are all set up to estimate Introduction

Conclusions References
Tables Figures

Back Close
Full and melt during the year.In all four 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 parameterised as: where the seasonal signal is given by: with amplitude D = α 2 + β 2 , period T (365 days), and a phase φ.
Each of the dH dt estimates from the four methods are associated with a variance from the regression procedure applied.dH dt values associated with a large variance are not used in the mass balance calculation.

Method 1
In principle, a Digital Elevation Model (DEM) could 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 DEM's available which cover the entire GrIS.Following Slobbe et al. (2008), we choose the DEM generated from the first 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 WGS 84 ellipsoid.
In order to subtract the DEM from the ICESat data, the DEM is linearly interpolated to estimate the value in 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 ).Introduction

Conclusions References
Tables Figures

Back Close
Full The measurements are catagorized according to the ICESat track (i ) and 500 m along-track segment denoted j .The mean of the ∆H M1 values of each ICESat campaigns is calculated in each segment, creating time series of ∆ HM1 values along-track.
where A i j = dH dt i j , B i j is the offset between the DEM and the ICESat elevations in the segment, and t 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.

Method 2
The second method (M2) is similar to the one presented by Pritchard et al. (2009).In each of the along-track segments, a reference surface is created from elevation measurements from two ICESat campaigns.The reference surface is represented by a centroid point (x 0 ,y 0 ,H 0 ) and slopes dH dx , dH dy .The choice of the two campaigns which are used to generate the reference surface is based on two criteria.The first criterium 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 are minimized.The second criterium is 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, it is considered the Figures 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 in 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.

Method 3
The third method (M3) is similar to the one presented in Howat et al. (2008) and Smith et al. (2009).In each along-track 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 i j = dH dt i j , dH dx is the along-track slope, dH dy is the cross-track slope, and B i j is an estimate of the topography underlying the elevation changes.(x 0 ,y 0 ) is the Introduction

Conclusions References
Tables Figures

Back Close
Full centroid point of the area spanned by all of the measurements in the track segment.In each segment, a least squares linear regression is performed to estimate the elevation change.This method is sensitive to track geometry, since the method assumes that the H dependence in x,y and t is independent.For certain track constellations this will certainly not be the case.
The surface elevation at a track crossover location is found by linear interpolation of the closest points on the two tracks, located at each side of the crossover.In order to secure a fair estimate of the elevation at the crossover, a crossover is rejected if the north-south distance between the two closest points are grater than 500 m.This rejection criterium results in a subset of approximately 266 701 crossovers accepted for further analysis.
In contrary to the other three methods M1-M3, the elevation change at the crossover locations only contains the seasonal signal and the actual change in elevation.The elevation change is estimated in the crossover location of track n and m by a simple least squares linear regression.
where ∆H M4 nm contains the elevation differences between track n and m, and ∆t nm contains the corresponding time differences.A nm = dH dt nm is the estimated elevation change in the location of the crossover between track n and m, s nm (t) is the seasonal signal and B nm is the offset.
The disadvantage of this method is the poor spatial coverage of elevation change results, especially in the southern part of Greenland.Introduction

Conclusions References
Tables Figures

Back Close
Full

Elevation change results
The elevation changes obtained by the four methods show that there is a good agreement between the patterns of elevation changes (see Fig. 1a-d).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 (Abdalati et al., 2001;Thomas et al., 2008Thomas et al., , 2009;;Slobbe et al., 2008;Pritchard et al., 2009).On the more local scale, the thickening of Flade Isblink (81.4A fixed threshold of 6 m 2 for the variance associated with the fit of the regression is applied, and the number of output values from each method is an indication of how well a given method performs.The number of dH dt estimates with variance below the threshold is 264 635 for M1, 257 241 for M2, 276 717 for M3, and 4457 for M4.
This result indicates that M3 is preferable, since the largest number of accepted output values is obtained with this method.

Deriving volume changes
In order to estimate the total annual volume change, a smooth surface is fitted through the dH dt estimates, which covers the entire ice sheet.For this purpose ordinary kriging is used.The uncertainty in 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 4 method results, an exponential variogram model with a practical range of 150 km has been used.The range and the choice of model are based on the experimental Introduction

Conclusions References
Tables Figures

Back Close
Full variogram.Due to the large number of the dH dt estimates, local neighborhood kriging is used.Cross validation analysis is applied in order to determine the sufficient number of closest points to be used in the interpolation.In order to pass on the variances from the regression analysis, from which the elevation changes are determined, these have been added to the variogram model (Pebesma, 1996).The R package gstat has been used for the kriging procedure (Pebesma, 2004).
The estimated volume changes are summarized 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 to quantify the uncertainty.

Bootstrapping
Bootstrap is a resampling method (Davison and Hinkley, 2006).The basic idea of this method can be explained by the following steps.
(1) Create a resample by drawing random samples with replacements from an original data set, where it is assumed that the observations are independent.In this way a new data set is obtained with the same length as the original data set.
(2) Estimate the wanted parameter from the resample, in this case the annual volume change.
These estimates represents a distribution of the wanted parameter, from which information of the uncertainty can be obtained.
Here, the original data set is the set of dH dt estimates.For each method 1000 resamples are created, from which an error estimate can be found.For method M1, M2, and M3 a resample is made by sampling between entire tracks contrary to individual Introduction

Conclusions References
Tables Figures

Back Close
Full In method M4 the dH dt values are independent at the crossover locations, hence a resample is made by sampling between the dH dt values.

Volume change results
The 1000 bootstrap resamples make up the distributions of the volume changes.For all methods these distributions are approximately normally distributed and centered around the point estimate of the volume change (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.In order to determine which method gives the best estimate, the four methods must be reevaluated.Method M4 gives the smallest volume change estimate of −147 ± 24 km 3 yr −1 .This was expected since the density of crossovers is clearly under-represented in the southern part of Greenland (see Fig. 1) where the largest thinning is found, and many of the outlet glaciers in these regions will then not be captured correctly.We believe that the volume estimate found from M2 of −179 ± 15 km 3 yr −1 is also 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.The fact that M2 most likely dampens the signal of areas with large elevation changes, is also reflected in the relatively low standard deviation of the bootstrap procedure.The volume change results from methods M1 and M3 are similar, with volume change estimates of −225 ± 23 km 3 yr −1 and −237 ± 25 km 3 yr −1 , respectively.We find that a larger number of accepted dH dt are obtained from M3 than M1, see Sect.3.5, and that the M1 estimates are associated with larger variances than those of M3.Furthermore, it is seen in Fig. 2 that the M1 distribution is wider than the M3 distribution.
From the above argumentation it is concluded that method M3 gives the most reliable estimate of the volume change.Introduction

Conclusions References
Tables Figures

Back Close
Full Firn compaction and surface density of the ice sheet must be taken into account, in order to relate the ICESat measurements of changes in surface elevation to mass changes.The firn compaction responds dynamically to changes in surface temperature and precipitation.This dynamic response will not contribute to the mass balance of the GrIS, and therefore it is subtracted from the observed elevation change before converting it into mass change.
In general the change in surface elevation can be parameterised by where ḃ is the surface mass balance, ρ is the density of the snow or ice, w c is the vertical velocity of the surface due to the changes 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 glacio-isostatic 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).A cartesian coordinate system with a vertical axis pointing upwards is used and we define accumulation positive and ablation negative.
Changes in w ice can be neglected over short time spans (Zwally and Li, 2002)  and ρ can only be estimated from models of the firn compaction and density and w br is estimated in accordance to Sect.6.1.

Firn compaction model
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 annual firn layer anomalies with respect to a steady state reference.The steady state reference is defined as the 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 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 r(t 0 ) is the amount of refrozen melt water, ρ i is the density of ice, τ is a time constant usually equivalent to one year and δ is the Kronecker delta function (Reeh et al., 2005).The firn density ρ f (t 0 ,t) can be derived from the Zwally and Li (2002) parameterisation of the Herron and Langway (1980) densification model

Conclusions References
Tables Figures

Back Close
Full 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 due to the overlaying pressure (Reeh, 2008).The Zwally and Li parameterisation differs from the original Herron and Langway densification model by the parameterisation of c, where the Zwally and Li parameterisation is more sensitive to the temperature (Reeh, 2008).This sensitivity is important when evaluating changes in firn compaction on short time scales.
The densification constant is given by where β(T ) is a scale factor accounting for changes in grain growth with temperature T and ρ w is the density of water.K G is the rate factor for densification adjusted for grain growth Here, K 0 G is the rate factor only for grain growth, E is the activation energy and R is the gas constant (Zwally and Li, 2002;Reeh, 2008).The effect of grain growth (β) was assumed to be eight by (Zwally and Li, 2002).Later empirical studies reported a site dependency of β between seven and three (Li et al., 2003) at sites with an annual mean temperature between −30 and −22 • C. We assume β = 8, since this study covers the entire GrIS, where annual mean temperature is exceeding temperatures of −22 • C in some areas.

HIRHAM5 -forcing of the firn compaction model
The annual mean temperature at two meter above the surface, runoff, snowfall and precipitation variables, required for the firn compaction model, are produced by dynamically downscaling the European Centre for Medium-Range Weather Forecast Introduction

Conclusions References
Tables Figures

Back Close
Full (ECMWF) ERA-Interim reanalysis with the HIRHAM5 regional climate model (RCM).The HIRHAM5 RCM (Christensen et al., 2006) is a hydrostatic RCM developed at the Danish Meteorological Institute (DMI).It is based on the HIRLAM7 dynamics (Eerola, 2006) and ECHAM5 physics (Roeckner et al., 2003).The ERA-Interim reanalysis (Simmons et al., 2007), provided by the ECMWF, is a comprehensive reanalysis of the state of the atmosphere, using measurements from satellite, weather balloons and ground stations.
A continuous simulation with HIRHAM5 at 0.05 deg.(∼5.55 km) resolution on a rotated grid is realized from 1989-2008 using the ECMWF ERA-Interim at T255 (∼0.7 • or ∼77 km) as lateral boundary conditions.The sea-surface temperature and seaice distribution, taken from ERA-Interim, were interpolated to the HIRHAM5 grid and prescribed to the model.The wind components, atmospheric temperature, specific humidity and surface pressure from ERA-Interim were 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 boxes.The high 5.5 km horizontal resolution data are appropriate to determine the precipitation distribution over the sharp edge of the ice sheet, where the ablation zone is located.The dynamical downscaling with a RCM allows to simulate climate variables, which are physically consistent, for every grid cell of the domain.
A comparison of the publicly available 1.5 • ×1.5 • ERA-Interim and the HIRHAM5 dynamical downscaling are shown in Fig. 3.It is clear how the high resolution HIRHAM5 RCM run is able to account for the complex coastal topography in Greenland.The coastal precipitation patterns propagate far inland to areas above the equilibrium line altitude (ELA), where the firn compaction is applied.This pattern is not captured by the ERA-Interim (see Fig. 3) and shows the need for the high resolution RCM's input to the firn compaction modelling.Introduction

Conclusions References
Tables Figures

Back Close
Full

Interpolated metric grid
In order to derive the mass change of the GrIS the area of each grid box has to be known.To ensure 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 neighbor interpolation.The snowfall of 2008 in the two different map projections is shown in Fig. 3.It is seen that the pattern of snowfall is preserved after the grid transformation.However, the interpolation becomes noisier the greater the distance is to the equator of the original HIRHAM map projection.The noise is seen in the high precipitation area near Station Nord in the Northeastern Greenland.Despite the noise induced by the transformation of map projections, the equal distance grid gives a good approximation of the precipitation and temperature field over the GrIS produced by the HIRHAM5 model.We will use the HIRHAM5 on the equal distance grid, to force the surface density and firn compaction.

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 annual layer in the firn, from which it originates, and the water is not allowed to penetrate deeper into the firn column.This assumption is in violation with observations from the Arctic snowpack where melt water is often seen to penetrate the snowpack until it reaches a hard layer where the melt water flows along until it refreezes or finds a crack to propagate downwards into the deeper firn (Benson, 1962;Bøggild, 2000;Jansson et al., 2003).In order to be able to model this behavior, sub-annual layering of the densification model and knowledge of grain growth in water-saturated firn would be required.Both of these are outside the scope of the present study of firn compaction, where the overburden Introduction

Conclusions References
Tables Figures

Back Close
Full 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 density of the snow/ice involved in the mass change of the GrIS in Eq. ( 10), is modeled in order to derive the mass change of the GrIS from the ICESat measurements.The density is assumed to be either the density of ice or firn, depending on the location on the ice sheet.The density of the surface firn is highly dependent on the temperature during the precipitation event.
In the ablation zone, defined here for simplification as the area below the ELA, all elevation change is assumed to be caused by ice.Above the ELA, in the accumulation zone, an elevation increase is assumed to be caused by an addition of snow/firn.However, an elevation decrease is assumed to be caused by the remote removal of ice in the ablation zone as a response to ice dynamics.The surface density is then parameterised by where ρ s is the surface density of firn including ice lenses, and is given by Here, r is the amount refrozen melt water inside an annual firn layer, ρ i = 900 kg m −3 and ρ 0 is the temperature dependent density of new firn before formation of ice lenses ρ 0 = 625 + 18.7T + 0.293T 2 (18) (Reeh et al., 2005).T is the temperature given in • C. The ELA is determined using the polynomial parameterisation described by Box et al. (2004), where the ELA is given by Introduction

Conclusions References
Tables Figures

Back Close
Full a 2nd order polynomial in West Greenland and a 5th order polynomial in East Greenland as a function of the latitude.
Based on the HIRHAM5 climatology for the period 1989 to 2008, the annual firn layer thickness has been computed according to Eq. ( 12).To derive the firn compaction velocity from Eq. ( 11) a steady state reference (λ ref ) has to be defined.The time span of the climate record is too short to define a robust steady state reference for the firn compaction.Moreover, the inter-annual variation in temperature and precipitation will bias a chosen reference to the climate pattern that is dominant in the time span of the reference period.To avoid defining the steady state reference layer thickness we have chosen to compare the thickness of the top firn layers in the period from 2003 to 2008.The maximum number of layers, which can be evaluated in 2003, is 15.Hence the thickness of the top 15 layers is compared from year to year in the period 2003 to 2008 at each grid point above the ELA.The change in the thickness is seen in Fig. 4a, along with the error in the linear fit in Fig. 4d.The change in the thickness of the 15 layers is a combination of changes in accumulation/surface melt and changes in the firn compaction.The change in the accumulation given in ice equivalent for the top 15 layer thickness is seen in Fig. 4b.By subtracting the change in the thickness of the 15 layers in ice equivalent from the 15 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. ( 11).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 15 layers for period 2003 to 2008, and is depicted in Fig. 4c.The error in the linear fit is seen Fig. 4f.
In mass loss of the ice sheet with 33-67 Gt yr −1 .This is a reduction of the mass loss of up to 30%, when compared to the direct mass loss estimate from the ICESat measurements without any firn compaction correction.The error induced by the HIRHAM5 RCM in the firn compaction model is difficult to account for.Further studies have to be conducted to compare the modeled firn densities with in situ measurements before it is possible to estimate the total errors in the firn compaction velocity.Hence, the only error estimate of the firn compaction model is from the error in the linear fit of the inter-annual variability of the firn column.The 2σ are seen in the lower panel of Fig. 4. As seen in the figure the error associated with the firn compaction velocity is most pronounced in coastal areas near large outlet glaciers, where the HIRHAM5 RCM has 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 seen in Fig. 4f has been summed over each of the 5 × 5 km grid boxes above the ELA, to give the resulting volume error.This volume induced by the error in the fitted firn velocities is then converted into mass by the surface density, resulting in a firn compaction induced error between 14-30 Gt yr −1 depending on which ice or firn density is assumed.

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 glacial isostatic adjustment (GIA), elastic uplift caused by the present-day mass changes and the ICESat intercampaign elevation biases.

Conclusions References
Tables Figures

Back Close
Full

Vertical bedrock movement
Elevation changes which are not related to ice volume changes will be detected by ICESat, and these must be removed from the estimated dH dt values in order to determine the mass balance of the ice sheet.A bedrock movement (w br ) caused by GIA and elastic uplift from present-day mass changes will be a part of the elevation changes observed by ICESat.
We use a GIA contribution, according to Peltier (2004).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 removed from the ICESat dH dt estimates.We find that 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 computed 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 with 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 results in Fig. 5 are 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 center of mass of the system (Earth+Load), and Introduction

Conclusions References
Tables Figures

Back Close
Full includes the harmonic component of degree one (Greff-Lefftz and Legros, 1997).We find that the elastic uplift correction correspond to −4 to −2 Gt yr −1 , dependent on the mass loss.The elastic vertical displacement based on the results from method M3 (Sect.3.3) is shown in Fig. 5.

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 intercampain 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 ± 4 Gt yr −1 to the mass balance.

Mass balance of the GrIS
Determining the mass change of the GrIS is a complex problem with multiple solutions, depending on the type of observation and/or the level of theoretical complexity applied to solve the problem.This complexity can explain the different estimates of the total mass balance of the GrIS, which appear in the literature.To summarize the results of our studies, the total mass balance estimates of the GrIS are listed in Table 2.We have chosen to derive the mass change with and without the firn compaction correction (Table 2) are derived with and without this remote mass loss of ice.In the calculation without remote ice loss, ρ s is applied for all elevation changes above the ELA.
Our best estimate of the present total mass balance of the GrIS is −210 ± 21 Gt yr −1 based on the comprehensive error analysis of the ICESat observations and theoretical treatment of the surface density and firn compaction modelling.The spatial distribution of the mass balance is seen in Fig. 6.This mass loss is equivalent to a global sea level rise of 0.58 mm yr −1 .The uncertainty estimate on the mass change is obtained from the bootstrap procedure.Each resample is transformed into a mass change estimates according to Sect. 5, hence the 1000 resamples will make up a distribution from which the error is obtained.
The mass loss of the major outlet glaciers is evident in the figure, along with the interior part of the GrIS showing no changes over the period.The western side of the South Greenland ice divide is appearing to gain mass, which may be caused by the increasing precipitation (cf.Fig. 4c).The most prominent area of mass increase is 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 four different methods to derive elevation changes of the GrIS from ICESat data during the period 2003-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, in order to derive a total annual vol-Introduction

Conclusions References
Tables Figures

Back Close
Full ume change of snow/ice together with uncertainties for all four methods.We find volume changes of −237 ± 25 km 3 yr −1 to −147 ± 24 km 3 yr −1 depending on the method used.We conclude that method 3 is preferable, corresponding to a volume change of −237 ± 25 km 3 yr −1 .In order to correct the observed elevation changes for processes not contributing to the mass balance, we have estimated the firn compaction, vertical bedrock movement caused by GIA and elastic uplift, and the ICESat intercampaign elevation bias.
The firn compaction model is forced by the HIRHAM5 RCM, and we find this correction to be the largest and that it contributes with approximately +57 ± 14 Gt yr −1 to the total mass balance.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 ± 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 The firn compaction model can, beside its application shown here, also be used to Finally, our total mass balance result is large compared to the ICESat derived mass loss of 139±68 Gt yr −1 found by Slobbe et al. (2009), based on data from 2003 to 2007.
We believe that we have improved the application of ICESat data to estimate the total mass balance of the GrIS, by using a novel approach including firn compaction and density modelling.Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | dH dt values, since these are highly correlated along-track.
Discussion Paper | Discussion Paper | Discussion Paper | 5 Modelling firn compaction and surface densities to be constant.Thus, over the short time span of the ICESat measurement Eq. (9) can be used to express the rate of mass change of the GrIS, derived from observations of the elevation change, ḃ = dH ICESat dt − w c − w br ρ, (10) where dH ICESat dt is the estimated elevation change observed from ICESat altimetry.w c Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Fig. 4c it is seen how the firn compaction velocity is mainly increasing in the central area of the GrIS, whereas, the firn in the coastal areas is becoming more dense.This pattern shows the importance of taking the firn processes into account, when relating an observed elevation change to a change in total mass balance of the GrIS.Depending on the assumed density of the volume changes the firn correction decreases the Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | of elevation change, to highlight the importance of this correction.The second key assumption in the derivation of the mass loss is ρ, 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 by ice dynamics, enhances the estimated mass loss of the GrIS.Therefore, the total mass balances estimates Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | validate the RCM forcing, by comparing the modelled stratification of the firn with in situ observation from the GrIS.However, a model comparison study of different RCs for the GrIS has not been within the scope of the presented work, but might be elaborated in the future.Modelled surface densities are used to convert the volume change into mass balance.Based on the preferred method M3, for deriving elevation changes, we estimate a mass balance of the GrIS for 2003-2008 of −210 ± 21 Gt yr −1 .This mass loss is equivalent to a global sea level rise of 0.58 mm yr −1 .This mass balance estimate is in good agreement with results by others.Based on GRACE data, Velicogna (2009) has estimated the mass loss to be 230 ± 33 Gt yr −1 during the period 2002-2009, and Wouters et al. (2008) find a mass loss of 179 ± 25 Gt yr −1 for the years 2003-2008.van den Broeke et al. (2009) find a total mass balance of −237 ± 20 Gt yr −1 for 2003-2008, from modeled surface mass balance and observed discharge.Discussion Paper | Discussion Paper | Discussion Paper |

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 3.The 2008 snowfall on a scale at 0 to 2 m of water equivalent (from blue to red).(Left) The ERA-Interim 1.5 • × 1.5 • resolution linear interpolated onto the equal distance 5 km × 5 km grid.(Middle) The regional HIRHAM5 dynamical downscaling of the ERA-Interim.HIRHAM5 applies a rotated 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 neighbor interpolation of the HIRHAM5 onto the equal distance 5 km × 5 km grid.The highly dynamic behavior of the precipitation from the HIRHAM5 model is preserved in the transformation of the map projections.