Time-evolving mass loss of the Greenland ice sheet from satellite altimetry

Introduction Conclusions References


Introduction
The mass balance of the Greenland Ice Sheet (GrIS) has been investigated using remote sensing data in numerous studies to date, using three distinct, and largely independent, approaches: the input-output method (IOM), which takes the difference between surface mass balance (SMB) and ice discharge (Rignot and Kanagaratnam, 2006;Rignot et al., 2008b), gravimetry (e.g.Luthcke et al., 2006;Schrama and Wouters, 2011), and satellite altimetry (e.g.Sørensen et al., 2011;Zwally et al., 2011).These approaches have been compared with varying levels of success (van den Broeke et al., 2009;Sasgen et al., 2012;Shepherd et al., 2012).For the GrIS, gravimetry and the IOM were shown to agree rather well (van den Broeke et al., 2009;Sasgen et al., 2012), both in magnitude and temporal variability.Altimetry-based estimates, on the other hand, are often presented as relatively long-term averages, for example 1992-2002 (Zwally et al., 2005) or 2003-2008(Sørensen et al., 2011)).It is therefore difficult to compare temporal variability across the three methods, and mass change magnitudes from altimetry are not always consistent with IOM and gravimetric results (Shepherd et al., 2012) or with each other (Helsen et al., 2008).Gravimetric measurements are only available from 2002 onwards (the launch of the GRACE (Gravity Recovery and Climate Experiment) satellites), whereas both altimetry and IOM can potentially provide estimates from 1992 (Rignot et al., 2011;Zwally et al., 2005).For the latter two, however, few analyses spanning such long periods have been published for the GrIS.Rignot et al. (2008b) and Rignot et al. (2011) compiled time series of GrIS mass change rates using the IOM, and radar altimetry was used by Zwally et al. (2005) for 1992-2002.Khvorostovsky (2012) examined a long dH /dt time series from satellite radar altimetry (SRA) but did not convert this to mass or volume changes.These results are, thus, not always consistent or comparable with each other, and, in a recent comparison of all three methods for both ice sheets (Shepherd et al., 2012), SRA was not included for the GrIS at all.Two factors complicate the use of SRA data for determining mass trends (as opposed to volume changes).These are inadequate sampling of the largest elevation rate areas around the margins, especially in areas of steep relief (Thomas et al., 2008), and the highly variable microwave properties of the snowpack, particularly in areas experiencing surface melting (Wang et al., 2007).Here, we address both these issues, using a novel interpolation method that accounts for unsurveyed sectors (Hurkmans et al., 2012b) and two different approaches for dealing with variable snow surface properties.We use a combination of SRA data from the second European Remote-Sensing Satellite (ERS-2, 1995(ERS-2, -2003) ) and laser altimetry (2003)(2004)(2005)(2006)(2007)(2008)(2009) to obtain a time series of GrIS mass change rate and compare the results with other published estimates.
As mentioned, over outlet glaciers, which are often steep and narrow, radar altimetry typically suffers from poor coverage due to problems with tracking over steeper slopes and a large footprint with respect to outlet glacier width (Bamber, 1994;Thomas et al., 2008).Zwally et al. (2005) attempted to resolve these issues by augmenting elevation change rates (dH /dt) from SRA with observations from the Airborne Topographic Mapper (ATM) laser altimeter instrument.Recently, an alternative interpolation method was applied and validated on Jakobshavn Isbrae, Greenland's largest outlet glacier (Hurkmans et al., 2012b), which improved dH /dt estimates from SRA over the glacier's most rapidly changing area.In this study, we extend this method to the entire GrIS and validate it on other major outlet glaciers.The same approach provides a means to interpolate the data in time as well as space and, thus, to increase the temporal resolution of the dH /dt estimates to near-annual resolution.To convert volume changes to mass, we employ a firn model that includes melt and refreezing (Reeh et al., 2005;Reeh, 2008) and is forced by output from the Regional Atmospheric Climate MOdel (RACMO2) (Ettema et al., 2009), which is in turn forced at the boundaries by European Centre for Medium-range Weather Forecasts (ECMWF) re-analyses.We obtain time-evolving mass change rates for the entire GrIS, as well as for individual basins, including a comprehensive error estimate.We then compare our results to previously published estimates to investigate the consistency across approaches.

