A leading-edge based method for correction of slope-induced errors in ice-sheet heights derived from radar altimetry

Satellite radar altimetry has been an important tool for cryospheric applications such as measuring ice-sheet height or assessing snow/ice anomalies (e.g., the extensive melt in Greenland in 2012). Although accurate height measurements are key for such applications, slope-induced errors due to undulating topography within the kilometre-wide pulse-limited footprint can cause multi-meter errors. Therefore different correction methods have been developed ranging from the slope method to the point-based method. Each of these methods have shortcomings as they either neglect the actual topography or the actual 5 footprint that can be estimated by a combination of the leading edge and topography. Therefore, a novel Leading Edge PointBased (LEPTA) method is presented that corrects for the slope-induced error by including the leading edge information of the radar waveform to determine the impact point. The principle of the method is that only the points on the ground that are within range determined by the begin and end of the leading edge are used to determine the impact point. Benchmarking of the LEPTA method to the slopeand point-based method based on CryoSat-2 LRM acquisitions over 10 Greenland in 2019 shows that heights obtained by LEPTA outperform the other methods when compared to ICESat-2 observations, both in the flat, interior regions of Greenland and in regions with more complex topography. The median difference between the slope-corrected CryoSat-2 and the ICESat-2 heights is almost negligible, whereas the other methods can have a 0.22 m and 0.69 m difference, and the Level-2 data provided by ESA have a 0.01 m difference. The median absolute deviation, which we use as an indicator of the variation of errors, is also the lowest in LEPTA (0.09 m) in comparison to the aforemen15 tioned methods (0.22 m and 0.13 m) and ESA Level-2 data (0.15 m). Based on that, we recommend considering LEPTA to obtain accurate height measurements with radar altimetry data, especially in regions with complex topography.


Introduction
Satellite radar altimetry is a key tool for assessing the status and dynamics of the cryosphere as it allows constructing digital elevation models (DEMs) (Slater et al., 2018), deriving height change of ice sheets (Helm et al., 2014a), understanding sea-20 sonal variations of snow (Adodo et al., 2018), and estimating snowpack properties (Lacroix et al., 2008). To obtain accurate information on heights, altimetry processing involves correction for instrument errors, atmospheric effects, tidal effects, and slope-induced errors (Bouzinac, 2012). Among the correction processes, correction for slope-induced errors has been of crucial importance as they can affect the results of height measurements significantly. For example, according to the error propagation in Brenner et al. (1983), the CryoSat-2 satellite at an altitude of 730 km, and measuring heights of a terrain with a 0.6 • slope, can give a vertical offset of approximately 40 m and a horizontal offset of 7.6 km.
To correct for the slope-induced errors, different methods have been developed (Brenner et al., 1983;Remy et al., 1989;Bamber, 1994;Roemer et al., 2007). The most widely used methods involve both a correction to the height as well as a relocation of the satellite measurement location from nadir to the expected impact point on the terrain. These correction methods are typically referred to as 'slope' and 'point-based' methods (Levinsen et al., 2016). The slope method assumes constant surface 30 parameters within the altimeter footprint and calculates the relocated longitude, latitude and height according to trigonometry (Brenner et al., 1983;Remy et al., 1989;Bamber, 1994). The point-based method takes the full height information within the satellite footprint and searches for the smallest range from the satellite to the terrain surface (Roemer et al., 2007).
Although both methods have been refined and applied with reliable results, they both show methodological shortcomings.
The slope method, for example, tends to ignore the local topography within the footprint and therefore may not be accurate 35 enough in undulating areas (Levinsen et al., 2016). The point-based method of Roemer et al. (2007), on the other hand, is more accurate in the undulating regions (Roemer et al., 2007;Levinsen et al., 2016) as it considers the detailed topography, but assuming a fixed footprint size neglects the actual footprint that illuminates the terrain. For example, by taking the averaged range within the assumed footprint, this method may ignore part of the terrain that actually contributes to the return signal, or assumes that part of the terrain not visible to the satellite could contribute to the return signal (See Fig. 1). The recent availability 40 of high-resolution DEM products, however, provides the opportunity to determine the part of the terrain contributing to the rise of the leading edge, and therefore helps determining the actual footprint of the radar altimeter.
To overcome the shortcomings of both methods, we present a novel Leading Edge Point-Based (LEPTA) method (Section 3) that exploits high resolution DEM information to correct for the slope-induced error by including the leading edge information of the radar waveform to determine the impact point. The principle of the method is that only the points on the ground that are 45 actually within range determined by the begin and end of the leading edge are used to determine the impact point.
The paper is organised as follows. The data used for radar altimetric processing and result assessment are described. Then, the different methods used for correction and the assessment workflow are introduced. The results and analysis/discussion will then be presented.
2 Data and pre-processing 50 To assess the performance of the LEPTA method, we apply it to all CryoSat-2 LRM acquisitions over Greenland in 2019 and benchmark it to the slope-and point-based methods by comparing it with laser altimeter ICESat-2 data.

