Corresponding author: R. T. W. L. Hurkmans, Bristol Glaciology Centre, School of Geographical Sciences, University of Bristol, Bristol BS8 1SS, UK. (firstname.lastname@example.org)
 Estimation of ice sheet mass balance from satellite altimetry requires interpolation of point-scale elevation change (dH/dt) data over the area of interest. The largest dH/dtvalues occur over narrow, fast-flowing outlet glaciers, where data coverage of current satellite altimetry is poorest. In those areas, straightforward interpolation of data is unlikely to reflect the true patterns ofdH/dt. Here, four interpolation methods are compared and evaluated over Jakobshavn Isbræ, an outlet glacier for which widespread airborne validation data are available from NASA's Airborne Topographic Mapper (ATM). The four methods are ordinary kriging (OK), kriging with external drift (KED), where the spatial pattern of surface velocity is used as a proxy for that of dH/dt, and their spatiotemporal equivalents (ST-OK and ST-KED). KED assumes a linear relationship between spatial gradients of velocity anddH/dt, which is confirmed for both negative (Pearson's correlation ρ < −0.85) and, to a lesser degree, positive (ρ = 0.73) dH/dtvalues. When compared to ATM data, KED and ST-KED yield more realistic spatial patterns and higher thinning rates (over 20 m yr−1 as opposed to 7 m yr−1for OK). Spatiotemporal kriging smooths inter-annual variability and improves interpolation in periods with sparse data coverage and we conclude, therefore, that ST-KED produces the best results. Using this method increases volume loss estimates from Jakobshavn Isbræ by up to 20% compared to those obtained by OK. The proposed interpolation method will improve ice sheet mass balance reconstructions from existing and past satellite altimeter data sets, with generally poor sampling of outlet glaciers.
If you can't find a tool you're looking for, please click the link at the top of the page to "Go to old article view". Alternatively, view our Knowledge Base articles for additional help. Your feedback is important to us, so please let us know if you have comments or ideas for improvement.
 Three main approaches exist to estimate the mass balance of the Greenland and Antarctic ice sheets, each with their own errors and uncertainties. First, mass balance may be estimated from the difference between the net surface mass balance and ice flux through a gate [Rignot and Kanagaratnam, 2006; Rignot et al., 2008b]. Secondly, variations in the gravity field can be used to estimate mass changes of the ice sheets using the GRACE satellites (Gravity Recovery and Climate Experiment) [Wouters et al., 2008; Luthcke et al., 2006]. Thirdly, mass changes can be estimated from elevation change (dH/dt) measurements using radar or laser altimetry [e.g., Zwally et al., 2005, 2011; Sørensen et al., 2011].
 To obtain an estimate of volume change, point-scaledH/dt measurements need to be integrated over the basin and thus require spatial interpolation. A method that is usually employed to achieve this is ordinary kriging (OK) [e.g., Zwally et al., 2005; Bamber et al., 2001]. Typically, the largest dH/dt values occur in the relatively narrow parts of outlet glaciers with high velocities [Krabill et al., 2004; Joughin et al., 2008a; Pritchard et al., 2009]. With conventional satellite radar altimetry, however, it is difficult to obtain sufficient measurements in these areas: because the radar beam width is about 20 km [Thomas et al., 2008], narrow outlet glaciers (such as Jakobshavn Isbræ) cannot be adequately sampled, and steep slopes cause a displacement of the signal, as the return signal does not originate from nadir but from the highest point within the radar footprint [e.g., Brenner et al., 2007]. While this slope-induced error can be corrected for to some degree, small-scale relief can result in noisy height estimates and loss of lock [Bamber and Gomez-Dans, 2005]. Another issue is the spatial coverage; the horizontal resolution is bounded by either crossover points (exact locations where ascending and descending satellite tracks cross each other) or across-track spacing [Thomas et al., 2008]. These problems degrade performance in the areas where the largest changes occur: the relatively narrow outlet glaciers along the ice sheet margin. For future studies, the recent launch of CryoSat-2 may partially alleviate these problems due to its higher resolution and latitudinal coverage (up to 88°N). Satellite laser altimetry can give an increased accuracy over outlet glaciers, but the laser instrument on the ICESat satellite (GLAS; the Geoscience Laser Altimeter System) suffered from operational problems [Abshire et al., 2005], resulting in reduced spatial and temporal resolution. Stereo photogrammetry can be used to determine volume change locally [e.g., Stearns and Hamilton, 2007], but its accuracy and applicability limits it to marginal areas with adequate contrast and large elevation differences. It is not suitable for ice sheet wide volume change estimation.
 Thus, given the limited number of valid dH/dtmeasurements in the high-velocity parts of the ice sheet (greater than 100 m yr−1), it is unlikely that an interpolation obtained by, for example, OK will reflect the true spatial pattern of dH/dt. The similarity of the spatial patterns of surface ice velocity and dH/dt [e.g., Rignot et al., 2008a; Joughin et al., 2008b; Pritchard et al., 2009] suggests that velocity, which is available from radar interferometry for nearly all of the major outlet glaciers [Joughin et al., 2008a; Howat et al., 2007], contains information that can aid the interpolation. In this paper, we employ an alternative interpolation method, kriging with external drift (KED), that takes into account auxiliary data [Deutsch and Journel, 1992]. KED has been applied in geo-sciences in several studies: for example bySnepvangers et al.  in interpolation of soil moisture values, where net precipitation was used as the external drift, or by Schuurmans et al.  for combining precipitation measurements from radar with rain gauges. An additional extension to the interpolation method is the use of spatiotemporal kriging, which handles irregular sampling by interpolating both in space and time [Snepvangers et al., 2003; Kyriakidis and Journel, 1999; De Cesare et al., 2001]. The resulting interpolation technique, spatiotemporal kriging with external drift (ST-KED), has two advantages compared to OK: first, it makes use of more measurements because it samples the past and the future with respect to the period of interest; and secondly, velocity gradients constrain the interpolation on fast-flowing regions where there are no, or limited, data.
 As a case-study area to test whether ST-KED indeed improves the estimation of volume change, we use Jakobshavn Isbræ, Greenland's most active outlet glacier. Jakobshavn Isbræ has been the subject of many investigations regarding its recent acceleration and thinning [Joughin et al., 2004; Thomas et al., 2003; Motyka et al., 2010, 2011; Luckman and Murray, 2005], and has been densely surveyed by airborne laser altimetry [Abdalati and Krabill, 1999; Krabill et al., 2000, 2004]. The measurements by the ATM, mostly located in the high-velocity region, are ideal for independently testing our approach. Using these validation data, we compare the interpolated patterns ofdH/dt, thinning rates and errors.