Data
As mentioned above, we use data from two spaceborne altimeters: (i) the RA-1 (Radar Altimeter) on ERS-2 and (ii) the Geoscience Laser Altimeter System (GLAS) on ICESat (Ice, Cloud, and land Elevation Satellite).ERS-2 was in orbit between 1995 and 2011, at an altitude of 780 km with a 35-day repeat cycle.In this study, data from ERS-2 have been used until 2003, when the altimeter stopped recording.We use two different SRA processing methodologies to assess the influence of different crossover sampling and variable surface microwave properties on the elevation rate estimates.This is particularly important in Greenland, where the presence of variable quantities of water within the snowpack and the heterogeneous nature of the firn both in space and time have a strong influence on the microwave properties of the surface.Markedly different approaches are used to remove this potential bias in the SRA-derived time series, which are described in detail elsewhere (Li andDavis, 2008, 2006;Khvorostovsky, 2012).For the first data set, we use monthly elevation time series that were computed at localised geographic regions (clusters) around the altimeter crossover points following the methods described in Li andDavis (2006, 2008).This approach requires exact repeat cycles, and we focus here, therefore, on the ERS-2 mission, which employed a fixed 35-day repeat cycle.In principle, the time series could be extended back to 1992, with ERS-1 data, if a different crossover scheme were used.A firstreturn retracking procedure was applied to these data, which reduces the impact of variations in volume scattering and waveform shape on dH /dt (Davis, 1997).The second SRA data set also comprises monthly averages of dH /dt for 0.1 • by 0.2 • grid cells but using a different approach for aggregating crossovers (Khvorostovsky, 2012).As a consequence the sampling in space and time of elevation changes differs from the first data set, and this can influence the volume change estimate (Sørensen et al., 2011).This second data set employs a similar retracking algorithm but, in addition, also utilises corrections, not only for variations in backscatter but also for other waveform properties, while taking into account the temporal evolution of the gradient between backscatter and dH /dt through the entire time series (Khvorostovsky, 2012).The aim, here, is not to undertake a detailed comparison of various SRA processing methods but to demonstrate that, with suitably mature methods, different approaches can provide consistent, robust results.
The second altimeter mission used in this study is ICE-Sat, which was launched in 2003 and operated until the end of 2009.Due to rapid laser degradation, the system was switched on for only around 33 days of a full 91-day repeat cycle (Abshire et al., 2005) and therefore has limited across-track resolution compared to the original mission objectives.Along-track spacing is approximately 172 m.We use all available ICESat data (release 633).To validate the interpolation algorithm, NASA's ATM was used.Since 1993, repeated flights have been performed over Greenland, generally focusing on the faster-flowing outlet glaciers (Krabill et al., 2000(Krabill et al., , 2004)).Jakobshavn Isbrae, in particular, has been frequently overflown, but also a number of other major outlet glaciers have multiple repeat flights.
Because of the different sampling in space and time between the radar and laser data, different approaches were used to extract dH /dt.For ERS-2, least-squares linear regression was used to obtain annual dH /dt from the monthly elevation time series.For any given year, data from a sliding 3-year window were used in the regression.Prior to regression, elevation measurements outside the 2σ range around the mean were discarded.Furthermore, regression was only carried out if (i) at least five data points were available, (ii) they span at least one year, and (iii) have a standard error on the resulting dH /dt of less than 0.4 m yr −1 .A problem for SRA is the slope-induced error: for a non-flat surface, the return does not come from the point directly beneath the satellite (nadir), but from the point in the radar footprint that is closest to the satellite.The horizontal difference between these points can be considerable: for a 1 • slope, it is about 14 km.This does not affect the value of a crossover difference, but its location.Using slope and aspect from Bamber et al. (2001), the ERS-2 crossover clusters were corrected for slope-induced error as described in Hurkmans et al. (2012a).
ICESat data were preprocessed using the standard quality flags supplied by NSIDC, and in addition the same set of geophysical filters that was used by Bamber et al. (2009) was applied.During the various ICESat campaigns different lasers were used, with variable power (Abshire et al., 2005), which may be the cause of inter-campaign biases that have been identified.The biases are also partly related to errors in the Gaussian fit to the return waveform Borsa et al. (2014).A bias correction for each campaign based on ocean levels was applied to elevations measured by each campaign prior to regression.For ICESat and ATM, dH /dt's were calculated along repeat tracks according to the so-called "plane method" (Howat et al., 2008;Moholdt et al., 2010): a vertically moving plane was fitted through all points within a 1 km 2 area (in this case), and a regression simultaneously solves for the slope and dH /dt.Regression was applied iteratively until the maximum residual elevation was less than 5 m, making sure there are at least 10 footprints available from at least four different tracks spanning at least 2 years.Similar to ERS-2, data from sliding 3-year windows were used for every annual dH /dt estimate.Figure 1 shows dH /dt as derived from ERS-2 and ICESat.It is important to note that dH /dt is estimated for both data sets separately.We thus assume that, by using elevation change rates, biases in estimates of absolute elevation between the data sets cancel out.The two records are merged in the space-time interpolation (Section 3).
For ice velocities, we use the mosaic that was described by Moon et al. (2012).They used a combination of interferometric synthetic aperture radar (InSAR) and speckle tracking methods (Joughin, 2002), derived from RADARSAT-1, TerraSAR-X, and Advanced Land Observation Satellite (ALOS) images.The mosaic is based on velocity maps for the winters 2000-2001, and 2005-2006 through 2010-2011.Lastly, the delineation of the major outlet glaciers is taken from Rignot and Kanagaratnam (2006).

