Data and pre-processing

The surface elevation of the ice-sheet periphery up to 20–50 km inland is well constrained due to about 3500 aerial photographs collected by the Danish Agency for Data Supply and Efficiency between 1978 and 1987 that Korsgaard and others (Reference Korsgaard2016) used to produce DEMs with a grid resolution of 25 m. In the following, we collectively refer to these as the AeroDEM. The elevation of the interior surface, on the other hand, is less constrained. Seasat and GEOSAT DEMs (Zwally and Brenner, Reference Zwally and Brenner1997) use observations from 1978 and 1985–86 but have a low resolution of 20 and 10 km, respectively. Furthermore, they are geographically restricted to the southern half of the ice sheet, up to $72^{\circ}\,\mathrm{N}$. Thus, we instead constrain the interior elevations through laser altimetry data from 1993 and 1994 that NASA collected using the airborne topographic mapper (ATM; Krabill and Thomas, Reference Krabill and Thomas2010). These altimetry flight lines sparsely cover the entire ice-sheet interior with a maximum distance of around 200 km between data points (Fig. 1). We ignore the time difference between the AeroDEM and ATM acquisitions and assume that they collectively represent the state of the ice sheet around 1980. This assumption is addressed further in the Discussion section.

Figure 1.

The field of standardized elevation differences $\Delta h_\mathrm{z}$ of the GrIMP DEM version 2 (Howat and others, Reference Howat, Negrete and Smith2022) minus historic elevations (either AeroDEM (Korsgaard and others, Reference Korsgaard2016) or ATM (Krabill and Thomas, Reference Krabill and Thomas2010)), which are the input data for the Gaussian process regression. See the Gaussian process regression Methods section and Eqn (2) for details on the standardization. The dense raster data at the margins correspond to the AeroDEM, while the flight line point data in the interior come from the ATM. Note that positive values do not necessarily mean that GrIMP elevations are higher than historic elevations due to the detrending. The shown projection is NSIDC Sea Ice Polar Stereographic North (EPSG 3413).

Figure 1 Long description

Standardized observations delta h sub z mapped over the outline of Greenland. Map elements and orientation: - Title text: Standardized observations delta h sub z. - Horizontal axis label: Easting (km), with tick labels -600, -300, 0, 300, 600. - Vertical axis label: Northing (km), with tick labels -1000, -1500, -2000, -2500, -3000. - Northing increases upward and easting increases to the right. Legend: - A vertical legend bar is shown at the right. - Tick labels on the legend: -3, -2, -1, 0, 1, 2, 3. - The legend indicates that mapped observations take values from -3 to 3. Mapped observations and symbols: - The interior of Greenland contains many thin, straight and gently curving track-like lines that cross and intersect, forming a network. Several tracks include loops and bends. - Along parts of the coastline, there are narrow, continuous bands of mapped values that follow the coastal outline. - No place names, region names, or landmark labels are visible on the map. Spatial distribution described without place names: - Track-like lines span much of the interior from the northern portion to the southern portion. - Coastal bands appear intermittently around the perimeter, with multiple segments visible on the western, eastern and southern coasts. - Values vary along both the interior tracks and the coastal bands across the full legend range from -3 to 3.

We only consider observations for grid cells that 1) are within the ice sheet according to the mask used by the IMBIE Team (2020), manually modified using the AeroDEM to adjust to the larger extent of some outlet glaciers in the 1970s and 1980s, and 2) are not classified as bedrock according to the GrIMP Land Ice and Ocean classification Mask, Version 1 (Howat, Reference Howat2017). For the AeroDEM, we filter out values where its reliability mask provided by Korsgaard and others (Reference Korsgaard2016) has a value below 40, as recommended. Then we co-register it to the GrIMP DEM version 2 on a grid resolution of 150 m using the Nuth and Kääb (Reference Nuth and Kääb2011) algorithm, before downsampling it to the target grid size of 600 m with bilinear interpolation. We apply the same co-registration to the ATM point data.