CryoSat-2 observations
On the interior of Greenland ice sheet, data acquired by CryoSat-2 are in Low Resolution Mode (LRM). LRM is the conventional pulse-limited mode that requires correction for slope-induced errors, and the footprint of this mode is approximately 1.65 km in diameter (Bouzinac, 2012). Our evaluation employs all acquisitions from Jan.-Dec. 2019, resulting in approximately 2.4 × 10 6 valid acquisitions, in order to ensure abundant spatial and temporal coverage.
To process the waveform information and obtain height estimations, Level-1b (L1b) (European Space Agency, 2019a) waveforms were retracked using the offset centre of gravity (OCOG) method (Wingham et al., 1986) documented in Bamber (1994, because of its precision and robustness (Bamber, 1994;Schröder et al., 2019). According to Davis (1997), a 10% 60 threshold is ideal for detecting ice-sheet height change, a 20% threshold is the most proper for estimating the absolute icesheet height, and a 50% threshold is the most appropriate for estimating the absolute height when the waveform is dominated by surface scattering e.g. Antarctic ice shelves. A 10% threshold is also applied when volume scattering is strong (Aublanc et al., 2018). In this study, a 20% threshold is selected to obtain accurate absolute height estimations of the interior of the Greenland ice sheet, as a compromise between surface scattering and volume scattering (Aublanc et al., 2018). The 20% threshold is close to the 25% threshold use by ESA L2 processing (Bamber, 1994). This is also a more realistic threshold when the performance is assessed with ICESat-2, which measures the snow-air surface. In addition, in this process, waveforms that do not have a distinguishable noise and beginning of leading edge are dropped.
Additionally, Level-2I (L2I) height data obtained with the OCOG retracker from European Space Agency (2019b) were used as benchmark dataset. In the L2I products the slope-induced error is corrected with the Helm et al. (2014b) DEM, of which the 70 resolution is 1 km×1 km (Helm et al., 2014a).

ArcticDEM
Within the slope correction methods, a reference DEM is required to determine the impact point. The slope method therefore uses a low resolution DEM (or a downsampled version) as it assumes a constant slope within the radar footprint. On the contrary, the point-based methods (i.e., LEPTA and the point-based method proposed by (Roemer et al., 2007)) require DEMs 75 with higher resolution, to provide the full information of the local terrain.
In this study, ArcticDEM was used as reference DEM as it is constructed from recent stereo satellite imagery and is available in high resolution (2 m×2 m) (Porter et al., 2018). The systematic error of ArcticDEM is less than 5 m (Noh and Howat, 2015) and the DEM has been updated since 2016. ArcticDEM was downsampled to 900 m resolution for the slope-based method and to 100 m resolution for the point-based and LEPTA methods used as a compromise between computational efficiency and the 80 demand for high resolution. We also downsample it to various resolutions (200 m to 900 m with a 100 m interval) to assess the impact of DEM resolution affects on the correction methods (Subsection 3.3).

ICESat-2 observations
For validation of the different slope correction methods, ICESat-2 L3A Land Ice Height products (Smith et al., 2020a) were used as they provide independent accurate processed height measurements. ICESat-2 uses the Advanced Topographic Laser

85
Altimeter System (ATLAS) which emits green light pulses and counts the received photons (Abdalati et al., 2010). The laser beams are configured in a 2 × 3 array. The distance between and within beam pairs is ∼3.3 km and ∼90 m, respectively . The along-track resolution of land ice height products is ∼20 m (Smith et al., 2020b). The ICESat-2 data are available at the National Snow and Ice Data Center (NSIDC) website (https://nsidc.org/data/atl06).

Slope correction methods
The geometry of different slope-induced error correction methods is briefly illustrated in Fig. 1.

Slope-based correction method
The slope method uses the slope of the low resolution DEM at the nadir point to compute the impact point. As such it assumes that the slope within the CryoSat-2 footprint is constant, and is defined by direction θ and magnitude Φ (Cooper, 1989;Bamber, 95 1994). The corrected height, represented by h C , is the height of the impact point P s , can then be obtained by (Bamber, 1994): where Γ is the central angle between the satellite and P s , R I is the range between P s and the centre of the curvature at P s , and R s is the range between the centre of the curvature and the satellite at latitude and longitude λ. R α represents the radius of Earth's curvature at P s . The corrected location of the impact point in latitude φ C and longitude λ C (in radians) are computed where X and Y define the position of P s in Cartesian coordinates and Application of the slope-method in Fig. 1 shows that the impact point will be assumed at the position P s . Inaccuracies usually occur when this method is applied to complex terrains, due to the simplification of the complex topography to a constant slope

Point-based correction method
The point-based method directly uses the topographic information from the a-priori DEM to find the impact point (P p ). It does so by minimising the mean distanceR P to the satellite over a pre-defined fixed-size rectangular footprint area (e.g., 2 km×2 km in Roemer et al. (2007)). Assuming the pre-defined rectangular footprint (area A) consists of n DEM grid cells,R P is 115 computed by (Roemer et al., 2007): where A P j andR P j are the area and range of each grid cell j. The position (defined by φ c and λ c ) of the footprint that minimisesR P is then obtained as P p . The range between the satellite and P p is referred to as r p . The corrected height h C is computed as (Roemer et al., 2007) 120 where h N is the surface height of the nadir point relative to the reference ellipsoid (i.e., h S minus the retracked range), h S is the ellipsoidal height of the satellite, and h I is the DEM height of P p . It also shows, however, that this approach can take DEM points into account that actually do not contribute to the rise of the leading edge (i.e., points that fall outside the footprint).

125
The LEPTA method is similar to the point-based method as it also uses the topographic information from the a-priori DEM to find the impact point (P l ), but differs in the search method of the impact point. Instead of pre-defining a fixed rectangular pulse-limited footprint size, the LEPTA-method identifies the parts of the terrain that contribute to the rise of the leading edge.
To identify these parts, we define the beginning of the leading edge as the point where the normalised waveform power (values are between 0 and 1) is greater than 0.05. The end of the leading edge is more difficult to define as there might be multiple 130 peaks before the waveform reaches its maximum power. Here, we defined it as the point that is located at a distance ∆r from the range obtained by applying the OCOG retracker mentioned in Subsection 2.1. In this study, we used ∆r =3.5 m. The robustness of the results regarding the choice of ∆r will be further assessed in Subsection 4.3.
For each point, the distances are computed between the satellite and all DEM grid points within an area of 8 × 8 km centred around the nadir point. Thereafter, the DEM points are identified for which the range is within the interval defined by the 135 beginning and end of the leading edge. In case no DEM grid points are identified, the interval is adjusted by the minimum difference between the computed distances and the retracked range. Next, the location of P l is computed as the average of all K identified DEM grid points. Finally, the corrected height h C is computed by where h i DEM and r i DEM are the height of and the range between the satellite and DEM grid point i.
One of the advantages of the LEPTA method compared to the point-based method is that it includes points that contribute to the rise of the leading edge signal but outside the fixed footprint.
The difference in footprint and impact points between LEPTA-method and point-based method is illustrated in an example in Fig. 2. The brown areas indicate the areas on the surface that (is assumed to) contribute to the return signal for both the point-based and LEPTA method. Contrary to the point-based method, for LEPTA this area may take any shape. This example 145 is an illustration of the theoretical advantage of LEPTA. Assessment of the performance of different methods will be described in Subsection 3.2.

Performance assessment
To assess the performance of the LEPTA method, we benchmark the different methods by comparing their accuracy relative to reference data. First, we directly compare the corrected height measurements for each method with the reference height requires weighting functions based on the surrounding points. The differences between the CryoSat-2 and ICESat-2 heights are referred to as ∆h. Similar as before, we compute for each tile the median and median absolute deviation of ∆h.
When benchmarking the methods, two aspects of accuracy are assessed. First, we determine the difference between the slopecorrected CryoSat-2 measurements and the reference heights (h DEM or h ICE2 ) by means of standard statistical parameters (median, median absolute deviation, mean, and standard deviation). Second, we assess the spatial differences between the methods. The 165 statistical parameters are computed with and without removing outliers. Here, we consider any difference between h C and h DEM /h ICE2 outside the 10th and 90th percentiles as an outlier.

Sensitivity analysis
The LEPTA method is potentially sensitive to i) the way the 'end' of the leading edge is defined (which in turn determines the satellite footprint), ii) a bias in the DEM, and iii) the resolution of the used DEM. Another aspect that may impact the results 170 is the adopted OCOG threshold (see Section 2.1). To assess how our choices impact the results, we conducted a number of sensitivity analyses in which we: -Varied ∆r from 2-5 m in steps of 0.5 m.
-Added a bias to the DEM of -7.5-2.5 m in steps of 2.5 m.
-Varied the DEM resolution from 200-900 m in steps of 100 m.