Interpolation
The interpolation method we employ is space-time Kriging with external drift (ST-KED), which was used and validated for Jakobshavn Isbrae in a previous study (Hurkmans et al., 2012b).There are two main differences between ST-KED and the more "standard", ordinary Kriging (OK), sometimes referred to as optimal interpolation.First, ST-KED interpolates in both space and time.Spatial and temporal data characteristics are used by fitting a semi-variogram in space and one in time, and combining them using the product-sum method (De Cesare et al., 2001;Gething, 2006).Second, ST-KED, as implemented here, uses the spatial gradient of velocity to constrain the interpolation (the external drift component) in places where data coverage is poor.The method assumes a linear relationship between elevation change rate and velocity, although the coefficients of the relation are implicitly solved for by the Kriging equations (Deutsch and Journel, 1992).A given gradient in velocity can thus produce different gradients in dH /dt, as constrained by the available altimetry data.For Jakobshavn Isbrae, the relationship had a Pearson correlation coefficient (r) of about 0.9, and ST-KED yielded more realistic dH /dt patterns and rates when compared to ATM (Hurkmans et al., 2012b).Prior to interpolation, the relationship was verified on a selection of other major outlet glaciers.Figure 2 shows scatter plots of velocity versus dH /dt (from all ICESat data, 2003ICESat data, -2009) ) and results of linear regressions between velocity and dH /dt for seven major drainage basins: Petermann (indicated by (P) in the Fig. 2), Nioghalvfjerdsbrae (N), Storstrømmen (S), Kangerdlugssuaq (K), Helheim (H), Jakobshavn (J), and Upernavik (U).The slope and r for all drainage basins are listed in Table 1 for both the ERS-2 period (1995ERS-2 period ( -2003) ) and the ICESat period (2003)(2004)(2005)(2006)(2007)(2008)(2009).
The inset of Fig. 2a shows the relationship between velocity and dH /dt for the entire GrIS.It suggests a bimodal distribution: one cloud is more or less horizontal, representing areas in which no (detectable) dynamic thinning is taking place, and another one that indeed suggests a near-linear relationship.The plot is rather noisy, most probably caused by the superposition of SMB variability onto the dH /dt trends.A similar scatter plot for some selected basins is shown in Fig. 2a, more clearly showing the bimodal behaviour.Basins (J), (H), and (K) -and, to a smaller degree, (U) -have high (strongly negative) regression slopes and correlation coefficients, indicating a dominant dynamic signal, whereas (P) and (N) show insignificant slopes and low r.(S) is a special case as it is known to be stagnated and dynamically www.the-cryosphere.net/8/1725/2014/The Cryosphere, 8, 1725-1740, 2014  (Rignot et al., 2008b) and the entire GrIS.N is the number of valid data points, r the Pearson correlation coefficient, and A the regression slope.The first three columns are for the ERS-2 data and the second three columns for ICESat.thickening (Mohr et al., 1998).This is confirmed by the scatter plot and the positive regression slope.The slope of the linear regression is, thus, an easily obtained indicator for whether or not a glacier is thinning dynamically.Figure 2b shows the slope for all basins, with the associated r plotted in Fig. 2c.The values of both are listed in Table 1.Again, Storstrømmen clearly stands out with its positive slope.SMB variability, the size of the drainage basin, gaps in velocity coverage (especially in the south), and poor dH /dt sampling all add noise to the result, mainly affecting the smaller basins.
In the interpolation, the external drift only plays a role in areas with sparse dH /dt data, or strong velocity gradients.In other areas, results from OK and KED are similar.In the interior, however, the velocity map shows a relatively low signal-to-noise ratio.For this reason, we decided to use external drift only in areas with ice velocity higher than 70 m yr −1 , and for a regression slope lower than −2.Other areas are interpolated using OK, with the same Kriging parameters.To assess the performance of the interpolation algorithm over the main outlet glaciers, dH /dt's along their centrelines are plotted in Fig. 3.The same basins (except Storstrømmen) are shown as in Fig. 2a.
Figure 3 shows interpolation results averaged for the ICE-Sat period (2003)(2004)(2005)(2006)(2007)(2008)(2009), for six selected drainage basins.For these basins, interpolated dH /dt along the glacier centreline is shown along with velocity.For comparison, four interpo-lation methods (OK, KED, and their spatiotemporal counterparts: ST-OK and ST-KED) are shown.Without velocity to aid the interpolation, the maximum thinning rate is essentially the highest value in the altimetry data set, which can be far from the grounding line, especially for ERS-2.For dynamically thinning glaciers, therefore, elevation change rates are much more realistic after interpolation with KED/ST-KED, as is confirmed by ATM measurements closer to the grounding lines of Jakobshavn Isbrae and, to a lesser extent, Helheim glacier.Because the spatiotemporal interpolation uses more data, noise that is occasionally present in OK and KED is generally smoothed out in the spatiotemporal interpolation results.