For the ATM, we use the dataset Pre-IceBridge ATM L2 Icessn Elevation, Slope, and Roughness, Version 1 (Krabill and Thomas, Reference Krabill and Thomas2010), which is an already resampled and smoothed product of the original measurements. However, these data points are still heavily over-sampled in direction of the flight line with less than 100 m spacing between most points. In order to reduce the number of data points and noise further, we apply a median block reduction to the co-registered data using the Python package verde (Uieda, Reference Uieda2018) that averages the data over 400 m, a length scale slightly smaller than the grid size of the resulting DEM. ATM points that are within a 600 m radius of any observation in the AeroDEM are removed to avoid overlaps and to prioritize the AeroDEM. For Gaussian process regression, these data points are used directly, whereas the SVD method requires them to be projected back onto the target grid with 600 m horizontal resolution. In the following, we use the names of the different datasets to refer to the downsampled versions.

Gaussian process regression (GP)

In Gaussian process regression (from now on referred to as GP), each interpolated point is a weighted average of the input data, where the weights depend on the distance to each data point and the chosen kernel function describing the spatial covariance of the data. The underlying assumption of this method is that the input data are second-order stationary, meaning that the mean and the variance are constant over the entire domain, and that the spatial covariance only depends on the distance between points (Webster and Oliver, Reference Webster and Oliver2007). Because of the large variability in absolute elevation values, and because increments of a variable are often closer to second-order stationary than the variable itself (e.g. Chiles and Delfiner, Reference Chiles and Delfiner2012), we perform the interpolation based on elevation differences $\Delta h$ of the observed elevation data to a reference DEM:

(1)\begin{equation} \Delta h = h_\text{ref} - h_\text{source}, \end{equation}

where $h_\text{ref}$ are the elevations of an arbitrary reference DEM, for which we choose the MEaSUREs Greenland Ice Mapping Project (GrIMP) Digital Elevation Model from GeoEye and WorldView Imagery, version 2 (Howat and others, Reference Howat, Negrete and Smith2022), and $h_\text{source}$ are elevation data of either the AeroDEM or the ATM flight lines. The fact that both elevation datasets are co-registered to the reference DEM reduces extreme elevation differences in steep terrain due to slight misalignment. We choose the GrIMP DEM version 2 as a reference DEM because it is derived using more recent data and processing techniques compared to the GIMP DEM version 1, resulting in less artifacts in the calculated elevation differences.

Here, we do not assume that second-order stationarity already holds in the raw elevation differences. Instead, we standardize the data such that both the bias and variance are approximately constant in space (e.g. as in Hugonnet and others, Reference Hugonnet2022). We assume that after this step, the dataset can be represented by a constant mean function of zero, a constant variance of 1 and a representative spatial covariance function. In Supplementary Material, we qualitatively evaluate the latter by showing empirical variograms calculated separately for different subsets of the data (Supplementary Fig. S3). We derive standardized elevation differences $\Delta h_\mathrm{z}$ by removing the median $\overline{\Delta h}$ and the standard deviation $\sigma_{\Delta h}$ that both depend on the elevation $h_\text{ref}$ and the slope $s_\text{ref}$ (Horn, Reference Horn1981) derived from the reference DEM GrIMP version 2 at 600 m grid resolution:

(2)\begin{equation} \Delta h_\mathrm{z}(h_\text{ref},s_\text{ref}) = \frac{\Delta h - \overline{\Delta h} (h_\text{ref},s_\text{ref})}{\sigma_{\Delta h}(h_\text{ref},s_\text{ref})}. \end{equation}

