1 Introduction

Satellite altimetry is one means of determining mass changes in ice sheets (Shepherd et al. 2018), which are affected by climate change and affect the global sea level. Mass changes are derived from estimations of volume changes combined with firn and ice densities (Shepherd et al. 2012; Hurkmans et al. 2014; Khan et al. 2015; McMillan et al. 2016).

Deriving height changes from satellite altimetry usually involves two steps: local height change determination from repeat altimetry, and subsequent spatial interpolation at unobserved areas, possibly involving smoothing. As the repeated altimeter measurements do not refer to exactly the same position, the local topography has to be accounted for. The repeat altimetry analysis (RAA) approach has been widely used to solve this problem (Legrésy et al. 2006; Flament and Rémy 2012b; Nilsson et al. 2016; Sørensen et al. 2018a; Schröder et al. 2019). However, it is subject to a number of methodological choices in the processing chain. They include the size and arrangement of RAA cells, modeling of seasonal height changes (Sørensen et al. 2011), the use of signal parameters such as leading-edge width or backscatter to model changing signal penetration (Simonsen and Sørensen 2017), and outlier elimination.

The spatial coverage of height change estimates from RAA depends on the orbit geometry and the success of RAA estimates, and is neither homogeneous nor complete. Therefore, subsequent interpolation and filtering is commonly applied. Hurkmans et al. (2012b, 2014) apply ordinary kriging (OK) with spatiotemporal modeling of the underlying process. They also introduce external data of higher resolution to improve the interpolation via external drift. OK is an exact interpolator (Cressie 1993). Exact interpolation is desired when the input values are free of error. This is not the case for height changes derived by RAA. They are the output of a fitting algorithm and are provided with individual uncertainties of the estimate. Uncertainties can be used to refine interpolation, for example, by using heterogeneous measurement-error-filtered kriging (HFK) (Christensen 2011).

Note that the preprocessing of altimetry data, which is not the subject of this investigation, can differ in slope correction and radar waveform retracking algorithms (Hurkmans et al. 2012a; Helm et al. 2014; Nilsson et al. 2016; Sørensen et al. 2018b). This affects positioning and height measurements as well as further derived values such as height change estimates. As a matter of fact, the irregular spatial coverage remains unaffected, and coping with it is the focus of this study.

Various approaches in RAA and interpolation lead to differences in the final height change estimates and their correspondingly derived uncertainties. In RAA, altimetry measurements are jointly processed in a defined area (called a cell). The effect of cell size and related modeling of the underlying local topography on trends in height change is investigated. For the subsequent interpolation OK, inverse distance weighting (IDW), filtered kriging (FK) and HFK are applied. These four different interpolation approaches are investigated with respect to the accuracy of the interpolated height changes and the reliability of their uncertainty estimates.

The coastal areas show the highest changes in the elevation of the Greenland Ice Sheet (GIS) (Sørensen et al. 2018a). These areas are covered by CryoSat-2 measurements in interferometric synthetic aperture radar (SARIn) mode. While the low-resolution mode (LRM) has a beam-limited footprint of 20 km and a pulse-limited footprint of 1.5 km, SARIn mode enhances along-track resolution to 300 m and uses the across-track angle for measurement attribution (Wingham et al. 2006). Obtaining full spatial coverage of height change in an area with sloped topography in combination with the irregular spatial distribution of RAA height change estimates is challenging. This simulation study is conducted at the eastern margins of the Northeast Greenland Ice Stream (NEGIS), where these circumstances can be reproduced without incorporating additional difficulties such as narrow fjords.

Section 2 introduces the mathematics of RAA and the interpolation algorithms used. Section 3 describes the simulations, including the synthetic data set and processing details. The results are analyzed and discussed in Sect. 4. Section 5 summarizes the results and highlights the effect of RAA parameter choices and the benefits of HFK for subsequent interpolation.

2 Theory

2.1 Repeat Altimetry Analysis

RAA is a fitting method based on least-squares regression. It uses height measurements \(h_i\) and the corresponding locations and times in a cell to estimate parameters that describe the underlying spatial and temporal variation in the measured elevations.

$$\begin{aligned} h_i = \text {Ft}_i + \text {Fl}_i + \text {Fs}_i + \text {res}_i. \end{aligned}$$
(1)

The components \(\text {Ft}_i\), \(\text {Fl}_i\) and \(\text {Fs}_i\) describe the dependence on time, location and radar signal return characteristics, respectively. Parameters of these components (as specified below) are estimated by RAA. The resulting residuals between model and observation are depicted in \(\text {res}_i\). The final choice of the RAA parameter set depends on the satellite track configuration and the measurement properties of the mission used.

The design of the cells can vary (Sørensen et al. 2018a). If the cells are arranged along the subsatellite tracks, rectangular cells that span several consecutive shots along track and the repeat corridor across track are well established (Legrésy et al. 2006; Ewert et al. 2012). In the case of regular grids, both rectangular (McMillan et al. 2014) and circular (Simonsen and Sørensen 2017) cell shapes are commonly used. The latter permit a constant maximum measurement distance to the cell center and are used in this article.

In this study the time dependence in Eq. (1) is characterized by a linear trend \(\frac{\mathrm{d}h}{\mathrm{d}t}\).

$$\begin{aligned} \text {Ft}_i = \frac{\mathrm{d}h}{\mathrm{d}t}(t_i-t_0). \end{aligned}$$
(2)

Additionally, seasonal elevation variations modeled by a combination of sine and cosine terms can be introduced here (Sørensen et al. 2011).

The location-dependent component models the local topography inside the analyzed cell. A common approach is the fit of a plane (Smith et al. 2009; Sørensen et al. 2015; Schröder et al. 2019)

$$\begin{aligned} \text {Fl}_i = a_0 + a_1 x_i + a_2 y_i, \end{aligned}$$
(3)

where \(x_i\), \(y_i\) are horizontal Cartesian coordinates with their origin in the cell center, and \(a_0\), \(a_1\), \(a_2\) are the parameters of the plane. Other local topography models exist, such as the biquadratic model with six parameters,