2. Study Area and Data
2.1. Study Area: Jakobshavn Isbræ
 Jakobshavn Isbræ is Greenland's largest outlet glacier, draining about 5.4% of the ice sheet (92,080 km2) [Motyka et al., 2011]. Its location and size are shown in the inset of Figure 1. The glacier converges into a 4 km wide, deep (over 1000 m below sea level) trunk, flowing into a deep fjord where, until recently, it formed a floating ice tongue of about 15 km length. Between 1850 and 1950, the position of the calving front of this tongue retreated some 25 km along the fjord. After 1950, the calving front was more or less stable until about 1998 [Thomas et al., 2009; Csatho et al., 2008]. After 1998, rapid melting of the floating ice tongue caused a retreat of the calving front, until it finally disintegrated in 2003 [Joughin et al., 2004]. This collapse was attributed to enhanced inflow of warm water into the fjord and subsequent increased basal melting of the tongue [Holland et al., 2008; Motyka et al., 2011].
 Along with the retreating calving front, measurements of surface velocity [Joughin et al., 2004, 2008a] and dH/dt [Krabill et al., 2004] indicated acceleration and thinning of the high-velocity part of the glacier. While this area had been thickening prior to 1997 [Abdalati et al., 2001; Thomas et al., 2003], it has been thinning since then, with rates in excess of 15 m yr−1 near the grounding line. Over the same period, Joughin et al. measured ice velocity, averaged over some 20 locations in the high-speed area of Jakobshavn Isbræ. Between 1985 and 1992, velocities were slowing down from 6.7 km yr−1 to 5.7 km yr−1. After 1999, however, the glacier has been speeding up, to 9.4 km yr−1 in 2000 and 12.6 km yr−1 in 2003 [Joughin et al., 2004].
 As input for our interpolation, we employ data from two spaceborne altimeters: i) the RA-1 (Radar Altimeter) on ERS-2 (European Remote Sensing satellite); and ii) the laser altimeter (GLAS) on ICESat. In addition, we use data from an airborne laser altimeter (NASA's ATM) for validation purposes. ESA's ERS-2 was launched in 1995 and has been in an orbit at an altitude of 780 kilometer with a 35-day repeat cycle until September 2011 when it was decommissioned. In this study, ERS-2 data between 1995 and 2003 have been used.dH/dtvalues derived from ERS-2 are based on crossover points that were averaged into clusters [e.g.,Li and Davis, 2008], shown as blue triangles in Figure 1. ICESat was launched in 2003 and carried a laser altimeter system (GLAS). Due to rapid laser degradation, the system was switched on for only two to three campaigns of 33 days per year [Abshire et al., 2005]. The spacing between tracks is therefore relatively large, about 20 km for Jakobshavn Isbræ (Figure 1). Along-track spacing is approximately 172 m, with a footprint size of about 60 meters [Zwally et al., 2002]. We use ICESat data (release 31) from 14 epochs: from February/March 2003 until February/March 2008 [Zwally et al., 2010]. The same preprocessing was carried out as in Sørensen et al. . ICESat repeat tracks are shown in red in Figure 1. ATM is also a laser altimeter, but airborne. Since 1993, repeated flights have been performed over Greenland [Krabill et al., 2000, 2004], with sampling over parts of Jakobshavn Isbræ [Abdalati and Krabill, 1999] shown in Figure 1. ATM measures ice-surface elevation with an accuracy of 10 cm or better [Krabill et al., 2000], within a swath of 150–200 m at a typical aircraft altitude of 500 m [Thomas et al., 2009]. Here, a processed version of the data was used, in which footprints (with size 1–3 m) were resampled to 70 m overlapping “platelets”, one on each side of the aircraft [Thomas et al., 2009].
 Ice velocity is derived from a combination of interferometric synthetic aperture radar (InSAR) and speckle tracking [Joughin, 2002]. In the winter of 2000/2001, Canada's RADARSAT-1 acquired near-complete coverage of the Greenland ice sheet. From 2005/2006 this coverage was repeated almost annually [Joughin et al., 2010]. We use a mosaic of RADARSAT-1 derived velocities averaged over the years 2000, and 2005–2008. To define the outline of the drainage basin of Jakobshavn Isbræ, finally, the delineation ofRignot and Kanagaratnam  is followed.
3.1. Obtaining dH/dt Values
 To be able to combine all data sets described before, we define a grid in polar-stereographic projection (central latitude 71°N; central longitude 39°W), at a 1 km resolution. The processing approach for laser altimetry (ATM, ICESat) differs from that for radar altimetry. The radar altimetry data were already preprocessed as crossover clusters, but not corrected for the slope-induced error, i.e., the fact that the return signal does not come from the point at nadir, but from the point within the rather large footprint that is closest to the antenna [e.g.,Brenner et al., 1983]. Near the ice sheet margin, the difference between these points can be considerable (about 14 km for a 1 degree slope). This is corrected for using slope and aspect from a digital elevation model (DEM) [Bamber et al., 2001] as described in Hurkmans et al. .
 For ICESat, several approaches have been used to extract dH/dt [Howat et al., 2008; Pritchard et al., 2009; Slobbe et al., 2008]. Because ICESat tracks almost never exactly repeat, a regression approach can be used in which slope (both across-track and along-track) anddH/dt are simultaneously solved for. To estimate dH/dt from ATM, footprints are resampled into 70 m platelets as mentioned in Section 2.2, and elevations at platelet centers of overlapping flightlines are compared [Krabill et al., 2004; Thomas et al., 2009]. However, in all these studies only one average dH/dt was obtained, where temporal resolution was sacrificed in a way to increase spatial coverage [Pritchard et al., 2009]. Because we are interested in the change of thinning or thickening rates through time, we need a higher temporal resolution as well.
 First, for each repeat track (or flightline for ATM), we average all footprint elevations within a given 1 km grid cell. We again use slope and aspect from a DEM [Bamber et al., 2001] to project all elevation footprints (for ICESat) or platelet centers (for ATM) toward the center of the pixel. This approach is likely to be less accurate than a method that estimates slope directly from the data as was done in, for example, Pritchard et al. . However, we expect this to add noise to the results rather than a bias. Because we employ rather strict data culling criteria (see below), we do not believe our results are significantly influenced by this. Time-averageddH/dt values are also calculated for each satellite (using regression including all available data). For ICESat, the spatial pattern of overall dH/dt for the entire Greenland ice sheet (2003–2007) compares very well with those reported by Pritchard et al.  and Sørensen et al.  (not shown).
 For both radar and laser altimetry, we then calculate annual values for dH/dtat each grid cell for every year between 1997 and 2008 using linear regression. For each year, data are used from a three-year period centered on the year under consideration, because the temporal sampling, especially before 2003, is too poor to allow reliabledH/dt calculations within single years. Regression is only carried out in grid cells where i) data from at least five overpasses are available that ii) span at least one year and iii) have a standard error on the resulting dH/dt rate of less than 0.40 m yr−1. The standard error on the regression coefficient SEcoef (in this case dH/dt rate) is calculated by:
where e is the vector of residuals, n is the sample size, and x is the input with mean . It should be noted that this standard error is not equivalent to measurement error, but takes into account sample size, and variance of both input and residuals of the regression as well. The threshold is selected by trial and error to avoid a noisy spatial pattern of points that are close together and opposite in sign, usually because the regression was based on a small subset of overpasses. In addition, prior to regression, measurements that are outside the range of average plus or minus two standard deviations (based on the entire available time series) are removed. The number of data points that was removed in the culling process is variable between data sets and time periods. For ICESat, typically about 25% of data points are removed, mostly because less than five overpasses were available or the spanned time range was too short. Figure 2 shows estimated dH/dt for different periods, based on linear regression. Figure 3 shows the associated standard errors of the dH/dt.
3.2. Kriging Theory
 A widely used interpolation method in geosciences is kriging, sometimes referred to as “optimal interpolation”. It is a form of linear regression minimizing an estimation variance [Deutsch and Journel, 1992], of which several varieties exist. In the following, the basic theory and the differences between the main types are briefly explained. For more information, see for example Deutsch and Journel , Herzfeld et al. , and Hengl et al. .
 In kriging, a value Z(x0) at an unsampled location can be estimated from N samples located at xi using:
where Z is treated as a random field with a trend component, m, and a residual component R = Z − m, and λi is the weight assigned to each sample location. R is treated as a random function with a mean of zero and a covariance C(h). C(h) can be estimated from the data in the form of a semi-variogramγ(h):
where N(h) is the number of data pairs at distance h, and Z(xi) and Z(xj) are the data values at points xi and xj, which are separated by a distance h.The semi-variogram thus indicates the expected squared increment of the data value over distanceh, and is characterized by three parameters [e.g., Herzfeld et al., 2000]: the nugget (point-scale variance), the range (largest distance at which there is spatial dependence), and the sill (the semi-variance at distances beyond the range). The covarianceC(h) can then be obtained by C(h) = C(0) − γ(h) where C(0) is the sill. The sill and range parameters are derived by fitting a model to the semi-variogram calculated from the data. Typically, an exponential, gaussian, or spherical equation is used for this [Deutsch and Journel, 1992].
 In most geostatistical applications, data to be interpolated do not only vary in space, but also in time. Therefore, several models have been developed for space-time variability [Kyriakidis and Journel, 1999]. One of these is the product-sum method [De Cesare et al., 2001], which is more flexible compared to other models, does not require an arbitrary space-time metric, and it is relatively straightforward to fit models to the sample semi-variogram [Gething, 2006]. In case of spatiotemporal kriging, a two-dimensional sample semi-variogram is obtained using a distance laghs and a time lag ht as dimensions, so that for every hta spatial semi-variogram is obtained and vice versa. A semi-variogram model is then fitted to both the spatial (average of allht's) and temporal (average of all hs's) sample semi-variogram. Using the product-sum model, the modeled spatial and temporal semi-variograms are then combined:
where γs(hs) and γt(ht) are the spatial and temporal semi-variogram, respectively;Cs(0) and Ct(0) are the spatial and temporal sills, respectively; and k1, k2, and k3 are given by:
where Cst(0) is the spatiotemporal sill, which is normally taken as the maximum value of the sample space-time semi-variogram surface [De Cesare et al., 2002]. The space-time covariance is then used for three-dimensional interpolation (x, y, and t), resulting in an interpolated surface at every pixel and at every desired time step, where also samples at earlier or later time steps are used in the estimate.
 The treatment of the trend m is the main difference between the various members of the kriging family. In this study, ordinary kriging (OK) and kriging with external drift (KED) are employed. The system of equations that needs to be solved for interpolation is derived by expressing the variance of the estimation error in terms of the covariances C(h), and then minimizing it by equating the partial derivatives with respect to λis to zero (see Deutsch and Journel  for details). In the simplest case, where m is assumed uniform, this yields:
where C(xi − xj) is the covariance function for lag h = xi − xj and C(xi − x0) is the covariance at lag (xi − x0). This case is called simple kriging (SK). When other types of kriging are employed the principle remains, but the form of the equation system in equation (6) changes slightly. In OK, the most widely used variant, m is assumed to be constant only for the local neighborhood of each estimation point (m(xi) = m(x0)). To be able to estimate m locally, equation (6) is extended with a constraint requiring the sum of weights to be one:
where μ(x0) is a Lagrange parameter associated with the constraint on the weights. Alternatively, a spatial trend can be assumed to exist in m.In kriging with a trend (KT), sometimes referred to as “universal kriging”, the trend is a linear or higher order function of the (x,y)-coordinates. In case of kriging with external drift (KED), a relation is assumed betweenZ and a secondary variable Y, so that the trend can be calculated from the secondary variable. Yis assumed to reflect the spatial variability of the variable to be interpolated and to be known at every location in the interpolation grid. In the case of KED, the trend is assumed to be linear (only a first-order trend):
 The system of equations to be solved is further extended, allowing local and implicit estimation of a0 and a1:
In this study, the FORTRAN programs GAMV and KT3D from the GSLIB geostatistical package [Deutsch and Journel, 1992] were used for calculation of the sample semi-variogram (GAMV) and the actual kriging (both OK and KED). An extension to GSLIB to allow for spatiotemporal kriging was provided byDe Cesare et al. , which we also employed here. It should be noted that all kriging types (OK, KED, KT, etc.) can be used in both spatial-only and spatiotemporal contexts. For example,Snepvangers et al. used space-time ordinary kriging (ST-OK) and space-time kriging with external drift (ST-KED) to model soil moisture contents. In the remainder of this study, we employ and compare four kriging types: OK, KED, ST-OK, and ST-KED.
3.3. Interpolation of dH/dt Values
 Many studies involving spatial interpolation of satellite measurements have used ordinary kriging (OK) [Zwally et al., 2005; Bamber et al., 2001; Herzfeld et al., 2000]. Considering the data coverage in Figures 1 and 2 however, it is unlikely that ordinary kriging, based on satellite altimetry alone, will yield a representative spatial pattern. To illustrate this, the maximum thinning rate obtained from ATM data (over the period 1999–2008) was about 17 m yr−1, at the location shown by the star in Figure 2. Note that this point is not located in the center of the main trunk, and some kilometers away from the grounding zone, so even higher values are possible. The maximum thinning rate as derived from ICESat, on the other hand, was only about 7 m yr−1 and values obtained by radar altimetry were even lower than that. A possible alternative interpolation method would be the hypsometric approach, in which a relationship is assumed between elevation and dH/dt [Arendt et al., 2002]. However, since dynamically induced dH/dt and elevation are not clearly related, the hypsometric approach is not further discussed here.
 Ice velocity measurements are available at nearly all locations on the ice sheet, especially over fast-flowing outlet glaciers. From visual inspection ofFigure 2, the spatial patterns of dH/dt and ice velocity seem very similar. Here, we investigate whether there is information in the ice velocity pattern that can aid the interpolation of dH/dt. As explained earlier, kriging with external drift (KED) uses auxiliary information to determine the underlying trend model for the interpolation. The secondary variable (velocity) is assumed to reflect the spatial trend in dH/dt up to a linear rescaling of units [Deutsch and Journel, 1992]. Spatial gradients in velocity should thus be related to spatial gradients in dH/dt.
 A spatial correlation between (dynamically induced) dH/dt and velocity has been noted before [Joughin et al., 2008b; Rignot et al., 2008a; Pritchard et al., 2009]. In fast-flowing areas, that are flowing faster than required to maintain mass balance, ice stretches longitudinally, which causes it to thin vertically [Rignot et al., 2008a]. The thinning increases surface slopes, which further increases velocity [Joughin et al., 2008b]. On outlet glaciers in both Antarctica and Greenland, significant thinning (and associated speed-up) is confined to the fast-flowing areas [Pritchard et al., 2009]. This is also the case for Jakobshavn Isbræ [Joughin et al., 2008b]. Therefore it is reasonable to assume that there is a relation between the spatial gradients of velocity and dH/dt. Such a relation would also depend on glacier width and thickness. Because Jakobshavn Isbræ is flowing through a deep trough and its thickness is >1000 m, and changes in width are also relatively small, changes in velocity are dominant in this respect. In addition, it should be stressed that such a relation is not fixed: its parameters (a1 and a0 in equation (8)) are determined implicitly from the kriging system of equations (equation (9)) and will be variable in space and time. Because of this flexibility, different spatial gradients of dH/dt can result from the same spatial gradient in velocity, constrained by the altimetry measurements. This allows us to use an average, static velocity field in order to maximize spatial coverage.