The functions $\overline{\Delta h}(h_\text{ref},s_\text{ref})$ and $\sigma_{\Delta h}(h_\text{ref},s_\text{ref})$ are obtained by binning of $\Delta h$ values with respect to elevation $h_\text{ref}$ and slope $s_\text{ref}$, calculating the median and standard deviation of each bin (Fig. 2a–c), and applying a linear interpolation and extrapolation to both the median and standard deviation fields (Supplementary Fig. S2). The dependence on elevation captures the reduced variability in $\Delta h$ of the ice-sheet interior compared to the margin. In addition, we calculate the variability as a function of slope to account for the large variability that comes from steep bedrock slopes surrounding the outlet glaciers, where cells may overlap with the bedrock and where the observational uncertainty is typically high. After standardization, we filter outliers by removing standardized $\Delta h_\mathrm{z}$ points whose absolute values are larger than ten times the normalized median absolute deviation (NMAD) of $\Delta h_\mathrm{z}$. Finally, for convenience of interpretation, we force the standard deviation of the whole dataset to be one by dividing by its actual standard deviation, and the mean to be zero by subtracting the actual mean. In the ideal case, the resulting dataset would be normally distributed around a mean of zero. Here, there is a slight excess of values close to zero (Fig. 2d).

Figure 2.

(a) Discrete, two-dimensional standard deviation $\sigma_{\Delta h}(s_\mathrm{ref},h_\mathrm{ref})$ and (b) median $\overline{\Delta h}(s_\mathrm{ref},h_\mathrm{ref})$ evaluated over bins, used for standardization of observations according to Eqn (2); (c) number of samples in each bin; (d) distribution of standardized observations, compared to a normal distribution; and (e) empirical variogram of standardized observations and fitted theoretical variogram (Eqn (3)).

Figure 2 Long description

Five graphs labeled a to e. Graph a shows elevation difference across slope and elevation bins, with higher differences at steeper slopes and lower elevations. Graph b displays standard deviation values, indicating increased variability at higher slopes and elevations. Graph c presents the number of samples per bin, with more samples at lower slopes and elevations. Graph d is a frequency distribution plot comparing standardized observations to a normal distribution, centered near zero with a slight skew. Graph e is an empirical variogram plot showing semivariance against distance, rising initially and leveling off, with empirical points and a fitted line. The graphs collectively illustrate relationships between elevation, slope, variability and sample distribution, highlighting trends and patterns across the dataset.

In spatial statistics, the kernel function is typically described through a variogram $\gamma$, which is closely related to the spatial covariance $C$ through $\gamma(d) = \sigma^2 - C(d)$, where $\sigma^2$ is the variance and $d$ the distance between points. Here, we evaluate an empirical variogram as defined by Cressie (Reference Cressie1993) from the standardized input data $\Delta h_\mathrm{z}$, which we then use to fit a theoretical variogram function. This step replaces the optimization of hyperparameters that is commonly performed in a GP workflow using the input data directly. For the theoretical variogram, we use the sum of two functions corresponding to two different ranges, which improves the representation of the structure at different length scales (e.g. Rolstad and others, Reference Rolstad, Haug and Denby2009). Specifically, we use the sum of two exponential variograms,

(3)\begin{equation} \begin{aligned} \gamma(d) = &a_1 \left[ 1 - \exp\left(-3 \left(\frac{d}{r_1}\right) \right) \right]\\ &+ a_2 \left[ 1 - \exp\left(-3 \left(\frac{d}{r_2}\right) \right) \right], \end{aligned} \end{equation}

where $a_1$ and $a_2$ are the two sill values and $r_1$ and $r_2$ are the two correlation length scales. We scale the function $\gamma (d)$ such that the sum of the two sills equals the variance of the data, that is, $\lim_{d \to \infty} \gamma(d) = \sigma^2 = 1$. We also constrain the nugget to small values as it can be erroneously high when fitting the variogram function freely. In practice, the nugget turns out to be essentially zero, meaning that $\lim_{d\to0^+} \gamma(d) = 0$, that is, $\gamma(d)$ starts at zero for small distances $d$. A nugget of zero is coherent for a continuous variable such as elevation since we do not expect any variance at near-zero range other than observational error, which we treat separately as is common in Gaussian process regression.

