“ Importance of slope-induced error correction in volume change estimates from radar altimetry ”

In deriving elevation change rates (d H /dt) from radar altimetry, the slope-induced error is usually assumed to cancel out in repeat measurements. These measurements, however, represent a location that can be significantly further upslope than assumed, causing an underestimate of the basin-integrated volume change. In a case-study for the fastflowing part of Jakobshavn Isbræ, we show that a relatively straightforward correction for slope-induced error increases elevation change rates by up to several metres per year and significantly reduces the volume change error with respect to laser altimetry for the area of interest.


Introduction
Several sources of uncertainty affect measurements of ice sheet surface elevation derived from satellite radar altimetry (SRA; e.g.Brenner et al., 2007;Bamber, 1994).Potentially the largest one, referred to as the slope-induced error (Brenner et al., 1983), is caused by regional surface slopes.For narrow, fast-flowing outlet glaciers, the lowest surface may not be sampled at all (see Fig. 9 of Thomas et al., 2008), however these are often the areas that show the largest elevation changes (Pritchard et al., 2009).The radar return signal does not originate from the point directly underneath the satellite (nadir), but from the closest point to the satellite, which can be significantly displaced upslope from nadir.For a 1 • slope and a satellite altitude of 800 km, which is not unusual at the ice sheet margin, the displacement between nadir and the actual measurement location is about 14 km, and the vertical error is about 120 m.Using a regional slope estimate, it is relatively straightforward to relocate the measurement to its correct location (Bamber, 1994).When elevation change is concerned, however, it is usually assumed that the slope-induced error remains constant and the effect cancels (Thomas et al., 2008).While the error in the vertical is indeed the same for repeating elevation measurements, the location of the measured elevation change rate will still be displaced from its true position.Here, we show that the integrated volume change can be significantly underestimated, because elevation changes that are measured at an upslope location are incorrectly located closer to the margin.

Data
To demonstrate the effect of correcting for slope-induced error, we use data from the radar altimeter (RA-2) on ESA's Envisat satellite, that was launched in 2002.It continues the SRA time series from ERS-1 and ERS-2 and is in a similar orbit with an altitude of about 800 km, a repeat period of 35 days, and a latitudinal coverage up to 81.5 • .We use Envisat cross-over clusters (Li and Davis, 2008) from which we derive average elevation change rates for [2003][2004][2005][2006].The selected study area is the fast-flowing region of Jakobshavn Isbrae, Greenland's largest outlet glacier, located on the southwest coast.Since about 1998, it has been accelerating and thinning significantly (Joughin et al., 2008), and has been densely surveyed by airborne laser altimetry (ATM; Krabill et al., 2004).
We use airborne (ATM) and spaceborne (ICESat) laser altimetry data as a validation dataset.Elevations from these data sources do not suffer from the same slope-induced error as they have a small footprint with known pointing: ≈60 m for ICESat (Zwally et al., 2002) and 1-2 m for ATM (Krabill et al., 2002).Elevation change rates for 2003-2006 (consistent with Envisat) were derived from ICESat using the "plane" method (Howat et al., 2008).A plane is fitted through data from several near-repeat tracks.For each plane, which is typically about 700 m long and a few hundred metres wide, two-directional slopes and a temporal elevation change rate dH /dt are fitted using multivariate linear regression (Moholdt et al., 2010).A regression is only performed if a plane has at least 10 points from four different tracks that span at least a year.Prior to the regression, outliers (outside 2σ ) are removed, and only elevation changes with an associated standard error on dH /dt of less than 0.50 m yr −1 are considered.ATM footprints were preprocessed and averaged over two 150 m platelets on each side of the plane (Krabill et al., 2002).We used a similar method as outlined above to derive 2003-2006 elevation change rates from all available flight lines, but instead of platelets, 1 km pixels were used.
Other datasets that are employed here are a 1 km digital elevation model (Bamber et al., 2001) from which slope and aspect are derived, and an ice sheet velocity mosaic used for interpolation.The velocity field was derived by Joughin et al. (2010) from a combination of radar interferometry and speckle tracking using RADARSAT-2 data from the winters 2000-2001, 2005-2006, and 2007-2008.Figure 1 shows the slope, as well as locations of Envisat cross-over points, ICE-Sat tracks, ATM flight lines, and velocity contours for the fast-flowing part of Jakobshavn Isbrae.