dH/dt is not only caused by ice dynamics but also by changes in surface mass balance, which will cause noise in any relationship between velocity and dynamically induced dH/dt.In the high-velocity area of Jakobshavn Isbræ, however, thinning rates are >15 m yr−1, whereas maximum SMB values are about −4 m ice yr−1 [e.g., Thomas et al., 2003]. Dynamic changes are therefore much larger than SMB-induced changes and we can assume that ice dynamics are the dominant signal.
 To assess the correlation between velocity and dH/dt, the two are plotted against each other in Figure 4. The plots and the highly negative Pearson's correlation coefficients (<−0.85) suggest a near-linear relationship, from which we conclude that velocity can be a suitable secondary variable for KED.Figure 4 shows correlation coefficients for the period 1999–2008 (2003–2008 for ICESat), but they were also calculated for individual years (Table 1). It was, however, not always possible to obtain a robust correlation coefficient for each year due to paucity of data. After 1999, correlation coefficients are reasonably high for all sensors.
Table 1. Correlation Between Velocity and dH/dt From All Instrumentsa
Correlation coefficients between velocity and dH/dt (ρ) and the number of measurements (N), for ATM, ICESat, ERS-2, and all data combined. Values are shown for individual years, based on data from moving 3-year periods, and averaged over the period 1999–2008 (ALL). Data before 1999 are excluded to limit the influence of opposing signs due to the transition from thickening to thinning in 1998 [Csatho et al., 2008; Thomas et al., 2003].
 Interestingly, positive correlation coefficients show up in the early nineties, in the period that Jakobshavn Isbræ was thickening [Abdalati et al., 2001] suggesting a dynamic origin. Apparently, the areas with the highest velocity after 2000, are also the areas that thickened most before 1997. Figure 5shows the corresponding scatterplots for ATM measurements in 1994 and 1996. For both years, essentially only one repeated flight line was available, so the number of samples is limited, and the significance of the trend is low. The presence of some warm summers in the mid-nineties [Thomas et al., 2003] and a simultaneous slow-down of the glacier [Joughin et al., 2004] suggest that the thickening was dynamic, and these positive correlation coefficients seem to confirm that. Between 1997 and 2003 correlations are lower because the glacier was in a transition between thickening and thinning, and very few valid dH/dtestimates are present before 2003. Because a 3-year window was required to calculatedH/dt, the period at the end of the ERS-2 era and the beginning of the ICESat era has a limited number of measurements. After 2003, when many ICESat epochs can be used, there is a strong anti-correlation (ρ < −0.80).
 To assess the value of including velocities in the interpolation, both OK and KED were used to interpolate the dH/dt data obtained in Section 3.1. In addition, to overcome the irregular temporal sampling by the different sensors, ST-OK and ST-KED were also used. One approach to do this would be to interpolate elevations prior to calculatingdH/dt. Due to the very different spatial coverage of ERS-2 compared to ICESat, many areas have coverage for only part of the period, causing artifacts indH/dt around these transition periods and noisy results. We therefore chose to calculate dH/dt prior to interpolation as described in Section 3.1 and interpolate the resulting annual values. Besides, as there is no clear relation between elevation itself and ice velocity, incorporating the external drift is not straightforward or sensible if elevations are interpolated directly.
 As a first step of the interpolation process, semi-variogram models were estimated (seeSection 3.2) using the GAMV routine in the GSLIB package (for OK/KED) or the extension to GSLIB by De Cesare et al. , after which range and sill parameters were estimated using standard curve fitting routines. Figure 6shows all resulting sample semi-variograms and models. Corresponding parameters and variogram model types (gaussian, spherical or exponential; see for exampleHerzfeld et al.  for details) are shown in Table 2.