To estimate the empirical variogram, we only use a random 20% of the AeroDEM data points since they have a much larger data density than the ATM data and would otherwise dominate the variogram and distort the theoretical fit at large ranges (see Supplementary Material for more details). We achieve the best fit for parameters $a_1\approx 0.2$, $r_1\approx 8~\mathrm{km}$, $a_2 \approx 0.8$ and $r_2\approx 140~\mathrm{km}$. At 50 km, roughly the radius of most data gaps, the signal is still correlated to around 25% ( $\gamma \approx 0.75$, Fig. 2e). The fact that the signal is correlated over a range larger than the size of the data gaps is also visible in the map of the standardized input data (Fig. 1), where the signal has a similar value across multiple ATM flight lines. This is an important fact that makes this type of interpolation possible.

Because we standardize the input data, we also need to account for heteroscedasticity, that is, the spatial variability of the error. To this end, we assume measurement uncertainties of 10 m for the AeroDEM based on the NMAD of the difference between AeroDEM and ATM elevations over stable terrain, and conservative 0.5 m for the ATM (cf. Krabill and others, Reference Krabill, Thomas, Martin, Swift and Frederick1995b), which are then standardized in the same way as the input data (Eqn (2)). We choose these uncertainties conservatively since the observations correspond to different acquisition periods.

These estimates correspond to the spatially independent part of the error, that is, the diagonal entries of the noise covariance matrix. Here, we also evaluate the spatial covariance of the error from the differences between AeroDEM and ATM elevations over stable terrain (cf. Hugonnet and others, Reference Hugonnet2022) and find that the error is correlated over roughly 5 km, which corresponds to around eight grid cells. We account for this spatial covariance by fitting another variogram function to populate nondiagonal entries in the noise covariance matrix (see Supplementary Material).

We use the Julia package GeoStats.jl (Hoffimann, Reference Hoffimann2018) to estimate the variogram and AbstractGPS.jl (Widmann and others, Reference Widmann2025) for the GP, and we only predict points in the ice-sheet interior, that is, excluding grid points occupied by the AeroDEM. To avoid the large memory requirements of constructing the full covariance matrix, we split the domain into blocks and make the predictions for each block at a time, taking into account the observations in bands of 140 km width on each side, a distance at which the signal is essentially decorrelated. In a last step, we de-standardize the predicted elevation differences and subtract them from the elevations of the reference DEM.

SVD method

The goal of the SVD method is to extract the spatial characteristics from a set of physically plausible DEMs and to recombine these in a way that is consistent with the AeroDEM and ATM observations. The way we extract the spatial characteristics is equivalent to characterizing images of human faces through so-called ‘eigenfaces’ (Sirovich and Kirby, Reference Sirovich and Kirby1987), which are the eigenvectors obtained from a principal component analysis and are equivalent to the singular vectors of the corresponding SVD. In this framework, each image of a face can be represented as a weighted sum of these characteristic ‘eigenfaces’, where the weights characterize the individual face well enough so that they can be used for face recognition (Turk and Pentland, Reference Turk and Pentland1991). Similarly to the gray values in an image of a human face, the elevations in a DEM of the Greenland ice sheet follow certain characteristics. We adopt the working assumption that these eigenvectors form a basis for the set of physically plausible surface elevations, and we aim to find the corresponding weights that are consistent with our sparse observational dataset. We refer to this approach as the SVD method, and to the respective reconstruction as the SVD reconstruction.

The set of physically plausible DEMs needs to represent a wide range of ice-sheet configurations that goes beyond what is presently observed, which is why we rely on model-generated elevation fields. For this purpose, we use the Parallel Ice Sheet Model (PISM; Khrulev and others, Reference Khrulev2025) to produce an ensemble of historical simulations at a horizontal grid resolution of 600 m. This grid resolution is the motivation for generating a new set of simulations rather than extracting the DEMs, for instance, from the ISMIP6 multi-model ensemble (Goelzer and others, Reference Goelzer2020) that is posted on a 1 km horizontal grid.