$$\begin{aligned} \text {Fl}_i&= a_0 + a_1 x_i + a_2 y_i\nonumber \\&\quad + a_3 x_i^2 + a_4 y_i^2 + a_5 x_i y_i, \end{aligned}$$
(4)

used by Nilsson et al. (2016), Simonsen and Sørensen (2017) or the nine-parameter model

$$\begin{aligned} \text {Fl}_i&= a_0 + a_1 x_i + a_2 y_i\nonumber \\&\quad + a_3 x_i^2 + a_4 y_i^2 + a_5 x_i y_i \nonumber \\&\quad + a_6 x_i^2 y_i + a_7 x_i y_i^2 + a_8 x_i^2 y_i^2, \end{aligned}$$
(5)

used by Ewert et al. (2012), Wouters et al. (2015). In contrast, the local topography model is reduced to only one parameter, \(a_0\), when a digital elevation model (DEM) is subtracted beforehand (Sørensen et al. 2011; Helm et al. 2014; Simonsen and Sørensen 2017).

Additional parameters may be useful depending on specific characteristics of the altimeter return signal, which is affected by characteristics of the reflecting surface and volume of firn or ice. For CryoSat-2, by the parameter \(dBS\), the effects of time-variable signal penetration with anomalies of backscattered power \(bs_i-\overline{bs}\) are described.

$$\begin{aligned} \text {Fs}_i = dBS(bs_i-\overline{bs}). \end{aligned}$$
(6)

This modeling can be further expanded by involving leading-edge width or trailing-edge slope of the signal waveform (Flament and Rémy 2012a; Simonsen and Sørensen 2017), or a bias between ascending and descending satellite tracks (McMillan et al. 2014; Simonsen and Sørensen 2017). However, such interactions between the radar measurements and the uppermost firn layer are complex (Simonsen and Sørensen 2017; Adodo et al 2018) and are not yet fully understood. This study focuses on the analysis of different aspects of spatial sampling. In regions where the topography is sufficiently flat to allow for reliable analysis of such waveforms, the influence of the scattering characteristics acts on significantly larger scales. Therefore, we expect a negligible influence of such parameters on our results and do not further analyse these types of parameters.

For each cell, all selected parameters are solved for simultaneously in a joint least-squares adjustment to all measurements \(h_i\) lying inside this cell. The elevation change parameter \(\frac{\mathrm{d}h}{\mathrm{d}t}\) is the target of this analysis. The RAA approach also provides an a posteriori standard error for \(\frac{\mathrm{d}h}{\mathrm{d}t}\), based on the statistics of the residuals \(\text {res}_i\).

The choice of the local topography model depends on several considerations. The actual topography of Greenland, which is smooth in the ice sheet interior and rugged at the margins, has to be taken into account. The ability to properly model the topography is also linked to the cell size and the number of observations available. Because of the limited number of observations in smaller cells, the number of topography parameters is restricted, while in larger cells a simple model might not be able to depict the actual local topography. The effect of different cell sizes and local topography models is investigated in this study.

2.2 Interpolation

The estimated height changes \(\frac{\mathrm{d}h}{\mathrm{d}t}\) have to be interpolated to obtain values at places where data are missing or RAA could not be solved successfully. The interpolation methods used in this study are all based on the same principle: observations \(Z\) at locations \(r_i\) are used to calculate a new value \(Z^*\) at a certain location \(r_0\) by weighted summation (Myers 1991; Cressie 1993; Chilès and Delfiner 2012). In this section, \(r\) denotes a two-dimensional position vector.

$$\begin{aligned} Z^*(r_0)=\sum _{i=1}^n \lambda _{i} \cdot Z(r_i). \end{aligned}$$
(7)

To derive the weights \(\lambda _{i}\), IDW uses only geometric information, while kriging uses a geostatistical approach. Different kriging methods have been developed, based on formulations by D. Krige and G. Matheron (Chilès and Delfiner 2012). This study focuses on OK, FK and HFK. Detailed information about the used kriging methods are given for example by Cressie (1993), Chilès and Delfiner (2012), Christensen (2011).

The number of points used for interpolation depends on the data set and the user’s decision. Irregularly distributed observations may lead to distorted results, for example due to unwanted screening (Cressie 1988; Chilès and Delfiner 2012) or spatial biases. Therefore, the surroundings are often divided into a certain number of sectors, selecting observations in each of the sectors to obtain a more uniform distribution. The same applies for variograms (Stosius and Herzfeld 2004).

In the field of geodesy, besides kriging, different least-squares collocation methods are commonly applied, which may also model uncertainties (Nilsson et al. 2015; Sørensen et al. 2018a). These include different treatment of measurement and interpolation errors. This study focuses on interpolation error, although methods such as IDW and OK would need additional uncertainty assessment. Basic agreement between kriging and collocation methods is confirmed (Dermanis 1984), although the individual requirements can differ. Distinguishing these two statistical approaches is outside the scope of this article.

Previous studies prove the general suitability of the selected methods (Rühaak 2015; Christensen and Berrett 2016; Kang et al. 2017) and provide some comparison results (Chaplot et al 2006; Li and Heap 2011). To the authors’ knowledge, HFK has not yet been applied to height changes derived by satellite altimetry.

2.2.1 Inverse Distance Weighting

The weights of IDW depend on the normalized distances \(d_i\) between the locations of the new point \(r_0\) and the surrounding observations \(r_i\).

$$\begin{aligned} \lambda _i = \frac{d_i^{-k}}{\sum _{i=1}^n {d_i^{-k}}}. \end{aligned}$$
(8)

The power \(k\) of the distance can be adjusted to any positive value (Webster and Oliver 2007). In this study, a value of 1 is used to model linear dependence.

The uncertainty of interpolation \(\sigma _{\mathrm{IDW}}\) can be estimated by error propagation as