Comparison with ArcticDEM
Benchmarking the different methods to the 100 m ArcticDEM (Table 1) shows that outliers have an important impact on the performance of all methods with multi-metre outliers for the 1st, 5th, 95th and 99th percentiles. We see moreover that the ESA 180 L2I and slope method results include larger outliers than the point-based and LEPTA methods. Since these outliers have a large impact on the mean and standard deviation, we repeated our analysis with all values outside the 10-90th percentile interval removed.
Removing the outliers significantly reduces the standard deviation of h C − h DEM for all methods and brings the mean closer to the median. Comparison of the mean and median values (Table 1)  and median, thus are the least ideal. An additional note is that the mean and median from all methods become above 0, which means that the heights obtained by these methods are generally higher than ArcticDEM heights. It is also noticed that the ESA L2I products result in more valid outputs than the self-implemented methods.
Comparison of the spatial patterns of median and median absolute deviation (Fig. 3) shows large spatial differences in both pattern and magnitude between the different methods. In general, the largest median and median absolute deviation values  Table 1 and  Fig. 3, it can be concluded that LEPTA has the most ideal performance, when compared with ArcticDEM.

Validation with ICESat-2 observations 205
Comparison of CryoSat-2 and ICESat-2 heights (Fig. 4) shows again the importance of outliers on the slope-correction methods although the outliers are typically lower than for the ArcticDEM comparison. For all methods, the ∆h values within the visualised range in the probability distribution plot have a longer tail on the right side than the left side of the median. Very large outliers (1st and 99th percentiles) occur at -39.7 m for ESA L2I and at 22.72 m for the slope method.  Comparison of the height differences between CryoSat-2 and ArcticDEM and ICESat-2, respectively, shows moreover that the height differences with ICESat-2 are smaller, probably due to the longer time gap between CryoSat-2 and ArcticDEM, as satellite imagery data for generating ArcticDEM were gathered since 2007 (Noh and Howat, 2017;Howat et al., 2019), whereas ICESat-2 measurements were obtained in the same month as CryoSat-2 data.
The spatial distribution of the mean and median absolute deviation of ∆h relative to ICESat-2 ( Fig. 4) shows clear spatial 220 patterns for ∆h, where the LEPTA method again outperforms the other methods with lower height differences compared to the point-based method (i.e. median values of 0 versus 0.69 m, respectively) and more spatially homogeneous patterns compared to ESA's L2I products and the slope method. The latter show for example large errors over the steeper areas near the edges.
To sum up, the point-based method generally has the highest median values, while the slope method has the largest contrast between the interior and the margins. To assess the sensitivity of the methods to biases in the DEM (e.g. due to surface mass balance changes), we performed the analysis introduced in Subsection 3.3 with ICESat-2 as validation data. Figure 6 shows that the slope and the point-based methods are not affected by biases in the DEM, while for LEPTA the median changes about 10 cm when the bias in the DEM is 2.5 m. The median absolute deviation is also around 3 cm higher. This is understandable as the presence of a bias in the DEM does not affect the slope or the relative differences between the DEM points, which are the key information to the slope 240 method and the point-based method respectively. However, in the case of LEPTA, when the DEM heights are biased, the search range is different from using the heights of the original DEM heights. Using the same search range may ignore some points on the terrain that actually contribute to the waveform leading edge, resulting in a slightly larger bias. This result indicates that LEPTA is relatively sensitive to the bias in DEM heights, compared to other methods.   The DEM resolution also has little impact on the method, although a high-resolution DEM is recommended. On the contrary, the choice of the retracker affects all the methods to a similar extent. This analysis also agrees with the study of Brenner et al. (2007) that the difference between radar and laser altimeters is retracker-dependent.