Table 2. Semi-variogram Parameters for the Semi-variograms Shown inFigure 6a
Parameters (Range and Sill) as well as type (“Exp.”, “Gau.” and “Sph.”, are, respectively, exponential, gaussian and spherical), are shown for the isotropic case (1–D) as well as two directions (east-west, E–W and north-south, N–S) in the anisotropic case.
 For all kriging methods, we use the sample semi-variogram calculated from the satellite data only (ICESat and ERS-2), as we intend to use ATM data for validation only. The resulting semi-variogram and model are shown inFigure 6a. For comparison, the semi-variogram resulting from ATM data only is shown inFigure 6b. Because most of the ATM data are concentrated in the high-velocity area of the glacier (roughly a 50 × 50 km area), semi-variance for ATM is only considered and shown for spatial lags up to about 50 km. Because ice flow is mainly east-west, and most data are organized in (almost) north-south oriented tracks,dH/dtis expected to be anisotropic. Therefore, specific variograms were calculated for east-west and north-south directions, shown as gray and orange lines inFigures 6a and 6b. However, for lags up to about 60 km (all of the ATM data), semi-variograms for specific directions are very similar. Only for larger spatial lags do the east-west and overall semi-variograms diverge, but not significantly. Because distances in east-west direction within the basin are much larger than north-south, the north-south variogram should only, like ATM, be considered at smaller lags. Based on these results, we use the overall semi-variograms in the remainder of this study. Moreover, when the anisotropic semi-variograms were used for kriging the resulting anisotropy in the interpolated pattern was greatly exaggerated (not shown).
Figures 6c and 6dshow the spatiotemporal semi-variograms. From the sample semi-variogram surface (black surface inFigure 6d), the marginal space and time semi-variograms are derived, and models are fitted to them (Figure 6c). These are combined using the product-sum model [De Cesare et al., 2001] and the red square, representing the modeled semi-variogram surface (equation (4)) is obtained. The marginal temporal variogram is not particularly structured and not a very good fit is obtained. It should be noted, though, that the temporal resolution was probably not sufficient to capture a complete covariance structure because of the regression in moving three-year periods.
 Based on the parameters from Table 2, annual dH/dtestimates for the entire Jakobshavn Isbræ area were derived, using all four kriging variants. In addition, longer-termdH/dtvalues were interpolated using all data after 1999, i.e., 1999–2003 for ERS-2 and 2003–2008 for ICESat. Data before 1999 were excluded here because the change from thickening to thinning took place in around 1998 and assuming a linear change rate would thus be invalid.
Figure 7 shows temporally averaged dH/dt along a transect following the approximate line of maximum velocity from the grounding zone toward the accumulation zone. dH/dtvalues are shown for all four kriging methods, as are input (ERS-2 and ICESat) and validation (ATM) measurements along the transect. It should be noted that OK and KED use the temporally averageddH/dtvalues (which are also shown) as input, whereas for ST-OK/ST-KED interpolation results were averaged. InFigure 7b, therefore, the transects for ST-OK/ST-KED do not pass exactly through the input data. In addition,Figure 8 shows the spatial patterns of dH/dt for the same period. Input data are shown as well (Figure 8a), and Figure 8bshows interpolated ATM data (using ordinary kriging), which can be considered as the “true” spatial pattern because of the high sampling density in the area of interest. To get an estimate from ST-OK and ST-KED for the same period, the averagedH/dt over the years 1999–2008 was calculated (Figures 8e and 8f).
 Both Figures 7 and 8 show that interpolated thinning rates are never higher than measured ones (about 7 m yr−1) if no velocity data are used. This is clearly unrealistic (Figure 8a). As was mentioned in Section 3.3, ATM data show much higher thinning rates. The highest measurement, close to the grounding line in 2004 (Figure 2e), is about 25 m yr−1, which is an indication of the maximum thinning rate near the grounding zone.
 The interpolated spatial patterns using KED and ST-KED are evidently more realistic than those using OK/ST-OK compared to the interpolated ATM data (Figure 8). Also the interpolated thinning rates along the transect closely follow those of the ATM validation data (Figure 7). It should be noted that for the ATM data shown in Figure 7 the threshold on the standard error was slightly relaxed to be able to show more ATM data points (0.6 m yr−1 whereas it is 0.4 m yr−1 in the rest of this study).
 The velocity along the transect, which is also shown in Figure 7, rapidly increases from about 4 km yr−1 to 12 km yr−1 in the 10 km nearest to the grounding zone. ATM measurements along the transect show that the interpolated dH/dtvalues for KED and ST-KED are accurate upstream of that point. Closer to the grounding zone, however, thinning rates as predicted by KED/ST-KED follow a similar gradient to velocity (as would be expected from the external drift algorithm), rapidly reaching improbable values (up to −80 m yr−1). Two explanations are possible for this: first, while we excluded the area west of the approximate grounding line in 2004, the grounding line has been retreating in recent years [Csatho et al., 2008; Thomas et al., 2003] and might have migrated several kilometers further inland by 2008. Part of the area with very high velocities may thus have become ungrounded. Second, the sudden increase in velocity gradient suggests a change in flow regime. In both cases, the relation between velocity and dH/dt breaks down and we thus decided to exclude the area with velocities >4000 m yr−1, and the area within 5 km east of the 2004 grounding line. The combined area that is removed is relatively small (Figure 8).
 Statistics of dH/dt averaged over the same period are shown in Table 3. As we are focusing here on the high-velocity area, statistics were calculated using pixels associated with velocities ≥100 m yr−1 only (Figure 8). However, for the entire basin, average dH/dt values are also negative (not shown). This is in line with the net mass loss reported by other studies [e.g., Rignot et al., 2008b; Sørensen et al., 2011].