Methodology
The slope-induced error is schematically illustrated in Fig. 2a by a range measurement R to an inclined surface with slope α and aspect β.The measurement location is displaced from nadir by a horizontal distance D. Three correction methods for the range exist (Bamber, 1994).First, the di-rect method corrects the range measured at nadir (corrected range R c = R/ cos(α)).The second method, the relocation method, corrects R to R c , where R c is now the range to the point closest to the satellite (now R cos(α)), and displaces the location by R sin(α).The intermediate method, finally, finds the location where R = R c , and relocates the measurement to that point (Remy et al., 1989).When one is interested in elevation changes, errors in the vertical cancel out, but measurements are still located at the wrong locations.This can lead to underestimation of area-integrated volume changes as dH /dt values obtained by radar altimetry are actually located further upslope than assumed (Fig. 2).For cross-over analyses, where dH /dt is derived for locations where ascending and descending tracks cross each other, a two-dimensional correction should be applied.The displacement D is given by D = E sin(α) cos(α), where E is the satellite altitude, equivalent (but not necessary equal) to R c in Fig. 2a, and 800 ± 20 km for Envisat.Sensitivity to variations in the orbit altitude is small: a sensitivity experiment indicated that, at 1 • slope, a 40 km (5 %) range in altitude causes a ≈5 % range in horizontal displacement with respect to the displacement for an altitude of 800 km.Sensitivity to slope angle is about the same (a 5 % range in slope causes a ≈5 % range in horizontal displacement); therefore ideally contemporaneous estimates of slope should be used.The direction of displacement is opposite to the slope aspect β, which is the direction of the steepest downward slope.If β is defined as 0 radians for north and increasing clockwise to 2π radians, the relocation in x and y directions are given by dX = D sin(β − π ) and dY = D cos(β − π ).Although small-scale undulations can also cause error (Bamber and Gomez-Dans, 2005), we use average slope and aspect over a 100 km 2 area centered on the nadir location to correct for regional slope, as 10 km is the approximate length scale of the expected displacement.
We illustrate the effect of the correction on volume change using a hypothetical testcase.We simulated a 100 × 100 km surface with a slope increasing linearly toward lower elevations from 0.5 to 1.5 degrees (Fig. 2b).Synthetic dH /dt data, ranging linearly from 2 m yr −1 at 1.5 • slope to 1 m yr −1 at 0.5 • slope, are evenly spaced at 10 km intervals.dH /dt data coverage is assumed to extend linearly beyond the domain, so the correction displaces data points "into" the domain.A full 1 km resolution dH /dt field is obtained using inverse distance interpolation.In Fig. 2c, all data locations are corrected for slope-induced error, enlarging the area with the largest dH /dt.For this particular testcase, the relative difference in volume change (i.e. between Fig. 2b and c) is 10.4 %.It should be noted that dH /dt values are relatively modest compared to Jakobshavn Isbrae.

Results
Figure 3a shows a scatterplot of dH /dt from Envisat cross-over clusters versus dH /dt from interpolation of ATM/ICESat, for both the uncorrected and the corrected Envisat locations.Interpolation was conducted using kriging with external drift (KED), which uses the spatial pattern of ice velocity as a proxy for that of dH /dt, i.e. a steep spatial gradient in velocity leads to a similarly steep gradient in dH /dt in the absence of local dH /dt measurements.The method is described in detail elsewhere (Hurkmans et al., 2012), where they found that for Jakobshavn Isbrae, KED results in more realistic dH /dt patterns (with respect to ATM) than other methods investigated.The sparsity of Envisat data is illustrated by the fact that there are only 23 Envisat crossover clusters in the study area.Uncorrected values are generally corrected towards the ATM/ICESat values effectively, sometimes with elevation change corrections of several metres per year.There is still, however, considerable noise in the corrected scatterplot, because (i) the correction only corrects Envisat data for the footprint-average slope and not for smaller scale undulations, (ii) interpolated values from ATM/ICESat were used because Envisat and ATM/ICESat footprints never exactly overlap, and (iii) various differences between radar and laser altimetry, such as footprint size, orbit errors, and the backscattering correction necessary for radar altimetry.The effectiveness of the correction is illustrated by the correlation coefficient which increases from 0.35 to 0.88 after the correction for slope-induced error.
Interpolated dH /dt values are shown in Fig. 3b and c.In Fig. 3b, a transect is shown constructed by calculating the average north-south dH /dt within the 300 m yr −1 velocity contour for each 1 km pixel moving east from the grounding zone.The difference in thinning rates between corrected and uncorrected Envisat data increases from about 0.4 m yr −1 at 80 km from the grounding zone to about 2.5 m yr −1 at 5 km.A three kilometre zone adjacent to the presumed grounding line was not taken into account because of uncertainty in its location (Hurkmans et al., 2012).After correction for slope-induced error, interpolated thinning rates are both larger and more widespread.This can be seen in between the 100 and 300 m yr −1 velocity contours (Fig. 1), with a more complete Envisat coverage, the equivalent error reduces from −55 % to −44 % (with respect to a volume loss of 6.6 km 3 yr −1 for ATM/ICESat).One explanation for the much larger volume changes for ATM/ICESat compared to Envisat is the better sampling by ATM/ICESat of regions affected by increased surface melt.This is not captured by the KED interpolation, which only accounts for dynamically induced dH /dt (Hurkmans et al., 2012).

Fig. 1 .
Fig. 1.Slope from Bamber et al. (2001), plotted over a MODIS image from June 2002.Also uncorrected (large black dots) and corrected (black triangles) locations of Envisat cross-over clusters, ICESat tracks (red dots), and ATM flight lines (small black dots) are shown.Velocity is shown as 100, 300, and 1000 m yr −1 contours.
Fig. 2. (A) Geometry of the two-dimensional slope-induced error correction.α and β are, respectively, the slope and aspect (both in radians), R is the range (≈800 km), R c the corrected range and R the vertical correction.D is the horizontal distance over which a cross-over point is displaced, and dX and dY are its components in the x-and y-direction.(B) shows dH /dt on a synthetic slope without correction for slope-induced error, and (C) with correction.Black dots show dH /dt data locations without (B) and with (C) correction.

Fig. 3 .
Fig. 3. Comparison of dH /dt from Envisat, before and after slope correction.(A) shows dH /dt from Envisat cross-over clusters versus interpolated dH /dt from ATM/ICESat at the corresponding locations.Red dots show uncorrected, blue dots corrected data.ρ is the Pearson correlation coefficient for both cases.(B) shows interpolated dH /dt along a transect inland from the grounding zone.All values along the transect are calculated as the north-south average within the 300 m yr −1 velocity contour, and the asterisks indicate the (corrected and uncorrected) location of the only Envisat measurement in this area (Fig. 1).Also dH /dt from ATM/ICESat and velocity are shown.(C) shows the difference between the interpolated dH /dt based on Envisat, before and after correction.Negative values indicate higher values after correction, and velocity is shown as 100, 300, and 1000 m yr −1 contours.