275
By showing the importance of impact points over steeper margin areas, our results confirm earlier work of Levinsen et al. (2016) in the margin regions, where they also showed that the point-based method outperforms the slope based methods. The improved performance of the point-based method and LEPTA methods can be explained by assumption of a constant slope within the footprint in the slope-based method, which results in a biased impact point further away from the satellite than the optimal location (Levinsen et al., 2016). An explanation for this improved performance of LEPTA over the point-based method 280 can be found in the design of the method which only takes into account areas that are within the footprint of CryoSat-2 ( Fig. 1 and Fig. 2).
Comparison of the ESA L2 products with the slope-based correction method shows finally that the ESA product outperforms standard slope correction method. This agrees with Levinsen et al. (2016) who attributed the different performance between the ESA L2 products and the self-implementation to the Doppler slope correction implemented in ESA L2 products (Blarel 285 and Legresy, 2012) and DEM differences. The ESA L2 products also provides more valid outputs, compared to the selfimplemented methods. This could also be attributed to the Doppler slope correction as it has the advantage of being valid directly for all continental surfaces (Blarel and Legresy, 2012). On the contrary, the self-implemented methods dropped the waveforms that are not similar to the standardised waveform (as in Fig. 1 of Simonsen and Sørensen (2017)). LEPTA also leaves out the points that do not have valid DEM values within the search range. More invalid points are therefore dropped in 290 this process. This phenomenon leaves LEPTA potential for improvement.