Table 3. Statistics of Data and Interpolations of dH/dt Values Over the Period 1999–2008a
For space-time kriging, the value for this period was obtained by averagingdH/dt over individual years. Shown are the mean, standard deviations, minimum and maximum, as well as the 10% and 90% percentiles. Note that for input and ATM data, the numbers (between brackets) are very much smaller than the total number of pixels in the area of interest (14,348) and thus heavily influence the statistics.
ATM data (450)
 To quantify the errors of all four interpolation methods, a cross-validation approach was carried out [Deutsch and Journel, 1992]. Using a separate run for each interpolation algorithm, every available altimetry measurement is left out and estimated from its neighbors. The error is then the difference between the interpolated and the actually measured value. Statistics of the errors, for individual years and the period 1999–2008, are shown in the form of box-plots inFigure 9a. For spatiotemporal kriging, this exercise could only be carried out for individual years. For cross-validation only the input data (i.e., no ATM) were used. These were located almost exclusively outside the high-velocity area. The influence of the external drift should, therefore, be non-existent or very small.Figure 9ashows that this is indeed the case, except for some years prior to 2002 where only (very few) ERS-2 measurements are available and one observation close to the high-velocity area has a noticeable effect. When comparing spatiotemporal kriging versus spatial-only kriging the differences are again small, although in most cases the errors are slightly smaller for ST-OK/ST-KED, because more (past and future) samples are taken into account. Exceptions are 1998 and 2003. 2003 is also the only year with a bias (i.e., the average of the errors is not zero). Because neither ERS-2 nor very many ICESat measurements are available for 2003 this can also be attributed to the low sampling density.
 In Figure 9b, box plots are shown for a comparison of interpolated values and ATM measurements where those are available. Because the number of ATM measurements changes dramatically over the years, interpretation of these errors is difficult. To compare results for the period 1999–2007, dH/dtfrom KED/ST-KED were averaged over these years. When many input measurements are available (typically the middle of the ICESat era and in the case of the averaged data), using the external drift reduces both the bias (mean error) and the maximum errors. If only few input data are available, i.e., only ERS-2, either all methods produce nearly the same errors, or maximum errors are slightly larger when using the external drift. There are two cases where errors are larger for KED (2004) or ST-KED (2002). This period marks the end of the ERS-2 coverage and the beginning of ICESat coverage. Because 3-year running windows were used to select data for the regression, in this period very fewdH/dt estimates with sufficiently low standard error are available. Hence, the external drift “overshoots” the ATM measurements in some of those occasions. Furthermore, it should be noted that in both years the available ATM measurements were clustered in a small area (Figures 2c and 2e).
 To further assess the temporal evolution of dH/dt, after interpolation, Figure 10shows the difference in time evolution between the methods. Six selected years from the period 1997–2007 are shown. It should be noted that from 2005 onwards the spatial patterns are similar. As mentioned before, prior to 2003 only data from ERS-2 was available, making it difficult to obtain a realistic pattern given the sparsity of ERS-2 (Figure 1). For KED and ST-KED, the sign of (sometimes very few) data points determines the direction in which the trend is calculated. For example in 2003, where the few satellite measurements that are available are predominantly positive (probably caused by a positive anomaly in the surface mass balance), values in the high-velocity area are predicted to be higher than the remaining part of the basin. This is not necessarily unrealistic, given the correlation between dynamic thickening and velocity that was discussed before (Section 3.3 and Table 1).
 Using spatiotemporal kriging, where measurements from adjacent years are employed, smooths the temporal variations considerably. Figure 10 shows that the transition between thickening and thinning is around 2000, which is about two years later than reported elsewhere [e.g., Thomas et al., 2009; Csatho et al., 2008; Joughin et al., 2008b]. However, as mentioned before, only ERS-2 data were available in this period, none of which is located around the high-velocity area. The dynamic thinning started near the grounding zone and then propagated upstream. It thus took some time to advance far enough upstream to be captured by ERS-2. The ST-OK and ST-KED patterns around this transition do not look realistic, as there is a sharp boundary between positive and negative values, caused by measurements with opposing signs. A higher sampling density in that period (e.g., additional sensors, or repeat-track processing of ERS-2 instead of crossovers) would alleviate this.
 In our results, kriging with external drift produces generally higher thinning rates that are closer to the ATM data. In addition, spatiotemporal kriging smooths effects of interannual variability in surface mass balance and it also improves estimates in periods with data gaps because it can use data from adjacent epochs. We believe, therefore, that ST-KED produces the best results. There are, however, some areas and periods where performance is less than optimal.
 Close to the grounding zone, where velocity rapidly increases from ∼4 km yr−1 to about ∼12 km yr−1in a distance of 10 km, thinning rates as interpolated by KED/ST-KED increase with a similar gradient, resulting in physically unrealistic values. In this area, the relation between velocity anddH/dtbreaks down because of uncertainty in the exact location of the grounding line. In addition, the large change in velocity gradient suggests that a different physical process becomes important. Furthermore, in the transition period from thickening to thinning, samples with opposing signs cause the trend from the external drift to have different signs, resulting in unrealistic patterns. This transition period, which in the high-velocity area happened around 1998, is delayed for some years in our interpolated time series because it took about two years for the thinning to propagate far enough upstream to be captured by radar altimetry.
 Because KED and ST-KED yield higher thinning rates over the fast-flowing part of Jakobshavn Isbræ, the resulting mass change might be different as well. The calculation of actual mass changes requires modeling density profiles and firn compaction, which is outside of the scope of this study. It is, however, straightforward to calculate volume changes by integratingdH/dt over the basin. Table 4provides volume changes for individual years between 1996 and 2007, as well as longer-term averages. In these estimates, the area with velocity >4000 m yr−1 and within 5 km of the 2004 grounding line was not taken into account, as discussed previously. These areas were assigned a dH/dt of −20 m yr−1, which is a reasonable value given the minimum value in ATM of −24 m yr−1. However, because this area is small, its impact on the integrated volume change is also small. As an extreme example, not assigning a value (i.e., 0.0 m yr−1) yields volume losses that are about 2 km3 yr−1lower. In general, volume change estimates for Jakobshavn Isbræ are about 10–20% higher for ST-KED than for ST-OK, with an average volume loss over 2003–2007 of 20.5 km3 yr−1for ST-KED versus 17.5 km3 yr−1for ST-OK (Table 4). Thus, the interpolation technique used can have a significant impact on the estimated volume changes at this scale. However, Jakobshavn Isbræ is one of the largest glaciers in Greenland with high thinning rates over a large area. When the entire ice sheet is considered, the differences can be expected to be smaller because not all outlet glaciers are thinning. To assess the difference in volume or mass changes between interpolation methods for the entire ice sheet will be part of future work.