$$\begin{aligned} \sigma _{\mathrm{IDW}}^2 = \frac{1}{n-1} \sum _{i=1}^n \lambda _i (Z(r_i)-Z^*(r_0))^2. \end{aligned}$$
(9)

Points that have observations retain their observed value, and their interpolation error is set to zero.

2.2.2 Ordinary Kriging

Kriging uses a variogram for the calculation of weights \(\lambda _i\). Variograms describe the inquired process in terms of the second moments of value differences in their dependence to distance \(h\). The sample variogram value \(\hat{\gamma }\) for distance \(d\) can be calculated from the observations by

$$\begin{aligned} \hat{\gamma }(d)=\frac{1}{2n}\sum _{i=1}^n(Z(r_i)-Z(r_i+d))^2. \end{aligned}$$
(10)

These values are calculated for several distance classes representing an interval of discrete width. A specific variogram model \(\gamma \) is fitted to this sample variogram \(\hat{\gamma }\). The characteristic parameters describing it are sill (representing the variance of the process), range (corresponding to the maximum distance at which correlation between the values can be observed) and nugget (a discontinuity at the origin). This discontinuity at distances near zero is caused, for example, by limitations of the sampling density (Chilès and Delfiner 2012).

In the formulation of the kriging system, this modeled function is applied to the distances \(d\) between the points involved according to

$$\begin{aligned} \begin{bmatrix} \gamma (d_{11}) &{}\quad \dots &{}\quad \gamma (d_{1n}) &{}\quad 1 \\ \vdots &{}\quad \ddots &{}\quad \vdots &{}\quad \vdots \\ \gamma (d_{n1}) &{}\quad \dots &{}\quad \gamma (d_{nn}) &{}\quad 1 \\ 1 &{}\quad \dots &{}\quad 1 &{}\quad 0 \\ \end{bmatrix} \begin{bmatrix} \lambda _{1} \\ \vdots \\ \lambda _{n} \\ m \\ \end{bmatrix} = \begin{bmatrix} \gamma (d_{10}) \\ \vdots \\ \gamma (d_{n0}) \\ 1 \\ \end{bmatrix}. \end{aligned}$$
(11)

This equation system is then solved for the weights \(\lambda _i\). \(m\) is the Lagrange parameter, which completes the system. Interpolation is then applied according to Eq. (7).

The resulting kriging variance \(\sigma _{\mathrm{OK}}^2\) at a certain point equals the minimized mean squared error on which the kriging formulation is based (Cressie 1988).

$$\begin{aligned} \sigma _{\mathrm{OK}}^2= \sum \limits _{i=1}^n \lambda _i \gamma (d_{i0}) + m. \end{aligned}$$
(12)

For OK, the value for \(\gamma (d_{ii})\), that is the variogram value for zero distance, is zero, so that the main diagonal in the matrix of Eq. (11) consists of zeros.

2.2.3 Filtered Kriging

In the case of no errors, the nugget of the variogram consists only of microscale variation. For noisy data, an error component has to be considered. The aim is to derive values of the error-free component \(T\) from measurements \(Z\) corrupted with noise \(\epsilon \) (Christensen 2011).

$$\begin{aligned} Z(r) =T(r) +\epsilon (r). \end{aligned}$$
(13)

If the variance \(\sigma ^2\) of the error \(\epsilon \) is known and assumed to be homogeneously distributed, for example due to known measurement errors, it can be introduced into kriging. Different notations can be found for example at Delhomme (1978), Cressie (1993), Rühaak (2015). Based on Christensen (2011), FK can be expressed as

$$\begin{aligned} \begin{bmatrix} \gamma (d_{11}) &{}\quad \dots &{}\quad \gamma (d_{1n}) &{}\quad 1 \\ \vdots &{}\quad \ddots &{}\quad \vdots &{}\quad \vdots \\ \gamma (d_{n1}) &{}\quad \dots &{}\quad \gamma (d_{nn}) &{}\quad 1 \\ 1 &{}\quad \dots &{}\quad 1 &{}\quad 0 \\ \end{bmatrix} \begin{bmatrix} \lambda _{1} \\ \vdots \\ \lambda _{n} \\ m \\ \end{bmatrix} = \begin{bmatrix} \gamma _{10} \\ \vdots \\ \gamma _{n0} \\ 1 \\ \end{bmatrix}, \end{aligned}$$
(14)

with