Our ensemble is generated by uniformly sampling from five parameters that control atmosphere and surface processes: a positive degree-day ice melt factor (0.5–12.5 mm K $^{-1}$ d $^{-1}$), a positive degree-day snow melt factor (0.5–5.5 mm K $^{-1}$ day $^{-1}$), a standard deviation of daily temperature (1–6 K), an additive scalar temperature offset ( $-$5 to +5 K) and a percentage increase in precipitation per 1 K increase in air temperature ( $-$3% to 6%). We note that our goal is to create a diverse set of DEMs rather than hindcasts that prioritize agreement with observations, and thus we sample from a parameter space that is larger than commonly used. Each simulation starts from the same initial state, with the subglacial topography and ice thickness of BedMachine version 3 (Morlighem and others, Reference Morlighem2017b) and the surface elevation of the GIMP DEM version 1 (Howat and others, Reference Howat, Negrete and Smith2015). The simulations are evolved 40 years forward in time, forced by climatic conditions of the years 1980–2020 provided by the regional climate model RACMO (Noël and others, Reference Noël2018).

Finally, we extract each of the 40 annual ice surface elevation fields of each simulation as a separate member of that ensemble, that is, without considering their temporal relation. We collect these modeled geometries in a matrix $X \in \mathbb{R}^{n\times m}$ containing the elevations $h$ for $n$ grid cells (rows) and $m$ different model geometries (columns), each representing the modeling output at a different time step or from a different simulation:

(4)\begin{equation} X = \begin{bmatrix} h_{1\,(M_1)} & h_{1\,(M_2)} & ... &h_{1\,(M_m)}\\ h_{2\,(M_1)} &h_{2\,(M_2)} & ... & h_{2\,(M_m)}\\ ... & ... & ... & ... \\ h_{n\,(M_1)} & h_{n\,(M_2)} & ... &h_{n\,(M_m)} \end{bmatrix}, \end{equation}

where $M_i$ is the $i$th dataset generated by the model. For each pixel, we subtract the mean elevation across all modeled geometries to arrive at centered values $\widehat{h}_{j\,(M_i)}$ comprising the matrix $\widehat{X}$:

(5)\begin{equation} \widehat{h}_{j\,(M_i)} = h_{j\,(M_i)} - \frac{1}{m} \sum_{i=1}^m h_{j\,(M_i)}. \end{equation}

The matrix $\hat{X}$ is then decomposed using the SVD, truncated to $r$ singular values, where $r \lt m$:

(6)\begin{equation} \hat{X} \approx U \Sigma V^T. \end{equation}

The columns of $U \in \mathbb{R}^{n\times r}$ are the left singular vectors, which are equivalent to the eigenfaces in concept, but here we refer to them as the spatial modes. They are vectors of length $n$, that is, they assign a value to each grid cell. Each of the $r$ modes represents a certain spatial characteristic prevailing across the modeled geometries (Fig. 3). The matrix $\Sigma \in \mathbb{R}^{r\times r}$ is a diagonal matrix containing the corresponding singular values, and the columns of $V \in \mathbb{R}^{m\times r}$ are the right singular vectors that, for each of the $m$ modeling datasets, assign coefficients to the $r$ modes.

Figure 3.

Modes from the SVD (Eqn (6)), where $u_i$ represents the $i$th mode, with decreasing associated singular value. We omit showing absolute values as the information content of the spatial modes lies mostly in their spatial distribution.

As depicted in Fig. 3, the first mode mostly represents the contrast between the margins and the interior of the ice sheet, as well as between the western and eastern margins. In the second and third modes, more local effects and the signatures of the largest outlet glaciers become apparent. Later modes become increasingly focused on small-scale structures, as visible in the $100$th mode.

Here we make the assumption that a set of new coefficients comprising a vector $\mathbf{v}_\mathrm{rec}$ can be used to generate a different but also physically possible state of ice-sheet elevations, which will be the reconstruction $\mathbf{h}_\mathrm{rec}$:

(7)\begin{equation} \mathbf{h}_\mathrm{rec} = U \Sigma \mathbf{v}_\mathrm{rec}. \end{equation}

We choose the coefficients $\mathbf{v}_\mathrm{rec}$ such that the corresponding reconstructed elevations $\mathbf{h}_\mathrm{rec}$ give the best fit to the AeroDEM and ATM observations, $\mathbf{h}_\text{obs}$. Specifically, we formulate the problem as regularized least squares, as it allows us to obtain an analytical solution. Therefore, we find $\mathbf{v}_\mathrm{rec}$ through