Table 4. Volume Change Estimates for Jakobshavn Isbræa
Volume change estimates are shown in km3 yr−1 for individual years 1996–2007, as well as averaged over the periods 1999–2007 and 2003–2007.
 Furthermore, it would be useful to assess whether ST-KED is suitable for use on other outlet glaciers, or whether Jakobshavn Isbræ is particularly suited to this approach. An important requirement is that dynamically induceddH/dt values are dominant with respect to the SMB induced dH/dt. If that is not the case, it could be possible to use SMB from a regional climate model [e.g., Ettema et al., 2009] to correct for the influence of SMB. If ice dynamics dominate, the assumption that significant dynamic thinning occurs mainly in the high-velocity area would allow the use of KED, as the exact relation betweendH/dt and velocity (which is different for each glacier) is determined implicitly and is locally variable. In areas where spatial gradients in velocity are very low, the spatial trend in dH/dt will also be very small so that interpolated values will be very similar to those obtained by OK. Another requirement is that ice velocities are available at every location on the ice sheet. However, nearly all outlet glaciers are covered by InSAR measurements [e.g., Moon et al., 2012] and the gaps that occur are mainly located in the interior where gradients are small.
6. Summary and Conclusions
 Interpolation of dH/dt values derived from satellite altimeters is hampered by sparse sampling over, often, narrow outlet glaciers, where the largest changes occur due to dynamic thinning. To mitigate this problem, we have used velocity data from InSAR and speckle tracking as secondary data to inform the interpolation. The approach was tested over Jakobshavn Isbræ, where extensive airborne laser altimetry data (ATM) are available for validation.
 First, the relation between ice velocity and dH/dtwas examined, revealing a near-linear relation with highly negative correlation coefficients (<−0.85). Interestingly, when the glacier was thickening in the mid-nineties a correlation was found as well, suggesting the thickening had a ice-dynamical origin (i.e., a slowdown at some earlier time). Using ATM data as reference, interpolated spatial patterns of averagedH/dt are more realistic when velocity is used as a background field. Maximum thinning rates, that are about 7 m yr−1 for ordinary kriging, become similar to ATM data, at about 20 m yr−1, when velocity is considered in the interpolation. In addition, spatiotemporal kriging smooths effects of interannual variability in surface mass balance and it also improves estimates in periods with data gaps because it can use data from adjacent periods. Errors with respect to ATM (mainly located in the high-velocity area) are generally reduced significantly, except for some cases with very few available observations. For these reasons, interpolation with ST-KED was considered to yield the best results. However, very close to the grounding line, spatial gradients ofdH/dtbecome physically unrealistic due to a change in flow regime in this region. In addition, the transition from thickening to thinning (around 1998 in the high-velocity area) was delayed in our interpolated time series because it took time for the thinning to propagate into the areas with higher elevation where the ERS-2 radar altimeter was sampling.
 Volume changes that were calculated from the four resulting interpolation methods showed that volume loss derived from ST-KED was about 2–4 km3 yr−1, or 10–20%, higher compared to OK. When applied to the entire ice sheet ST-KED may therefore result in increased mass loss estimates. However, as Jakobshavn Isbræ is one of the largest and most dynamically active glaciers in Greenland, the significance of this percentage increase will likely be smaller and remains for further research.
 Previously, in order to estimate the Greenland mass balance, radar altimetry was merged with ATM data over outlet glaciers [Zwally et al., 2005]. The interpolation method proposed here yields similar results with radar altimetry only. Although for newer and future sensors such as CryoSat-II and ICESat-2 coverage over outlet glaciers will likely be less of a problem, our interpolation will improve ice sheet mass balance reconstructions from existing and past satellite altimeter data sets.
 This work was supported by funding to the ice2sea programme from the European Union 7th Framework Programme, grant 226375 and Ice2sea contribution 28. We thank the three reviewers and the editor for their thoughtful and constructive comments.