Firn density modelling
In order to convert measured dH /dt to mass change rates, several processes have to be taken into account.Firn compaction, the compression of older snow by subsequent accumulation, affects the volume but not the mass and therefore needs to be corrected for.In addition, changes in mass can be caused by variability in SMB, ice dynamics, or a combination of the two, and the effective density required to convert the volume change to mass change is, therefore, generally some (time-varying) combination of the surface and ice density.Firn compaction and the surface density are obtained from a model based on annual values of climate variables (Reeh, 2008;Reeh et al., 2005).The firn model is forced by output from the regional climate model RACMO2 (Ettema et al., 2009).RACMO2 was run at 11 km spatial resolution for the period 1958-2010, using lateral boundary conditions from ERA-40 re-analysis until 1989, and ERA-Interim after that.When validated with surface mass balance observations, the model was found to simulate SMB well (van den Broeke et al., 2009;Ettema et al., 2009Ettema et al., , 2010)).
The original model by Reeh (2008) only uses annual temperature and accumulation to force a degree-day model (Reeh, 1991) to calculate snow melt, a part of which (60 % of the annual accumulation) refreezes as a layer of superimposed ice remaining (SIR) at the end of the melt season.Every annual layer is thus composed of a fraction of SIR and a fraction of firn.Meltwater is assumed to not percolate deeper into the firn profile.However, because RACMO2 estimates snow melt and refreezing, annual SMB can be calculated directly from these data, negating the need for the degree-day model.In this way, SMB is calculated as b = S − M + R, where S is annual accumulation, M snow melt, and R refreezing, all in units kg m −2 yr −1 .SIR is then equal to refreezing, or SMB, whichever is lower.For each year for which elevation change rates are available, the ice sheet is then divided into three zones (Reeh, 2008): (1) the accumulation zone, where SMB > SIR ≥ 0; (2) the superimposed ice zone, where SMB = SIR ≥ 0; and (3) the ablation zone, where SMB < 0. In the ablation zone, measured elevation change rates are assumed to be caused by ice loss (gain), so in both of these cases the density of ice is assumed.In the accumulation zone, the surface layer consists of a SIR frac-tion with the density of ice and a firn fraction from which the surface density ρ s is calculated using an empirical relation The Cryosphere, 8, 1725-1740, 2014 www.the-cryosphere.net/8/1725/2014/based on firn temperature (Reeh et al., 2005): where T f is the firn temperature at 10 m depth in • C, which depends on SIR (units: metres ice yr −1 ) and the annual mean temperature TMA ( • C): With time, the thickness of the firn fraction of each annual layer decreases (and its density increases) due to firn compaction.The process is described using the formulation of Herron and Langway, Jr. (1980), where a rate parameter c determines the compaction speed, depending on accumulation and temperature.c is based on measured depth-density profiles in Greenland and Antarctica (Herron and Langway, Jr., 1980), where steady-state conditions are assumed.Zwally and Li (2002) base c on a relationship between temperature and activation energy fitted to measurements for Greenland, which they implemented in a time-dependent model, yielding results in the dry-snow and upper percolation zones that are consistent with observations (Zwally and Li, 2002;Reeh, 2008).In our model, therefore, we use the Zwally and Li (2002) parameterisation for c, which can be written as where b is the annual mass balance (m we yr −1 ), and β(T ) and K 0G are empirical factors related to the processes of densification and grain growth.T is absolute temperature, E(T ) is the activation energy, and R is the gas constant (8.314).K OG (T ) exp(−E(T )/(RT )) is parameterised as 8.36(273.15− T ) −2.061 (Reeh, 2008), and β(T ) is β = 139.21− 0.542T (Zwally and Li, 2002).The total thickness and density of each annual layer as it lowers down the firn profile is then easily calculated from the two fractions (Reeh, 2008).
To calculate the dH /dt due to compaction and accumulation variability, first anomalies with respect to a reference period (in which the ice sheet is assumed to be close to balance, here 1961-1990) are calculated, because in steady state the surface elevation does not change, whereas accumulation and firn compaction still occur.A firn profile of 100 annual layers is assumed, and the firn compaction velocity of a given year is calculated as the combined change in thickness of the lower 99 layers.Since RACMO2 data are available from 1958, and our dH /dt analysis started in 1995, at least 37 years of modelled layers are available.The profile is then completed with annual layers that have thicknesses and densities according to the "reference profile" .The bulk of the thickness changes occurs in the upper part of the profile, and this is, therefore, a reasonable approach.The anomaly is then this firn compaction velocity minus that of a complete reference profile.Similarly, the anomaly of the surface layer thickness that is deposited during that given year represents the surface elevation anomaly due to accumulation variability.
In order to obtain corrections that are consistent with the altimetry-derived dH /dt, which is based on linear regression of 3-year periods (see Sect. 2), a similar regression is conducted using the cumulative anomalies of SMB, surface layer thickness, and firn compaction velocity.From this regression, dH /dt due to firn compaction ( Ḣfc , where Ḣ is equivalent to dH /dt) and accumulation variability ( Ḣsmb ) are derived, as well as the Ṁ component due to SMB.
Figures 4 and 5 show these trends over the periods 2003-2009 and 1995-2002, respectively, along with the average surface density ρ s .Ḣsmb and Ḣfc will generally show opposite patterns because of our definition of positive and negative: a positive anomaly in accumulation will cause a thick layer of snow, which will compact relatively quickly immediately afterwards, causing a negative (because the surface lowers) Ḣfc .In Figure 4a and b, some areas appear to have slightly positive SMB but negative Ḣsmb anomalies.This is mainly in the percolation zone, where most, or all, of the annual accumulation melts and refreezes.SMB is not affected, but the replacement of snow by ice causes surface lowering.
To estimate the errors in the modelled dH /dt components, we start from the accuracy estimates for SMB and temperature as produced by RACMO2, and propagate the uncertainty through the model by running the model with high (plus error) and low (minus error) estimates of SMB and temperature.We chose a value for the random temperature error (root mean square error, RMSE) of 2.5 • C, which is slightly higher than the estimate by Ettema et al. (2010) over both ice sheet and land area.The accuracy of SMB, σ smb , as simulated by RACMO2 is (Ettema et al., 2009) σ smb = 15 + 0.01SMB + 0.0002SMB 2 . (4) In the accumulation zone, the maximum for σ smb is 30 % of the SMB.In the ablation zone, the quoted value is 20 % of SMB.We assumed the same value (20 %) for the SIR content in each annual layer.The errors in Ḣfc , Ḣsmb , and ρ s , given by σ H fc , σ H smb , and σ ρ , respectively, are obtained from the model runs with perturbed SMB and/or temperature and shown in Figure 6.