Conclusions
Reducing slope-induced errors is a key correction algorithm when applying the radar altimeter data. To correct for this error, different methods have been developed to determine the impact point, which all rely on footprint assumptions: e.g. slopemethod, which assumes a constant slope within the footprint, or the point-based method, which assumes a fixed footprint size and defines the reflecting point as the shortest mean range of points within each assumed footprint. Each of these methods has shortcoming as they either neglect the actual topography or the actual footprint that can be estimated by a combination of the leading edge and topography.
To overcome this shortcoming, we present a novel Leading Edge Point-Based (LEPTA) method that corrects for the slopeinduced error by including the leading edge information of the radar waveform to determine the impact point. The principle of 300 the method is that only the points on the ground that are within range determined by the begin and end of the leading edge are used to determine the impact point. This requires the assistance of a high-resolution DEM, e.g. 100 m resolution.
Statistics show that the LEPTA method outperforms all other methods with the smallest median and variability of errors. The median is almost identical to the ICESat-2 height measurements. Spatially, LEPTA has a good improvement compared to the traditional slope method on the margins of the ice sheet. Therefore, LEPTA is a method worth considering to obtain accurate 305 height measurements with radar altimeter, especially in regions with complex topography.
Author contributions. WL conducted data management, processing and analysis, produced the figures and provided the manuscript with contributions from all co-authors. CS designed the study and provided expertise and software for radar-altimetry processing. SL provided support on statistical analysis and data visualisation.
Competing interests. The authors declare they have no conflict of interest.