(8)\begin{equation} \min_{\mathbf{v_\mathrm{rec}}} \quad ||\mathbf{h}_\mathrm{obs}- A \mathbf{v}_\mathrm{rec}||_2^2 + \lambda ||\mathbf{v}_\mathrm{rec}||_2^2, \end{equation}

where the matrix $A = U_\text{obs} \Sigma $ comprises the rows of $U \Sigma$ corresponding to grid cells where observations are available, and $\lambda$ is a regularization parameter that penalizes large values of the coefficients $\mathbf{v}_\mathrm{rec}$. Note that we weigh every residual equally here as we have found that using weights—for instance based on different uncertainties for the two types of observations—does not improve the result in terms of the cross-validation error (see Supplementary Material). A detailed derivation of the solution of Eqn (8) is given in Supplementary Material.

In this setup, there are two parameters, $\lambda$ and the SVD truncation $r$, that have an impact on the smoothness of the solution in similar ways. The parameter $r$ determines the truncation in the decomposition of $\hat{X}$ that includes geometries of the entire ice sheet. Choosing a higher $r$ means that more spatial modes are used, representing variability on an increasingly small scale and corresponding to decreasing singular values (Fig. 3). The parameter $\lambda$ regulates how closely the observations are fitted in the least squares solution in Eqn (8). If the value for $\lambda$ is too low, the solution will be overfitting the observations and result in nonphysical solutions. In practice, these two parameters act in nearly equivalent ways (shown and discussed in Supplementary Material).

Cross-validation

To validate and compare the two reconstruction methods, we use the method of cross-validation: We leave out part of the observations, make the predictions in the same locations based on the remaining data and calculate the errors between predictions and observations. In the absence of independent datasets for validation, cross-validation is a commonly used method for approximating the prediction errors (e.g. Stone, Reference Stone1974; Hastie and others, Reference Hastie, Tibshirani and Friedman2003). Here, we follow a leave-block-out approach (Pohjankukka and others, Reference Pohjankukka, Pahikkala, Nevalainen and Heikkonen2017; Roberts and others, Reference Roberts2017) where we split the domain into approximately ten rectangular blocks measuring 200 km across and repeat the procedure described above, leaving out one block at a time. In such a way, we arrive at an error for each location where an observation is available. We denote this error by $\varepsilon_\text{SVD}$ and $\varepsilon_\text{GP}$ for the two respective methods:

(9)\begin{equation} \varepsilon_\text{method} = h_\text{rec} - h_\text{obs}, \end{equation}

where $h_\text{rec}$ is the elevation predicted by the respective reconstruction method without using the corresponding observation $h_\text{obs}$ as a constraint. We choose a block size of 200 km such that the points in the center of a block are roughly one correlation length scale away from the closest observation; therefore, we can validate the behavior of the GP, which is expected to perform more poorly at these large distances. For the SVD method, we split the entire dataset (both AeroDEM and ATM data) into these blocks. For the GP, we only exclude and reconstruct ATM data because, in contrast to the SVD method, the interior is the only part of the domain that is interpolated. However, we still use the AeroDEM data as constraints in each block interpolation.

In the case of the SVD method, we also use the cross-validation to determine optimal regularization parameters $\lambda$ and $r$. To this end, we calculate the map of errors $\varepsilon_\text{SVD}$ for a range of different values of both $\lambda$ and $r$ and choose the optimal parameter combination as the one with the smallest mean absolute error ( $r=1000$ and $\lambda = 10^7$, see Supplementary Material). In a strict sense, we are thus merging the processes of what is called validation (tuning of parameters) and testing (evaluation of the method) in the machine learning literature. However, in this case the regularization affects mostly the smoothness and small-scale structure of the solution rather than the larger-scale biases we are interested in detecting. Furthermore, we avoid having to split the dataset into three parts in a meaningful way between different iterations, rather than leaving out one block at a time.