$$\begin{aligned} \gamma _{i0} = \Bigg \{ \begin{array}{cl} {\gamma }(d_{i0})-\frac{\sigma ^2}{2}, &{}\quad d_{i0} \ne 0 \\ \frac{\sigma ^2}{2}, &{}\quad d_{i0}=0 \end{array}. \end{aligned}$$
(15)

The kriging variance is

$$\begin{aligned} \sigma _{\mathrm{FK}}^2= \sum \limits _{i=1}^n \lambda _i \left( \gamma _{i0}-\frac{\sigma ^2}{2}\right) + m. \end{aligned}$$
(16)

This method leads to a filtering of the input data set, so that the values of the observed points are modified, depending on the error variances used. In contrast to OK, the kriging variance at observed points is no longer zero.

2.2.4 Heterogeneous Measurement-Error-Filtered Kriging

HFK, which incorporates heterogeneous measurement errors into kriging, was developed by Christensen (2011) and successfully applied, for example, by Christensen and Sain (2012), Christensen and Berrett (2016), Kang et al. (2017). The error variances \(\sigma _i^2\) of the observations are used individually in Eq. (18), after modifying the variogram. The values of the original variogram \({\gamma }\) are reduced by the mean of the individual error variances \(\sigma _i^2\) of the observations.

$$\begin{aligned} {\gamma ^*}(d)=\Bigg \{ \begin{array}{cl} {\gamma }(d)-\frac{1}{n}\sum _{i=1}^n \sigma _i^2, &{}\quad h \ne 0 \\ 0, &{}\quad h = 0. \end{array} \end{aligned}$$
(17)

The arithmetic mean of the RAA-derived a posteriori errors used in Eq. (17) for HFK variogram adjustment is used in this study as homogeneous error for FK in Eq. (14). It differs between the different RAA calculations.

The newly modeled variogram \(\gamma ^*\) is used in the HFK equation.

$$\begin{aligned} \begin{bmatrix} \gamma ^*(d_{11}) &{}\quad \dots &{}\quad \gamma ^*(d_{1n}) + \frac{\sigma _1^2+\sigma _n^2}{2} &{}\quad 1 \\ \vdots &{}\quad \ddots &{}\quad \vdots &{}\quad \vdots \\ \gamma ^*(d_{n1}) + \frac{\sigma _n^2+\sigma _1^2}{2} &{}\quad \dots &{}\quad \gamma ^*(d_{nn}) &{}\quad 1 \\ 1 &{}\quad \dots &{}\quad 1 &{}\quad 0 \\ \end{bmatrix} \begin{bmatrix} \lambda _{1} \\ \vdots \\ \lambda _{n} \\ m \\ \end{bmatrix} = \begin{bmatrix} \gamma ^*(d_{10}) + \frac{\sigma _1^2}{2} \\ \vdots \\ \gamma ^*(d_{n0}) + \frac{\sigma _n^2}{2}\\ 1 \\ \end{bmatrix}. \end{aligned}$$
(18)

According to Eq. (17), the variogram value for zero distances is zero. Therefore, the main diagonal of the matrix in Eq. (18) is zero, just as for OK. All other elements of the matrix incorporate the individual error variances.

The kriging variance for HFK is defined analogously to OK.

$$\begin{aligned} \sigma _{\mathrm{HFK}}^2=\sum \limits _{i=1}^{n} \lambda _i \left( \gamma ^*(d_{i0}) + \frac{\sigma _i^2}{2}\right) + m. \end{aligned}$$
(19)

This method leads to a filtering of the input data set. In contrast to FK, the filtering considers the individual uncertainties at the observation points.

3 Simulation Setup

In order to assess the different RAA models and interpolation methods, simulations are performed on synthetic data sets. Figure 1 introduces the area investigated in this study. The area of approximately 23,000 km\(^2\) covers the lower part of the NEGIS drainage system, based on a slight modification of its delineation by Zwally et al. (2012) (Fig. 1). The main glaciers are Nioghalvfjerdsbræ (79N) and Zacchariæ Isstrøm (ZAC) (cf. Fig. 1a). They have shown substantial changes in the past 10 years (Mouginot et al. 2015; Mayer et al. 2018). The CryoSat-2 measurements in this area are mainly in SARIn mode, except for the southwestern corner. The low-resolution mode (LRM) measurements were not included in this study. All geographic data sets and figures are in polar stereographic projection, with 45\(^\circ \)W as the central meridian (European Petroleum Survey Group Geodesy (EPSG) Code EPSG:3413).

Fig. 1
figure 1

a Location of synthetic data (green) used in this article and the discussed glaciers Nioghalvfjerdsbræ (79N) and Zacchariæ Isstrøm (ZAC). b Simulated height changes. c Assumed true topography (TanDEM-X DEM). d Topography-derived slope. All with 100 m resolution

Full size image

3.1 Simulation Data

Several real data sets were combined to obtain an authentic simulated data set. A rate of elevation change (Fig. 1b) is simulated by summing up contributions related to position, flow velocity, elevation and surface mass balance patterns. For each position \(i\), height change is simulated as

$$\begin{aligned} \frac{\mathrm{d}h}{\mathrm{d}t}_i=b_0 + b_1 x_i + b_2 y_i + b_3 v_i + b_4 h_i + b_5 s_i. \end{aligned}$$
(20)

The terms \(b_0 + b_1 x_i + b_2 y_i\) simulate a component with a simple linear dependence on position. This reflects height loss from southwest to northeast, based on topography and location of the outlet glaciers. The term \(b_3 v_i\) creates a trend that depends on ice flow velocity provided by Joughin et al. (2010a, b), in order to mimic changes related to ice flow dynamics. \(b_4 h_i\) denotes topography-related changes, using the TanDEM-X (TerraSAR-X add-on for Digital Elevation Measurements) DEM (Krieger et al. 2007; Rizzoli et al. 2017). The term \(b_5 s_i\) introduces an additional spatial pattern. \(s\) is taken as a temporal snapshot of the cumulative surface mass balance anomaly at a certain time from RACMO 2.3 (Regional Atmospheric Climate Model) (Noël et al. 2015). The factors \(b_i\) balance the different components and adjust the annual height change between ± 2 m year\(^{-1}\). The spatial resolution of this data set can be adjusted. It is 100 m for simulated CryoSat-2 measurements, and is adapted to the RAA resolution (250 m to 2500 m) for comparisons with RAA results later on.

Locations and times of the actual CryoSat-2 measurements between December 2010 and January 2014 were used for simulation. The measured heights were simulated based on the TanDEM-X DEM (set as initial topography of June 2009) and the synthetic height change rates. A second data set of simulated height measurements is produced by adding simulated errors. The errors are generated as pseudo-random numbers from a uniform distribution between \(-c\) and \(c\), where the slope-dependent \(c\) is defined based on a study of altimetry precision by Schröder et al. (2019) with

$$\begin{aligned} c = 0.11 \text {m} + 0.79\cdot \mathrm{slope}^2 \frac{\text {m}}{\text {degree}^2}. \end{aligned}$$
(21)

In RAA, these errors are introduced as a priori standard errors for the noisy data.

3.2 Simulation Procedure

The simulated height measurements are used as input for RAA, where different choices of RAA cell size and local topography model are assessed. Figure 2 emphasizes the role of local topography in RAA applications. For a selected cell with 2,000 m diameter, the topography of TanDEM-X DEM and three differently parametrized local topography models are compared. As the local topography models are used to reduce the satellite observations to the cell center, model and reality should optimally match.

Fig. 2
figure 2

Illustration of RAA method and local topography modeling. a Distribution of cells (black circles) with 2000 m diameter and the CryoSat-2 measurement locations (gray dots) for a selected area with underlying TanDEM-X DEM. For the red highlighted cell, the modeled local topography using different numbers of parameters (indicated in the top left corner) is shown on the right side, namely b with three parameters, c with six and d with nine parameters. The synthetic CryoSat-2 height measurements are shown in the same color scale

Full size image

RAA was applied to both the error-free and the noisy data set, with cell diameters of 500 m, 1,000 m, 2,000 m, 3,000 m, 4,000 m and 5,000 m. Because of the rather short time span and the main focus on the trend of height change, seasonal parameters were not included. The RAA cells are distributed on a regular grid (Helm et al. 2014; Schröder et al. 2019; Sørensen et al. 2018a) with overlapping areas (Wouters et al. 2015; Sørensen et al. 2018a). The local topography fit is parametrized using three, six and nine parameters (cf. Eqs. (3–5)), and DEMs are used to subtract the topography in each cell before the height change trend is estimated. The TanDEM-X DEM represents the true surface topography used for simulation, which is usually not available in real data applications. Therefore, additional DEMs were introduced, namely the Greenland Mapping Project (GIMP) DEM (Howat et al. 2014) and the ArcticDEM (Porter et al. 2018). They are likely to be used in actual applications because of their Greenland-wide coverage. The use of TanDEM-X DEM is abbreviated with T in this article, the GIMP DEM with G and the ArcticDEM with A. The parametrized local topography models are distinguished by the number of parameters, three, six and nine.

Fig. 3
figure 3

ArcticDEM (a–d) and GIMP DEM (e–h), restricted to the area of interest. a, e Surface topography, c, g corresponding slope. b, f and d, h: differences in the TanDEM-X DEM surface topography and slope, respectively. For the height subtraction in b, f, the values are mean-adjusted with a mean of 2.687 m and \(-27.143\) m, respectively

Full size image

Figure 3 depicts the differences between the additional DEMs and the TanDEM-X DEM, which is introduced as the true topography. The DEMs differ in data source and nominal time, but are similar in spatial resolution (about 100 m). The mass loss between different height acquisitions leads to different surface heights. After offset removal, the influence of fast-changing surface heights at the outlet glaciers on the different DEMs becomes apparent.

For RAA, the slope is of main interest. The slopes of the ArcticDEM match well those of the TanDEM-X DEM, except for some distinct features. Significant discrepancies occur for the GIMP DEM in some regions, which will influence the RAA results. The reasons for these differences, involving different data acquisition methods and time spans, are not a focus of this article and are therefore not further discussed.

After RAA, an outlier detection was applied on the resulting parameters. The simulated height changes vary between \(-2\) m year\(^{-1}\) and \(+2\) m year\(^{-1}\). Therefore, absolute height changes exceeding 10 m year\(^{-1}\) are removed. Additionally those with an a posteriori standard deviation of more than 1 m year\(^{-1}\) are removed. This criterion is commonly applied to RAA results in Greenland (Simonsen and Sørensen 2017) and affects less than 0.1% of the results. Prior to interpolation, a bicubic function is removed from the original data, in order to reduce the influence of spatial trends on the variogram modeling. This bicubic function was re-added after interpolation.

The gridded elevation trend \(\frac{\mathrm{d}h}{\mathrm{d}t}\) (hereafter denoted \(\hat{g_i}\)), derived from synthetic data by applying RAA and further interpolation, is compared with the original synthetic (“true”) elevation trend (\(g_i\)). The calculation of true values \(g_i\) is adapted to the spatial resolution at which the estimates \(\hat{g_i}\) are calculated. That is, for each RAA cell size a data set with true values is calculated based on the data of Eq. (20) in the respective resolution. The differences \(g_i-\hat{g_i}\) are termed true errors.

The accuracy of the results is assessed by the root-mean-square error (RMSE) over all \(n\) grid cells.

$$\begin{aligned} \text {RMSE} = \sqrt{\frac{\sum _{i=1}^{n}{(g_i-\hat{g_i})^2}}{n}}. \end{aligned}$$
(22)

In contrast, the standard uncertainties are defined as the square root of the kriging and interpolation variances.

4 Results

4.1 RAA Performance

The various cell size and local topography model combinations lead to different RAA results. Figure 4 shows the ratio of grid cells with a successful height change estimate by RAA, dependent on cell size and local topography model and restricted by outlier criteria. Cells without valid RAA results need to be filled by interpolation. The use of noisy versus error-free measurements leads to only a negligible difference with regard to the spatial coverage with RAA results. The size of the RAA cells strongly affect the coverage ratio, as larger cells cover the gaps between subsatellite tracks, while smaller cells adhere more closely to the tracks.

Fig. 4
figure 4

Ratio of cells with valid RAA results after outlier detection depending on cell size and local topography model. Cell diameter and local topography model are indicated in the lowest and second lowest rows, respectively

Full size image

The choice of the local topography model significantly affects the coverage ratio for the 500 m and 1,000 m cell sizes. For a diameter of 500 m, increasing the number of parameters from three to six and nine decreases the number of successful estimates dramatically. Here, the number of observations restricts the quality of parameter estimation. As cell size increases, the influence of the local topography model on the success of the height change estimates decreases. For cells with diameters of 3,000, 4,000, and 5,000 m, the three-parameter model of local topography yields slightly less valid RAA results than the six- and nine-parameter models. The use of DEM subtraction generally leads to more successful height change estimates, especially for the smallest cell size.

The spatial pattern of a posteriori uncertainties reveals a reason for this local topography model-dependent coverage. Figure 5 shows the standard deviations for cells with 3,000 m diameter. The spatial pattern is related to slope, leading to larger uncertainties in the steep areas, for example near the steeply sloped zone leading into the glacier tongue of 79N. The more complex the topography modeled, the better the height change estimate and the lower its a posteriori uncertainty. In this investigation, values with an a posteriori RAA-derived standard deviation of more than 1 m year\(^{-1}\) are rejected (cf. Sect. 3.2). The rejections lead to gaps in the steep areas, depending on the local topography model.

The true errors, \(\hat{g_i}-g_i\), for the selection of Fig. 5 are shown in Fig. 6. Comparison with Figs. 1d and 5 indicates that, similar to the formal uncertainties, the actual errors depend on slope. Striking is the case where the GIMP DEM is used. There, the actual errors are comparably high, even in flat areas. This observation concerning the GIMP DEM subtraction is similar for all cell sizes. The use of TanDEM-X DEM or ArcticDEM leads to both low uncertainties and small actual errors, except for the steepest areas. For the cells with a diameter of 1,000 m and less, no apparent slope-dependent pattern of true error occurs (not shown). These observations apply for both the error-free and the noisy data set.

The standard uncertainties from RAA estimation are compared with the true errors. Figure 7 illustrates a nearly linear relationship between them for RAA cells with at least 3,000 m diameter, which confirms the visual assessment by Figs. 5 and 6. The standard uncertainties slightly overestimate the true errors.

Fig. 5
figure 5

Standard uncertainty of RAA-estimated height changes, using noisy data and cells with 3,000 m diameter. The local topography model is indicated in the top left corner

Full size image
Fig. 6
figure 6

True error, \(\hat{g}_i-g_i\), of RAA-estimated height changes from noisy data for cells with 3,000 m diameter. The local topography model is indicated in the top left corner

Full size image
Fig. 7
figure 7

Relationship between true error (ordinate) and RAA standard uncertainty (abscissa) for noisy data. For clarity, RMS values were calculated over 20 intervals between 0 and 1 m year\(^{-1}\) of sigma. The RAA cell size is indicated in the top left corner; the local topography models are color-coded

Full size image

For the cells with diameters of 500 m to 2,000 m, small standard uncertainties are significantly exceeded by true errors. This is similar for all local topography models, and most pronounced for the GIMP DEM. The difference of this DEM from the simulated topography has a very high effect when small RAA cells are used.

Figure 8 represents the color-coded RMSE information for the different combinations of local topography models, cell sizes and the noisy and error-free data sets. The RMSE spans from 0.13 to 0.85 m year\(^{-1}\). This is a difference of factor 6.5. The RMSE of the thee-parameter local topography model increases with increasing cell size. The six- and nine-parameter models have the lowest RMSE for cell diameters of 2,000 m and 3,000 m. In contrast, the use of a DEM leads to the best RMSE with large cell diameters. It is striking that the GIMP DEM-reduced RAA shows the worst results among all local topography models.

Fig. 8
figure 8

RMSE for the RAA-derived height change estimates for error-free (left) and noisy (right) data. Cell diameters are indicated on the left, topography models at the bottom

Full size image

The pseudo-random errors of the noisy data set affect the results and lead to higher RMSE values compared with the error-free data set. The effect of the added noise is strongest for cells with diameters of 2,000 m and less. Asides from this effect, the noisy and error-free data sets show similar results. The best results are obtained with 5,000 m cell diameter and the TanDEM-X DEM.

In conclusion, different combinations of cell size and local topography models can lead to satisfactory height change estimates. In the applied specific constellations, cells with diameters of less than 2,000 m do not cope well with perturbations, such as the additional random errors. Larger cells include more observations in the estimation process and therefore are more capable of determining the desired linear height change estimates. If no DEM is used, high parametrized local topography models should be preferred over a plane-fit model, as long as a sufficient number of observations are available.

The DEM subtraction leads to good results, as long as the DEM is close to the actual topography. As shown in Fig. 3, the ArcticDEM is close to the assumed true topography of the TanDEM-X DEM. Therefore it leads to better results than the GIMP DEM. The importance of finding suitable DEMs concerning resolution and matching time span has already been addressed by Sørensen et al. (2011) and is well demonstrated here.

As the difference between noisy and error-free height change estimates is not significant, and the assumption of noise is assumed to better reflect the actual situation, further analysis uses the noisy results only.

4.2 Impact of Variogram Models on Kriging Interpolation

Before interpolation is applied, the influence of different variograms is analyzed. This analysis is done based on an RAA data set with 2,000 m cell diameter and nine topography parameters. HFK with different variogram selections is applied. The selected RAA result leaves sufficient area for interpolation, so that the effect of different variograms on the kriging result can be properly studied.

Fitting a variogram is an essential part of kriging, as different spatial distances, class divisions, weighting schemes and variogram models have to be considered. The sample variograms were calculated in 30 distance classes (cf. Sect. 2.2.2) ranging from zero to the chosen maximum distance. To fit the models, the different classes were weighted with \(p\) depending on sample distance \(h\) and number of observations per class \(n\) (Pardo-Igúzquiza 1999) as

$$\begin{aligned} p=\frac{n}{h^2}. \end{aligned}$$
(23)

Variograms with three options for the maximum distance (10; 50; 100 km) and two options for the analytical model (Gaussian and spherical) are considered, which results in a total of six variogram model options. The choice of parameters for variogram modeling depends on the assumptions on the underlying physical processes. Changes in ice heights proceed on small and large scales. However, due to the radar footprint and the coverage with satellite data, as well as the restrictions of RAA, the actual spatial and temporal resolution is limited. The sample variogram can provide different solutions depending on the scales considered. In the process of interpolation, points are selected out of eight sectors, with a maximum of three points per sector. Therefore, correct modeling of the variogram on short distances (depending on the cell size) is of great importance.

The resulting variograms are illustrated in Fig. 9. The shapes of the sample variograms differ slightly according to the spatial scales. The sharp increase at distances up to 3 km (Fig. 9a) is less pronounced for large maximum distances (Fig. 9c). Further increase of variogram values at more than 70 km can be neglected, as it exceeds the maximum distance between observation points for interpolation. The largest data gap in the RAA results (at the heavily sloped region inland from the grounding line of ZAC) spans an area of approximately 15 km times 40 km.

Fig. 9
figure 9

Variograms of the noisy results with 2,000 m cell diameter and nine topography parameters, calculated for a maximum difference of a 10 km, b 50 km and c 100 km (note the different distance scales). Sample (black) and fitted Gaussian (blue) and spheric (green) variograms are shown

Full size image

For further investigation, HFK was applied to the selected RAA data set using the six differently modeled variograms. Here, the focus is on the differences between the results induced by different variogram models. These analyses show that the effect of the choice of variograms on the final HFK result is negligible. The RMSE values differ marginally (e.g. 0.129 m year\(^{-1}\) to 0.145 m year\(^{-1}\) for the complete interpolated result).

In contrast, the kriging uncertainty is significantly affected by the choice of the variogram, which thus affects the realism of the uncertainty characterization for the interpolation result. The kriging uncertainty should reflect, in a statistical sense, the true error. In Fig. 10 the kriging uncertainty is plotted against statistics of the true error. Additionally, the underlying number of points per ratio is illustrated as relative density. With a maximum distance of 10 km, a nearly linear relationship is obtained. As an exception, for the lowest uncertainty bin, a strong discrepancy is observed between uncertainty and error. This is caused by just five cells, where the RAA standard uncertainty, that is, the uncertainty of the input to HFK, is significantly underestimated.

Fig. 10
figure 10

HFK standard uncertainty (abscissa) is plotted against the true error (ordinate) for different variogram models, applied to the noisy RAA results with 2,000 m diameter and nine topography parameters. For clarity, RMS values (black dots) were calculated over 20 intervals between 0 and 0.3 m year\(^{-1}\) of sigma. The number of individual values are shown color-coded from many (purple) to fewer (yellow). The specific variogram model and maximum distance are indicated in the top left corner

Full size image

The meaningfulness of the kriging uncertainty is best achieved for the fit with a maximum distance of 10 km. The performance of the spheric and Gaussian variogram models is comparable. This investigation shows that the variogram affects mainly the kriging uncertainty, and less the interpolation result. Therefore, the variogram modeling for interpolation of height changes should be focused on short distances. For the following kriging, a spherical variogram model is fitted to sample variograms spanning a maximum distance of 10 km.

4.3 Interpolation Performance

Interpolation was applied to the noisy data sets of the RAA results using the four interpolation methods IDW, OK, FK and HFK. The process included the calculation of the sample variogram, the fitting of the model variogram, the calculation of weights and finally the interpolation itself.

The RMSE values for the interpolated height changes are illustrated in Fig. 11. In contrast to IDW and OK, FK and HFK change the values of grid cells that have valid observations. Therefore, grid cells that have no valid value before interpolation are marked as “interpolated” and the entirety of the grid cells is marked as “complete”, while the statistics of the filtered values are marked as “filtered”.

Fig. 11
figure 11

RMSE after interpolation with IDW (first row), OK (second row), FK (third row) and HFK (last row), distinguished between interpolation (middle) and complete result (right), as well as filtered (left) for FK and HFK. Cell diameters are indicated on the left, topography models at the bottom. The initial RAA result (copy from Fig. 8b) is additionally shown in the top left, to enable comparison

Full size image

For large cell sizes (4,000 m; 5,000 m diameter) the RMSE of the complete grid is determined predominantly by the RMSE of the filtered values, rather than by the RMSE of the interpolated values, because of the high coverage of RAA results (cf. Fig. 4).

For IDW, the RMSE of interpolation, shown in Fig. 11b, increases with increasing cell size for the parametrized local topography models. The use of DEMs does not lead to a clear advantage of certain cell sizes in interpolation. The complete result (Fig. 11c) shows increasing RMSE with cell size for the parametrized models and decreasing RMSE with cell size for the DEM subtraction. Compared with the underlying RAA result (see Fig. 11a), the RMSE values are comparable or slightly lower. Notably, for small cell sizes, where many grid cells are interpolated, the RMSE of interpolated values is smaller than that for the original RAA results. This indicates that the weighted averaging process of interpolation reduces noise in the RAA results. The best complete results are achieved with 5,000 m cell diameter and TanDEM-X DEM.

Similar to IDW, the interpolation by OK (Fig. 11d) performs best for small cells and parametrized models or with the TanDEM-X DEM or ArcticDEM and 4,000 m or 5,000 m cell size. But the RMSE values are higher than for IDW, especially when DEMs are used.

FK (Fig. 11f–h) interpolates similarly to OK, but considers a constant error that is used to filter the observation points. This leads to improved RMSE values at these points, as well as an improved complete result compared with OK and IDW. The pattern of RAA RMSE is thereby maintained, leading to the best complete results with TanDEM-X DEM and 4,000 m, 5,000 m cell diameters.

For HFK, the filtering of the RAA results improves the RMSE significantly, much more than for FK (cf. Fig. 11i, with Fig. 11a, f). In part, the RMSE decreases by more than 0.3 m year\(^{-1}\) after filtering with HFK (e.g. 500 m cell diameter and GIMP DEM subtraction). The height change rates best reflecting the simulated truth are achieved with cell diameters of 2,000 m to 4,000 m and a local topography correction with six or nine topography parameters or the TanDEM-X DEM or the ArcticDEM.

Table 1 shows the RMSE values for the observation points, interpolation and complete results for a chosen example: the noisy data set with cell diameter of 3,000 m and nine topography parameters. The quality of IDW and OK interpolation is similar, with IDW performing slightly better. FK interpolation performance is similar to OK, but the filtering included at observation points improves the complete result. HFK not only filters better than FK, but in many cases is even able to perform better interpolation, leading to significantly improved complete height change results. The accuracy improvement for this example (3,000 m cell diameter and nine topography parameters) is 72% between OK and HFK.

Table 1 RMSE [m year\(^{-1}\)] for the noisy data set with 3,000 m cell diameter and nine-parameter topography model
Full size table

An example of the spatial pattern of the interpolated height changes, differences and standard uncertainties can be seen in Fig. 12. The results of IDW and OK are very similar and show speckled patterns. Their standard uncertainties neglect uncertainties at observation points (value zero) and increase with distance to them. In the interpolated areas, the OK standard uncertainties are generally higher than for IDW, and with less variation. The standard uncertainties are further discussed in Sect. 4.4. FK application leads to a smoother height change result and less error compared with OK. The pattern of uncertainties is similar, but has a fixed value (not zero) at observation points. HFK leads to less error and a much smoother height change result than the other interpolation methods, which is due to the spatially varying filtering. The spatial pattern of the standard uncertainties does not simply reflect the existence of observations, but gives more reliable information about areas with higher uncertainties, which are mainly the sloped regions near the grounding lines of the two glaciers. Additionally, the uncertainties are no longer zero at observation points.

Fig. 12
figure 12

Resulting height changes (left), true error (center) and standard uncertainties (right) for the noisy data set with 2,000 m cell diameter and nine-parameter local topography model. The interpolation methods are indicated in the top left corner of the left-hand plots

Full size image

The southwestern corner of the study area is not observed by CryoSat-2 in SARIn mode, but only in LRM mode. The consequent data gap is filled via extrapolation by the different interpolation techniques. OK is not recommended for extrapolation because the extrapolated values approach the data mean (Chilès and Delfiner 2012). This is reflected in the corresponding higher standard uncertainties. IDW extrapolates slightly better than OK and FK. Although HFK is based on OK, it performs best in extrapolation as well.

The comparison of different interpolation methods shows that HFK is best suited for application to height changes derived from satellite altimetry. In particular, the filtering improves the results substantially. More simple geostatistical methods such as OK do not necessarily outperform other approaches such as IDW.

4.4 Uncertainties

The kriging standard uncertainties of the four interpolation methods are investigated to obtain more information about their reliability. As in Sects. 4.1 and 4.2, the relationship between true error and the standard uncertainties is analyzed. In Fig. 13, this relationship is shown for all cell sizes and interpolation methods. The focus here is on the distinction between the different interpolation approaches. Therefore, the different topography models are not plotted separately.

Fig. 13
figure 13

Relationship between standard uncertainty and true error for the noisy data of the complete result. For clarity, the RMS values were calculated over 100 intervals between 0 and 1 m year\(^{-1}\) of sigma as a mean of all topography models. The RAA cell sizes are shown in meters in the top left corner

Full size image

IDW has a rather linear relationship, except for the observed points, where the RAA value is maintained and the uncertainty set to zero. The standard uncertainties underestimate the errors.

The inconsistency between error and uncertainties of RAA results observed in Fig. 7 for cells with a diameter of 2,000 m and less propagate to inconsistencies for OK, FK and HFK uncertainty estimates. This can be seen in Fig. 13a–c, where no simple linear relationship between uncertainties and errors is visible. For cell diameters of 2,000 m and more, the plotted OK and FK relations are more scattered than those for HFK and IDW. While the uncertainties of OK and FK for cell diameters of 2,000 m and below overestimate the errors, the opposite happens for cell diameters of 3,000 m and larger. Similar to IDW, the observation points of OK are provided with uncertainty value of zero. FK and OK show very similar behavior in the relationship of error and uncertainty analyzed here.

The HFK uncertainties estimated with data based on cell diameter of at least 3,000 m represent the actual errors very well up to approximately 0.25 m year\(^{-1}\). Higher uncertainties underestimate the errors. The best accordance of errors and estimated uncertainties is achieved with HFK for 3,000 m and 4,000 m cell diameter.

This investigation shows that standard uncertainties should be handled with care. The values at observation points for OK and IDW in particular are not meaningful. HFK improves the reliability of the uncertainty estimate.

5 Conclusions

To obtain reliable information about the performance of different RAA configurations and interpolation methods, a synthetic data set was created in order to compare the derived height changes with known true values.

It was shown that the RAA results differ depending on the cell size and topography parametrization. In these investigations, the models with six and nine topography parameters lead to good results. The smallest analyzed cell size of 500 m does not cope well with the induced random errors, while cells with a diameter larger than 4,000 m can lead to larger errors than those with smaller cells. The best results using parametrized models are achieved with a cell diameter of 3,000 m and nine topography parameters. The results of DEM subtraction in RAA depend very much on the quality of the DEM; the more the DEM represents the actual sampled topography, the better the results. ArcticDEM and TanDEM-X DEM (assuming true topography), in conjunction with cells of 3,000 m and greater diameter, provide the most accurate RAA results among all options tested. As the agreement of the DEM with the topography is difficult to assess when real data is used, DEM subtraction should be applied carefully. The a posteriori standard errors of RAA are reliable for cells with at least 3,000 m diameter, and can be used for filtering with FK and HFK.

The variograms used for kriging focus on short distances, as this gives the best results for interpolation and reliable uncertainties. The subsequent interpolation was accomplished with IDW, OK, FK and HFK. OK and IDW performed comparably well. The resulting height changes are improved by the filtering included in the FK and HFK algorithm. The best results are achieved by incorporating heterogeneous errors with HFK. Additionally, the corresponding standard uncertainties are reliable, and their spatial patterns reflect actual errors.

In this study, linear height change is the parameter of interest. This is a simplification of the real process, as interannual signals are present. They can be resolved by RAA and included in spatiotemporal interpolation.

Further research based on other regions and satellite missions, especially pulse-limited radar data, could expand the applicability of these results. Additionally, the influence of outlier criteria and selection of points used for interpolation can be elaborated. Stacked variograms are another approach to cope with the different behavior apparent at different spatial scales that would be worth considering.

Based on these investigations, HFK can be recommended to achieve full spatial coverage of height changes from satellite altimetry measurements derived by RAA, as the results show the smallest error and least speckle, and provide meaningful and reliable uncertainties.