Conversion of volume to mass and error analysis
The total mass change rate, dM/dt, of a given grid cell with area A can be calculated by Here, SMB is expressed units of kg m −2 yr −1 , Ḣalt is the observed dH /dt (m yr −1 ), and ρ i is the density of ice (917 kg m −3 ).In A, a correction factor D is taken into account for scale distortion in the polar-stereographic projection grid cells away from the standard parallel.D is www.the-cryosphere.net/8/1725/2014/The Cryosphere, 8, 1725-1740, 2014  calculated by where φ and φ s are the local latitude and the latitude at the standard parallel (71 • N), respectively.As Eq. ( 5) indicates, total mass changes contains a component caused by changes in SMB ( Ṁsmb ) and a component caused by ice dynamics ( Ṁdyn ).The latter is calculated from altimetry by correcting the total observed Ḣalt for changes in firn compaction Ḣfc and SMB ( Ḣsmb ).Ḣsmb is directly related to SMB through the surface density.Both Ḣsmb and Ḣfc are derived from RACMO output (Sect.4).Because of the combination of modelled and observed quantities, any inconsistencies between the two will show in the results.Because neither model nor observations are perfect and the error margins are relatively large, some inconsistencies are expected.Figure 7 shows Ḣalt after correcting for Ḣfc and Ḣsmb .Indeed, there are discrepancies between the modelled SMB and Ḣalt .In particular, underestimation of accumulation anomalies or overestimation of melt anomalies by the model lead to positive dH /dt values that are associated with ice dynamics in Eq. ( 5) and assigned ice density.Especially in the ice sheet interior, this is not realistic.
As a consequence, we treat dH /dt as a composite of dynamically induced dH /dt and Ḣsmb only in areas where dynamically induced dH /dt is expected, i.e. areas with ice velocity above a certain threshold, and where dynamically induced dH /dt was found to occur in the analysis described in Sect.3. Here, we choose a velocity threshold of 70 m yr −1 .In other areas, we assume all dH /dt's to be caused by SMB.The results are slightly sensitive to the value of this threshold (a range of 30 to 100 m yr −1 gives a range of 14 Gt yr −1 ).Using 70 m yr −1 , the results are very similar to those obtained when simply the surface density (which is equal to ice density in the ablation zone) was used for all volume change, or when ice density was assumed for all areas below the equilibrium line altitude and surface density for the remaining areas.This is explained by the fact that the areas with relatively high velocity largely overlap with the ablation zone, so ice density is assigned to changes both due to SMB and ice dynamics.Using a relatively high threshold avoids ice density being erroneously assigned to SMB anomalies.
To quantify the error, we assume the total uncertainty in dM/dt (σ M ) is an uncorrelated combination of the error in the interpolated volume changes and errors in the firn model and its forcing (σ smb , σ H smb , σ H fc , and σ ρ s ; see Section 4).The error in the volume changes σ H + alt is a (again uncorrelated) combination of the error in the altimetry data, taken as the standard error on dH /dt from the regression, and the interpolation error.We calculate the interpolation error as the square root of the Kriging variance, minus that of the employed nugget (which represents point-scale variance).Formally, the Kriging variance is not an error estimate.However, as it is zero at locations containing data and increases with distance, it provides a measure of interpolation uncertainty.
So, the two components of mass change ( Ṁsmb and Ṁdyn ) and their errors (σ M smb and σ M dyn ) are calculated as follows.
If velocity ≥ 70 m yr −1 and the regression slope (Table 1) and otherwise Now we have a mass change rate and error at the grid cell scale, which then need to be aggregated to the ice sheet (or basin) scale.Ṁ for all grid cells can simply be added together, but aggregation of σ M is more complicated.Taking the RMS assumes all grid cells are independent, whereas a sum assumes they are all dependent.In reality, only grid cells within a certain decorrelation distance (D decor ) are dependent on each other.The decorrelation distance close to the ice sheet margin is expected to be different from that in the interior.Therefore, we follow Rignot et al. (2008a) and calculate D decor separately for all areas above and below an elevation of 2000 m, by taking a random subset (n = 5000) from all grid cells and fitting an exponential function to the decay curve of correlation with distance.In addition, to account for differences between the spatial structure of data based on radar and laser altimetry, D decor was estimated using all data from either 1995-2003 or 2003-2008 separately and then used for the appropriate years.Based on the two D decor estimates, the number of regions that can be considered to be independent N ind can be calculated (van de Berg, 2008): where A b is the area of interest, i.e. the area above or below 2000 m of a basin or ice sheet.The spatially aggregated error σ agg is then where N is the number of grid cells in the area of interest and A is the grid cell area.The two zones (above and below 2000 m) are assumed independent, and the resulting errors are added by taking the RMS.The resulting error is somewhat sensitive to the cut-off point of the correlation decay curve.The e-folding distance (the distance at which the correlation coefficient drops below 1/e, ≈ 0.368) is commonly used and produces a relatively low error estimate because of the large number of independent zones.On the other hand, at lower cut-off values the curve was found to approach an asymptote.As a trade-off between these two, a correlation coefficient of 0.25 was used as the cut-off to define D decor , resulting in a fairly conservative estimate and relatively large error bounds.

Results and discussion
According to the methods described in Sect.5, time series of mass change rate and their errors are calculated for 38 individual draining basins and the entire ice sheet.The temporal resolution is nominally annual, and the time series are relatively smooth because each annual dH /dt estimate includes data from adjacent years, and the spatiotemporal interpolation tends to smooth the time series.The resulting time series are shown in Fig. 8, and for 2 years (1998 and 2007) values for individual drainage basins (Rignot et al., 2008b) are shown as maps.In addition, mass change rates and errors for two 6-year time periods for the same drainage basins are shown in Table 2.
The impact of the different SRA processing on the estimated Ṁ time series can be seen by comparing the solid black and blue lines, which were derived from the Davis and Khvorostovsky processing approaches, respectively, using Table 2. dM/dt estimates for two 5-year periods (1996-2001 and 2003-2008) for all 38 drainage basins and the entire GrIS, including estimated errors.Units are Gt yr −1 .Basin  1996-2001  2003-2008 Basin  1996-2001  2003-  ST-KED.There are differences of a few tens of cubic kilometres for individual years due, we believe, to differences in sampling strategies and filtering criteria between the two data sets, but, importantly, the gradient of changes is consistent between the two methods.Differences are generally substantially less than the uncertainty bounds for the mass trends shown by the grey shading in Fig. 8 and are similar in magnitude to the differences between the various interpolation approaches.Dotted lines in Fig. 8 show the sensitivity of the resulting Ṁ to the interpolation method used.The difference between the interpolation methods is typically 10-20 %.The average difference between ST-OK and ST-KED over all individual years is 21 Gt yr −1 .This is mostly within the error bars of ST-KED, with the exception of the period between about 2001 and 2004.In this period, the sampling density is particularly low because it is the end of the ERS-2 period and the beginning of ICESat.OK and KED produce a positive mass change rate, based on very few data.ST-OK and ST-KED use data from adjacent years, resulting in a smoother time series.The same holds for approximately the end of the ICESat era in 2009.According to Fig. 8, the mass change rates were modest and mainly caused by SMB until about 2000/2001, when both SMB decreased and mass loss due to ice dynamics increased.From a study on Jakobshavn Isbrae (Hurkmans et al., 2012b), it was apparent that the onset of the strong thinning appeared delayed in the interpolated SRA results compared to airborne laser altimetry from the ATM.Because the thinning started close to the grounding zone, where SRA observations were largely absent, and propagated upstream, it took time for the thinning signal to extend far enough inland to be captured by SRA and thus appear in the mass change rates.This issue may occur at other major outlet glaciers as well, which is why the onset of the decrease in dM/dt in Fig. 8 appears delayed compared to IOM results (Rignot et al., 2011).Spatiotemporal interpolation only partly ameliorates this issue (Fig. 8).Mass loss increases after 2000 until about 2006, and then decreases slightly between 2006 and 2009.This is broadly consistent with gravity-derived mass trends which show a reduction in mass loss in 2007 (Rignot et al., 2011).Other studies (e.g.Khan et al., 2014) indicate that after 2010 ice loss has increased further.
As mentioned earlier, SRA provides one of the longest continuous time series of mass trends for the ice sheets, but there are relatively few comprehensive estimates from this approach published to date.It is instructive, therefore, to compare our results with a range of other published estimates for the GrIS.Table 3 provides a summary of other estimates in the literature, in combination with our results for the corresponding periods.In addition, Howat et al. (2011) estimated temporal mass variations from three major outlet glaciers (Helheim, Kangerdlugssuaq, and Jakobshavn Isbrae) at monthly resolution using the IOM. Figure 9 shows our mass variations for the same basins, together with resampled results from Howat et al. (2011), to match our temporal resolution.The temporal variability agrees reasonably well (see Fig. 2 in Howat et al. (2011)).However, for Kangerdlugssuaq and Jakobshavn, our mass change rates are less negative than those from Howat et al. (2011), whereas for Helheim our mass loss is more negative.For Jakobshavn, and probably also Kangerdlugssuaq, one of the causes is the delay in the onset of thinning that was discussed above and in Hurkmans et al. (2012b).Furthermore, as a third estimate for Jakobshavn Khan et al. (2010) gave a mass loss of about 20 Gt yr −1 between 2006 and 2009, which is comparable to our results.For Helheim, the variations in mass are more strongly influenced by SMB than discharge compared to the other two glaciers discussed here.Our results match the change in discharge better than the change in total mass (SMB discharge), which suggests that the regional climate model used in the IOM (RACMO) may be overestimating the SMB changes for this glacier.It should also be borne in mind that we are not comparing precisely the same quantities.In Fig. 9, the curves for the IOM method have been corrected for frontal retreat, which is not captured directly by altimetry.This can be a significant factor in the mass loss in some areas (see Fig. 3 of  2011) (red), resampled to our temporal resolution using a moving 36-month window, and SMB anomalies from RACMO2 (green).Howat et al., 2011) and is a potential weakness of altimeterderived mass trends.
Other estimates of GrIS mass change rates that are based on satellite altimetry have mostly used averages over multiple years.For instance, Zwally et al. (2005) used the period 1992-2002 (based on ERS-1 and ERS-2), and Sørensen et al. (2011) used the ICESat ERA (2003-2008).It is, therefore, not always possible to compare our results directly with other results in the literature, especially when those results are obtained from altimetry.In Table 3, we compare dM/dt estimates reported in the literature with our own estimate, where we calculated the average over the corresponding period.When the periods were not identical, we used various solutions to match the periods of interest (see footnotes of Table 3).For the period 1992-1994 we assume the ice sheet is close to balance, which seems reasonable considering the period 1995-1997 in Fig. 8, and we extend our time series for these years by 0 ± 50 Gt yr −1 .To extend the time series to the period after ICESat's lifetime, we used an extension and update of the GRACE results reported by We note, however, that GRACE solutions include some component of marginal glacier and ice cap mass loss (the precise proportion being dependent on the approach used for masking land-ocean mass exchange) and this may be as much as about 20 % of the ice sheet mass loss.
Our results agree rather well with most of the other results.Some of the exceptions are other altimetry-based results, and most of these differences can be explained.As mentioned before, Zwally et al. (2005) augmented radar altimetry data with airborne data (ATM) over outlet glaciers and also pre-sented results without using ATM.Both are shown in Table 3, and -while the estimate without ATM is, as expected, much higher (61 ± 3 versus 6 ± 38 Gt yr −1 ) -the updated estimate from Zwally et al. (2011) (7 ± 3 Gt yr −1 ) is close to ours for the period 1992-2002.Our interpolation algorithm, aided by velocity as a secondary variable, thus produces similar results without using additional (airborne) altimetry data.The estimates by Li and Davis (2008) do not utilise other data and are more positive compared to ours, as ice losses close to the margin are not adequately captured.Our results are similar to the largest loss value determined by Sørensen et al. (2011) from ICESat only, most likely due to the fact that we explicitly account for poorly sampled, high-elevationrate areas.The near doubling in mass loss that Khan et al. (2014) et al. (2008) produce lower mass losses compared to our results and to newer estimates that include more recent years (Sasgen et al., 2012;Chen et al., 2011).The ICE-Sat mass change rate estimate from Zwally et al. (2011) is also more positive than ours and another nearly contemporaneous ICESat-based estimate.Our result, however, does agree very well with Sasgen et al. (2012), where the ICESat estimate used is identical to that in Sørensen et al. (2011).A recent mass change rate estimate that combines various results discussed here (Shepherd et al., 2012) obtains a mass change rate over 1992-2011 for the GrIS of −142 ± 49 Gt yr −1 , whereas our estimate, adjusted to match the period, is −114 ± 42 Gt yr −1 .For the period 1992-2002, Shepherd et al. (2012) only use IOM results, in which the onset of mass loss occurs earlier than in our results (Fig. 8), possibly explaining the difference.The averaged estimate by Shepherd et al. (2012Shepherd et al. ( ) for 2003Shepherd et al. ( -2008, on the other hand, is nearly identical to ours.

Conclusions
Using ERS-2 data for the period 1995-2003 and ICESat data for 2003-2009, we have reconstructed the time-evolving mass change rates of the GrIS for the combined period.Due to our regression approach and the use of the spatiotemporal interpolation methodology, ST-KED, the effective temporal resolution is about 3 years.ST-KED also uses velocity as auxiliary information to constrain the spatial pattern of elevation trends over outlet glaciers where altimetry measurements are sparse.The underlying assumption is that the spatial patterns of velocity and elevation change rate are similar up to linear rescaling (Deutsch and Journel, 1992).This was shown to be the case for Jakobshavn Isbrae by Hurkmans et al. (2012a) and is now verified for other major outlet glaciers on the GrIS.It was found that the regression slope between velocity and elevation change rate was a useful metric for whether an outlet glacier is dynamically thinning, and it improved dH /dt estimates (with respect to airborne validation data) on glaciers where this is the case.ST-KED was then only used on glaciers where dynamic dH /dt dominates the SMB-related changes, and ST-OK was used for the remaining area.Elevation changes were corrected for firn compaction using RACMO2 data in combination with a firn model (Reeh, 2008), and the appropriate densities from the model were used for the volume-to-mass conversion, assuming that elevation changes in fast-flowing areas are caused by both ice dynamics and SMB, and by SMB only elsewhere.Considerable care was taken in assessing the error sources and how they are spatially correlated when determining the aggregated mass trend error.
Until about 2000, dM/dt is largely caused by SMB and agrees well with interannual variability of SMB as modelled by RACMO2.After that, both changes in SMB and ice dynamics contributed roughly equally to the increasingly negative mass balance.This is consistent with several earlier estimates using IOM and GRACE (van den Broeke et al., 2009).A comparison with IOM results that were truncated to our period of interest (from Rignot et al., 2011) and resampled to the same temporal resolution (Fig. 8) indicates that, according to the IOM, the onset of the mass loss is around 1998, followed by a gradual increase to about 200 Gt yr −1 in 2004.Our results produce a more sudden acceleration in mass loss in 2002.From 2004 onwards both methods are very similar.The apparent delay in the onset of mass loss was also observed at Jakobshavn Isbrae (Hurkmans et al., 2012b) and is caused by SRA sampling issues: it takes time for the dynamic thinning to propagate far enough upstream to be captured by radar altimetry.
A comparison of our results with other published estimates for the GrIS over various periods generally shows that satellite altimetry is in reasonable agreement with other methods and other laser-altimetry-based estimates.Mass change rate estimates based on SRA alone are generally more positive than our results (for reasons outlined earlier), but for the combination of SRA and ATM data results are consistent (Zwally et al., 2011).Some issues over outlet glaciers remain, for example the inability of altimetry to account for mass loss from grounding line retreat and the poor sampling close to the grounding zone that was discussed in the previous paragraph.
The good agreement between our results and other estimates, especially the IOM, which was the only estimate to be used in the averaged estimate by Shepherd et al. (2012) prior to 2003, provides confidence to extend our time series back to 1992 and forward in time with CryoSat-2.To overcome problems caused by the sparse crossover density in our results, it would be worthwhile to redo the analysis with a repeat-track approach for the SRA data, such as the approach that was recently carried out for Antarctica (Flament and Rémy, 2012).The disadvantage of this is that it requires exact repeats, precluding much of the ERS-1 mission.
Author contribution.R. T. W. L. Hurkmans carried out the work and wrote the paper, J. L. Bamber conceived the study and contributed to the writing of the paper, C. H. Davis and K. S. Khvorostovsky provided their SRA time series, B. S. Smith contributed to the analysis of the ICESat data, I. Joughin provided the velocity data, and N. Schoen interpolated the K. S. Khvorostovsky data using the techniques discussed.All authors reviewed and commented on the manuscript.

Figure 1 .
Figure 1.Velocity mosaic (a), and elevation change rates as derived from ERS-2 (b) and ICESat (c).The velocity mosaic is plotted logarithmically for clarity and represents the average for the years 2000 and 2005-2008.The data used to calculate the elevation change rates in (b) and (c) are from 1995-2003 and 2003-2009, respectively.

Figure 2 .
Figure 2. (a): dH /dt versus velocity for selected basins: Petermann glacier (P), Jakobshavn Isbrae (J), Storstrømmen (S), Kangerdlugssuaq (K), Helheim glacier (H), Upernavik (U), and Nioghalvfjerdsbrae (N), using all ICESat data (2003-2009).The slope of the linear regression is shown as well with the corresponding r in brackets.The inset shows the scatter plot for all of Greenland.Note that the x axis is limited at 3 km year −1 for clarity.The few points that have higher velocities are included in the relationship but not in the plot.(b) and (c) show the slope and Pearson's correlation coefficient (multiplied by 10), respectively, for all 38 Greenland drainage basins.

Figure 3 .
Figure 3. Interpolation results for selected test basins and averaged ICESat data (2003-2009).The map shows interpolated elevation change rates using ST-KED for the entire GrIS and the outline of six selected glaciers.The remaining plots show elevation change rate and velocity along a transect following the centreline of each glacier.In addition, ICESat points (brown) and ATM measurements (green stars) along the transect are shown.Note that the ATM data were not used in the regression.Surface velocity along the flowline is shown by the solid black line.

Figure 6 .
Figure 6.From left to right: σ smb , σ H smb , σ H fc , and σ ρ as obtained by propagating errors in SMB and temperature from RACMO2 (Eq.4) through the firn model.

Figure 8 .
Figure 8.(a) Time series of Greenland dM/dt from ST-KED using the Davis SRA data (solid black line) and Khvorostovsky SRA data (solid blue line) combined with ICESat.For comparison, the SMB anomaly from RACMO2 (resampled to the same temporal resolution as the altimetry data) is shown in green, and in dashed lines also total dM/dt resulting from the other three Kriging methods (OK/KED/STOK) are shown.The red line shows the time series produced by Rignot et al. (2011) using the IOM, smoothed over 3 years.(b) and (c) show temporal snapshots of the dM/dt per basin in 1998 and 2007.

Figure 9 .
Figure 9. Ṁ for three major outlet glaciers -Kangerdlugssuaq, Helheim, and Jakobshavn Isbrae-over the period of interest (blue).For comparison, also shown are the mass balance estimates from the IOM from Howat et al. (2011) (red), resampled to our temporal resolution using a moving 36-month window, and SMB anomalies from RACMO2 (green).

Table 1 .
Statistics of a linear regression between velocity and dH /dt, for 38 individual drainage basins (from

Table 3 .
Overview of mass change rate estimates in the literature compared to estimates resulting from this study.Estimates that are outside our estimated errors are printed in bold font.Periods were adjusted such that for every row both estimates consider the same period (see footnotes for details).Units are Gt yr −1 .

www.the-cryosphere.net/8/1725/2014/ The Cryosphere, 8, 1725-1740, 2014 to
find based on ICESat, i.e. 172 Gt yr −1 for 2003-2006 292 Gt yr −1 for 2006-2009, is a larger increase and larger absolute value than we obtain for the latter period.We note, however, that their GRACE-based estimate of 257 Gt yr −1 for 2006-2009 is identical to our value for this period.For the GRACE results,Luthcke et al. (2006)and Wouters