# Uncertainty analysis of digital elevation models by
spatial inference from stable terrain

Romain Hugonnet, Fanny Brun, Etienne Berthier, Amaury Dehecq, Erik Schytt Mannerfelt, Nicolas Eckert, Daniel Farinotti

<figure><img src="image_1.png" /></figure>

## To cite this version:

Romain Hugonnet, Fanny Brun, Etienne Berthier, Amaury Dehecq, Erik Schytt Mannerfelt, et al. Uncertainty analysis of digital elevation models by spatial inference from stable terrain. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2022, 15, pp.6456-6472. 10.1109/JSTARS.2022.3188922. [hal-03719624](https://hal.science/hal-03719624v1)

HAL Id: hal-03719624

https://hal.science/hal-03719624v1

Submitted on 11 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Uncertainty analysis of digital elevation models by
spatial inference from stable terrain

Romain Hugonnet, Fanny Brun, Etienne Berthier, Amaury Dehecq, Erik Schytt Mannerfelt, Nicolas Eckert and Daniel Farinotti

**Abstract**—The monitoring of Earth's and planetary surface elevations at larger and finer scales is rapidly progressing through the increasing availability and resolution of digital elevation models (DEMs). Surface elevation observations are being used across an expanding range of fields to study topographical attributes and their changes over time, notably in glaciology, hydrology, volcanology, seismology, forestry and geomorphology. However, DEMs frequently contain large-scale instrument noise and varying vertical precision that lead to complex patterns of errors. Here, we present a validated statistical workflow to estimate, model, and propagate uncertainties in DEMs. We review the state-of-the-art of DEM accuracy and precision analyses, and define a conceptual framework to consistently address those. We show how to characterize DEM precision by quantifying the heteroscedasticity of elevation measurements, i.e. varying vertical precision with terrain- or sensor-dependent variables, and the spatial correlation of errors that can occur across multiple spatial scales. With the increasing availability of high-precision observations, our workflow based on independent elevation data acquired on stable terrain can be applied almost anywhere on Earth. We illustrate how to propagate uncertainties for both pixel-scale and spatial elevation derivatives, using terrain slope and glacier volume changes as examples. We find that uncertainties in DEMs are largely underestimated in the literature, and advocate that new metrics of DEM precision are essential to ensure the reliability of future Earth and planetary surface elevation assessments.

***Index Terms***—Geostatistics, error propagation, remote sensing, variogram, spatial correlation, surface height.

I. INTRODUCTION

DIGITAL elevation models (DEMs) are gridded, numerical representations of surface elevation. DEMs have a long history of interpolation from point measurements and digitized historical maps [1], [2]. Nowadays, DEMs are mostly generated from radar interferometry [3], [4], optical stereophotogrammetry [5], [6] or laser scanning [7], [8] of a planetary surface. When produced from these remote sensing techniques, DEM grid cells essentially represent surface elevations timestamped to the date of instrument acquisition. With the ever-improving coverage and precision of satellite and airborne sensors [9], land surface assessments based on DEMs

R. Hugonnet and E. Berthier are with the LEGOS, Université de Toulouse, CNES, CNRS, IRD, UPS, F-31400 Toulouse, France (e-mail: romain.hugonnet@gmail.com).

R. Hugonnet, A. Dehecq, E.S. Mannerfelt and D. Farinotti are with the Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zurich, Zurich, Switzerland, and the Swiss Federal Institute for Forest, Snow and Landscape Research (WSL), Birmensdorf, Switzerland.

F. Brun and A. Dehecq are with the IGE, Université Grenoble Alpes, CNRS, IRD, Grenoble INP, Grenoble, France.

N. Eckert is with the Université Grenoble Alpes, INRAE, UR ETNA, Grenoble, France.

Manuscript received March 23, 2022; revised June 10, 2022.

are advancing towards estimates that are both more spatially
and more temporally resolved [10], [11]. Additionally, the
recent unlocking of historical optical archives has created
unprecedented potential for studying half a century of Earth's
surface elevation [12]–[14].

Studies that harness elevation observations can generally
be divided into two groups. The first group relies on single-
acquisition and often gap-filled DEMs to extract essential
topographic characteristics, e.g., in river discharge and flood
modelling [15]–[17], geomorphological terrain analysis [18]–
[21], tectonic monitoring [22]–[25], avalanche risk prediction
[26], land classification [27], [28], onshore inundation and sea-
level rise forecasting [29]–[31] and planetary surface char-
acterization [32], [33]. The second group requires multiple
acquisitions to study surface elevation changes over time, e.g.,
for landslide and rock avalanche detection [34]–[36], seasonal
snow depth assessment [37]–[39], lava flow volume quantifica-
tion [40], [41], canopy height evolution [42]–[44] and glacier,
ice sheet and ice shelf mass balance estimation [45]–[47].
In both groups, and for all applications, the interpretation of
results and its robustness are inextricably intertwined with the
accuracy and precision of the underlying DEMs.

Accuracy and precision are related to systematic and ran-
dom errors. In the case of DEMs, they have been the focus
of specific research [48]–[51], software development [52] and
questioning [53]–[56] since the beginning of the numerical
era. Yet, these efforts are dwarfed by the tremendous increase
of studies that rely on DEMs [57] and the processing of ever
larger data volumes [58]–[60]. Most critically, the analysis of
many modern studies is still confined to simplified metrics for
accuracy and precision (e.g., [61]–[63]) that mix systematic
and random errors and fail to describe the strong spatial
variations and correlations in errors observed in DEMs (e.g.,
[64]–[66]).

Here, we present a statistical workflow to robustly estimate
and propagate uncertainties in DEMs; most specifically, we:

* perform a literature review of analyses dealing with DEM accuracy and precision;

* propose a framework based on spatial statistics to consistently address DEM accuracy and precision;

* present robust inferential methods to estimate elevation heteroscedasticity and spatial correlation of errors;

* analyze the impact on the uncertainty of elevation derivatives, using terrain slope and glacier volume changes as examples;

* provide access to our methods through the open, tested and documented Python package xDEM.
## II. LITERATURE REVIEW

### A. Mitigating poor DEM accuracy before studying precision

The term accuracy has been used to describe either systematic errors or, in some instances, both systematic and random errors, leading to some confusion. In the present article, we define accuracy as the description of systematic errors only, also known as "trueness" [67], which is related to elevation biases. Poor accuracy is common in DEMs and has been a major source of error in elevation assessments, particularly during the advent of space-borne DEMs. Limitations in instrument positioning, orientation, or post-processing often lead to erroneous horizontal referencing [66], [68], vertical shifts [69], [70] or tilts [71], [72] that propagate into elevation biases (Fig. 1a). By utilizing terrain with elevation assumed stable over time, methods performing 3-dimensional alignment of DEMs have flourished, relying on either generic registration methods [73]–[75], least squares approaches [71], [76] or specifically-developed DEM registration based on terrain constraints [77], [78]. These methods proved robust for aligning a DEM either to an external reference DEM, or to accurate geolocated point elevation data such as space-borne laser altimetry [79], [80]. The above registration methods are only successful at correcting elevation biases common to the entire DEM grid, however. Other biases remain present once 3-dimensional alignment is attained and can arise from resolution [81], [82], specific image deformations and instrument biases [13], [66] or physical properties of the observed terrain such as radar penetration into snow and ice [83], [84] or into forest canopy [43]. Most of these biases are instrument- or application-dependent and, therefore, require specific considerations. Notwithstanding those, poor DEM accuracy has been largely addressed by the robustness of registration methods that have become increasingly widespread, thereby shifting the focus towards the next limiting factor: better quantifying DEM precision.

### B. The inherent variability of vertical precision

Precision describes random errors [67] and is related to elevation variance. One aspect of DEM precision consists of the pixel-scale dispersion of elevations that we refer to as "vertical precision". DEMs are generated from acquisitions that possess intrinsic, random measurement errors. At the pixel scale, instrument resolution, spectral range, and encoding depth of optical sensors directly affect the quality of stereo-correlation [5], [6], [85], radar slant angle and height of ambiguity play an important role in interferometric coherence [86], [87] while laser wavelength, sunlight background radiation, target reflectivity, and backscattering properties modulate laser signal-to-noise ratio [88], [89]. Many instrument- or processing-related metrics constitute quality indicators of the estimated elevations. These indicators have been almost exclusively used for the filtering of observations of lesser quality, however, and only occasionally as a tool towards improved modelling of sensor-specific variability in vertical precision (e.g., [90]). Besides, the geometry of instrument acquisition can exacerbate random errors depending on the relief of the observed landforms (Fig. 1a). Vertical precision has indeed been long shown to decrease with terrain slope [48], [91]–[94].

Several assessments account for this variability by partitioning the elevation variance into categories of flat and steep terrain (e.g., [59]). Most studies use a single metric to describe vertical precision, however, often reporting a standard deviation (e.g., ±2 meters). Such simple metrics are insufficient in describing the heteroscedasticity of elevation measurements, i.e. the variability in vertical precision. Although some studies quantified and modelled this heteroscedasticity [64], [95], [96], this modelling was generally performed without validation of the underlying methodology and, most critically, without considering the effect of spatial correlations.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Patterns of random and systematic errors in DEMs. a, Elevation differences of horizontal shifts (left) and after alignment (right) with terrain hillshade (top), extracted from the data in Table I. The horizontal shift between Pléiades and SPOT-6 DEMs is of 2 m east and 4 m north, creating large biases despite being relatively small (half a pixel). b-c, Noise owed to (b) along-track undulations in a Pléiades–Pléiades DEM difference and to (c) digitization artefacts in a KH-9–ArcticDEM DEM difference [13], after alignment.</figcaption></figure>

### C. The correlated noises that plague DEMs

Another aspect of DEM precision concerns the inter-pixel spatial dependency of random errors, here referred to as "spatial correlations". Spatial correlations describe structures of noise that show a location-dependent pattern, which can often be traced back to limitations during acquisition or post-processing. Along-track undulations have been observed in many DEMs generated from air- and space-borne sensors (Fig. 1b), including the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [47], [66], the Satellite Pour l'Observation de la Terre (SPOT) [97], [98], Pléiades [39], [99] and the Shuttle Radar Topography Mission (SRTM) [58], [100], [101] (Fig. S1). Processing noise is common in DEMs requiring image digitization including aerial photographs [102], [103] or historical satellite imagery such as Corona and Hexagon KeyHole-9 (KH-9) [13], [104] (Fig. 1c). To mitigate these correlated noises, DEM correction methods have emerged [82], [105] but are still burgeoning for specific
types of errors [66], [106], [107], and their performance is highly dependent on the type of terrain. Furthermore, nearly all DEMs contain structural short-range correlation of errors. The degree to which a DEM grid spacing represents its native resolution [108], [109] and how that resolution has possibly been degraded through interpolation [110] determine the severity of these short-range correlations. When upsampled to a larger grid spacing, vertical precision improves directly as a function of the underlying spatial correlations [111]. Spatial correlations are generally quantified using an empirical variogram [112], [113] estimated either on the basis of differences with independent elevation observations [2], [114] or those with simulated elevation surfaces [115], [116]. Many studies have used variograms, but have almost exclusively used short range models (i.e. 5 to 20 times the pixel size). Few studies modelled longer-range correlations, that is, correlations that persist over distances several orders of magnitude larger than the pixel size [13], [47], [65]. The widespread occurrence of long-range noise in DEMs thus constitutes a critical limitation in the analysis of DEM precision, and one that directly affects uncertainty propagation.

## D. Uncertainty propagation to elevation derivatives

To propagate elevation variance into uncertainties of elevation derivatives (i.e. variables that are derived from elevations), a large set of methods has been applied that generally relies on spatial statistics. Spatial statistics, also known as geostatistics [112], [113], [117], provide a large body of theories and methods that, among others, can address spatial uncertainty analyses [118]–[120] by characterizing spatial correlations that depend only on the distance between observations. These uncertainty propagation methods can be subdivided into two groups: (i) Monte Carlo techniques that simulate multiple random realizations of correlated error fields [121]–[123], notably including Sequential Gaussian simulation [117] and Fourier randomization [124]; and (ii) gradient techniques that analytically approximate the variance of a derivative through simplified equations, that can be either based on Taylor series expansion [121] for any derivative of elevation, or approximations of variogram integration [65] for spatial derivatives. The first group has been widely used for topographic variables, notably in hydrology [16], [57], [125], [126] and occasionally for spatial derivatives in glaciology [127]. The second group is used less frequently, both for Taylor series expansions developed in few applications [128]–[131], and for variogram integration implemented mainly in glaciology and geomorphology [65], [132]. Although both groups are expected to perform similarly, Monte Carlo techniques are computationally expensive, especially at fine resolution. Analytical approximations, instead, require a theoretical description of variance propagation that can reach a high degree of complexity for some derivatives [133]. To our knowledge, few studies [132] constrained these propagation methods with estimates of heteroscedasticity and spatial correlation of errors into a single framework for DEMs, and none tested the underlying assumptions of spatial statistics. In the following, we propose such a framework, and later describe methods to robustly estimate its key components.

## III. PROBLEM FORMULATION

### A. Elevation bias and variance at each location

We consider the elevation observation $\hat{h}(x, y, t)$ located at $(x, y)$ in space and $t$ in time, and pertaining to the DEM $\mathcal{D}$. Annotating the true unknown elevation at the same location $h(x, y, t)$, we can state that the elevation observation has a bias $\delta h(x, y, t)$ if, over a large number of repeated measurements $i \in \mathcal{I}$ of elevation $\hat{h}(x, y, t)_i$ at $(x, y, t)$, we have:

$$ \overline{\hat{h}(x, y, t)}_{i \in \mathcal{I}} - h(x, y, t) = \delta h(x, y, t). \quad (1) $$

The repeat elevation measurements around the bias $\delta h(x, y, t)$ are subject to random measurement errors $\epsilon_h(x, y, t)$ with variance $\sigma_h^2(x, y, t)$, whose distribution is not necessarily normal and might depend on time and location:

$$ \hat{h}(x, y, t) = h(x, y, t) + \delta h(x, y, t) + \epsilon_h(x, y, t). \quad (2) $$

In practice, acquiring a large number of repeat measurements at both the same location and time is not feasible, and we therefore turn towards inferential methods to estimate these biases and variance.

### B. Inference from stable terrain

DEMs benefit from a great asset, largely uncommon to other remote sensing data, which is that large proportions of planetary surface elevations remain virtually unchanged through time. In fact, elevation changes caused by erosion, short vegetation growth, or continental drift are typically small compared to the precision of the measurement. Terrains such as bare rock or grasslands – later referred to as “stable terrain” – thus provide the means of analyzing multiple elevation measurements acquired at different points in time as if they were acquired from simultaneous measurements $\hat{h}(x, y, t)_i$:

$$ \frac{dh(x, y, t)}{dt} \approx 0 \text{ for } (x, y) \in \text{stable terrain}. \quad (3) $$

While this temporal consistency unlocks the potential to analyze elevation acquisitions independently of time $t$, it is impeded by the number of required DEMs. For each location $(x, y)$, the number of samples to perform the statistical analysis would always be at best equal to the total number of independent acquisitions, requiring a large number of DEMs. Therefore, we investigate the spatial properties of elevation biases and variance.

### C. Spatial homogeneity after affine alignment

Elevation biases and variance are inherent to instrumental limitations, to the physical properties of the observed terrain, as well as its topography (see previous Sections II-A and II-B). Among many types of location-specific biases, a general exception is that of grid misalignment to the true elevations $h(x, y, t)$ that follows specific geometric distributions linked to the gridded nature of DEMs (Fig. 1a). In our framework, we therefore split elevation biases into two categories: affine biases $\delta h_A$ that are common to the entire DEM (e.g., translation, rotation, scaling), and non-affine “specific” biases $\delta h_S$
<figure><img src="image_4.png" /><figcaption>Fig. 2. A framework for uncertainty analysis of DEMs. Non-stationary spatial framework to analyze the accuracy and precision of DEMs based on elevation differences on stable terrain (Sections III-A–III-D), with accuracy divided into affine and specific biases (Section III-C), and precision divided into heteroscedasticity and spatial correlation of errors (Section III-F).</figcaption></figure>

that occur at the grid cell level and vary with instrumental and topographical effects (Fig. 2):

$$
\delta h(x, y, t) = \delta h_A(x, y, t) + \delta h_S(x, y, t). \quad (4)
$$

Once an alignment is attained by the affine transformation
A giving $\mathcal{A}(x, y, t) = \delta h_A(x, y, t)$, we assume that, for a
single DEM $\mathcal{D}$, specific elevation biases $\delta h_S$ and elevation
variance $\sigma_h^2$ have a spatial distribution that is homogeneous
with the properties of the instrument and the observed terrain
$\mathcal{P}$. We use this spatial homogeneity to substitute space for
time. For example, we consider that elevations $h(x_1, y_1, t)$ and
$h(x_2, y_2, t)$ of $\mathcal{D}$ acquired on the same surface type (e.g., bare
rock), and under the same topographical attributes (e.g., flat)
will have similar specific biases and variance:

$$
\left.
\begin{array}{l}
\delta h_s(x_1, y_1) \approx \delta h_s(x_2, y_2) \\
\sigma_h(x_1, y_1) \approx \sigma_h(x_2, y_2)
\end{array}
\right\} \text{for } \mathcal{P}(x_1, y_1) = \mathcal{P}(x_2, y_2). \quad (5)
$$

Combining the assumptions of Eqs. 3 and 5, and provided
that we describe all the properties $\mathcal{P}$ of spatial homogeneity,
a large sample size can be used to infer $\delta h$ and $\sigma_h$ at each
location $(x, y)$ from a single difference between a DEM and an
independent source of elevation data. The properties of spatial
homogeneity $\mathcal{P}$ could differ between biases and variance. In
the following, we assume that specific elevation biases, if they
exist, are independently corrected and focus on characterizing
the elevation variance $\sigma_h^2$.

D. Elevation difference with an independent source

After performing affine alignment of elevations $\hat{h}_1(x, y, t_1)$ from a first source $\mathcal{D}_1$ and elevation $\hat{h}_2(x, y, t_2)$ of a second source $\mathcal{D}_2$, we subtract them to derive elevation differences $dh_{1-2}(x, y)$. Assuming independence between the error of each elevation source, the variance of the difference is:

$$
\sigma_{dh_{1-2}}^2(x, y) = \sigma_{h_1}^2(x, y) + \sigma_{h_2}^2(x, y). \quad (6)
$$

By selecting a second source to observe $\hat{h}_2$ that is of higher
precision than the first source that observes $\hat{h}_1$, the analysis of
the differences $\hat{h}_2 - \hat{h}_1$ will largely capture the variance of the
first source. For example, if the second source is three times
more precise than the first, Eq. 6 implies that about 95% of
the variance of the elevation difference will originate from the
first source, yielding:

$$
\sigma_{h_1}(x, y) \approx \sigma_{dh_{1-2}}(x, y). \qquad (7)
$$

Alternatively, if $h_1$ and $h_2$ originate from independent acquisitions of the same instrument and processing, we have:

$$
\sigma_{h_1}(x, y) = \frac{\sigma_{dh_{1-2}}(x, y)}{\sqrt{2}}. \qquad (8)
$$

Thus, we use elevation differences to infer on $\sigma_h$, which
can be converted from either Eqs. 7 or 8.
E. Discriminating elevation bias from variance in spatial statistics

To further analyze elevation variance, we need to discriminate bias from variance. When analyzing elevation differences, what appears as a bias at the local scale could also be a form of long-range correlation at larger scales (Fig. 1b-c). This distinction is directly related to the assumption of second-order stationarity of spatial statistics. For elevation differences, second-order stationarity implies that the following assumptions should be fulfilled (see Supplementary Section II-A):

1) a first assumption of stationary mean, i.e. that the average of elevation differences dh(x, y) is constant over large areas;

2) a second assumption of stationary variance, i.e. that the variance of elevation differences $\sigma_{dh}(x, y)$ is constant over large areas; and

3) a third assumption of spatially consistent covariance, i.e. that the correlation between random errors of elevation differences only depends on the distance between observations.

Large areas here refer to areas slightly smaller than the size of the study domain, typically within an order of magnitude. As such, a correlated error with a correlation range that is orders of magnitude larger than the size of the study domain might be considered a vertical bias common to the entire DEM grid (Fig. 2). And, inversely, such a bias placed in the context of a larger study domain might be considered as a correlated error, if the elevation differences fulfill the above assumptions.

Thanks to the affine alignment of our elevation differences, we verify the first assumption of stationary mean. However, the heteroscedasticity of elevations (see Section II-B) invalidates the second and third assumptions, and therefore a non-stationary framework needs to be defined.

F. A non-stationary spatial framework for DEM analysis

To perform spatial statistics with a non-stationary variance, transformation of the data towards a stationary variance is necessary. The transformation depends on the nature of the spatial variability and correlations. In DEMs, we identify two types of correlation: short-range ones related to resolution, and long-range ones related to correlated noise or digitization artefacts. While the latter appear unrelated to the heteroscedasticity of elevation, the former are similarly linked to local instrument- and terrain-dependent variables (see Sections II-B and II-C). We thus subdivide elevation variance into elevation heteroscedasticity and spatial correlation of errors (Fig. 2) assuming that longer-range correlations are independent of elevation heteroscedasticity, which yields:

$$\sigma_{dh}^2(x, y) = \sigma_{dh_{sr}}^2(x, y) + \sigma_{dh_{lr}}^2, \quad (9)$$

where $\sigma_{h_{sr}}^2(x, y)$ is the variable short-range variance at $(x, y)$, $\sigma_{h_{lr}}^2$ is the constant long-range variance.

Using the variable spread $\sigma_{dh}(x, y)$, the elevation differences can be standardized into a standard score $z_{dh}$ with unit variance, which fulfills the second assumption of second-order stationarity:

$$z_{dh}(x, y) = \frac{dh(x, y)}{\sigma_{dh}(x, y)}. \quad (10)$$

Additionally, the spatial covariance $C_{z_{dh}}$ of $z_{dh}$, related to the variogram $\gamma_{z_{dh}} = 1 - C_{z_{dh}}$, is also free of the influence of heteroscedasticity and now fulfills the third assumption of second-order stationarity:

$$\gamma_{z_{dh}}^2(d) = \left( \frac{\overline{\sigma_{dh_{sr}}|_D}}{\sigma_{dh}} \right)^2 \gamma_{sr}(d) + \left( \frac{\overline{\sigma_{dh_{lr}}}}{\sigma_{dh}} \right)^2 \gamma_{lr}(d), \quad (11)$$

where $d$ is the spatial lag, i.e. the distance between two given observations, $\overline{\sigma_{dh_{sr}}|_D}$ is the average of $\sigma_{dh_{sr}}$ in the DEM $D$, and $\gamma_{sr}$ and $\gamma_{lr}$ are the short- and long-range variogram functions.

With all the assumptions in our framework fulfilled, we can now reliably use spatial statistics for uncertainty propagation. To this end, we require an estimate of the elevation dispersion $\sigma_{dh}(x, y)$ and of the variogram of the standard score $\gamma_{z_{dh}}(d)$, which describe the heteroscedasticity and the spatial correlation of errors, respectively. We also need to ensure that our assumption of spatial homogeneity remains valid when using stable terrain as an error proxy to infer heteroscedasticity and spatial correlations on moving terrain. In the following, we address these aspects by utilizing near-simultaneous data and implementing robust methods.

IV. DATA

A. Mont-Blanc case study: simultaneous DEMs

To demonstrate the methods associated with our proposed framework, we present a case study of two DEMs generated one day apart in the Mont-Blanc massif, French Alps (Fig. 3b, Table I). These DEMs were produced with a spatial posting of 5 m from SPOT-6 and Pléiades stereo images using the Ames Stereo Pipeline [134]. We utilize the temporal closeness of the two acquisitions to assess if stable terrain can be used as a proxy for moving terrain, considering a negligible elevation change on moving terrain.

We present an additional case study in the Northern Patagonian Icefield to illustrate the influence of the quality of stereo-correlation, a sensor-dependent variable, on elevation heteroscedasticity (Supplementary Section I-A with additional refs. [66], [135], [136]). This case study is based on simultaneously acquired ASTER [47] and SPOT-5 images. Furthermore, the DEMs used to illustrate noise patterns (Figs. 1 and S1) are described in the Supplementary Section I-B with additional refs. [137]–[139].

<table><caption>TABLE I NEARLY-SIMULTANEOUS PLÉIADES AND SPOT-6 DEMS USED FOR THE MONT-BLANC CASE STUDY.</caption><thead><tr><th>Instrument</th><th>Acquisition time</th><th>Resolution of stereo-pair</th></tr></thead><tbody><tr><td>Pléiades</td><td>24/10/2017, 12:00 CET</td><td>1.5 m</td></tr><tr><td>SPOT-6</td><td>25/10/2017, 12:30 CET</td><td>0.7 m</td></tr></tbody></table>
## B. Inventory and land cover products

We define moving terrain as glacierized, forested and seasonally snow-covered terrain, and exclude water bodies from our analysis. The remaining terrain is assumed to be stable. We mask glaciers using the Randolph Glacier Inventory 6.0 (RGI 6.0) outlines [140], which are delineated from images with a typical resolution of 15–30 m. We mask forests and water bodies using the ESA Climate Change Initiative Land Cover version 2.0.7 [141] which has a resolution of 300 m. Forested terrain corresponds to either broadleaved, needle-leaved, evergreen, or deciduous tree cover classes.

We identify specific elevation biases over forested terrain between the SPOT-6 and Pléiades DEMs – likely owing to different native resolution, orientation and spectral bands (Fig. S3a) – and thus exclude this terrain from our analysis. Our end-of-summer acquisitions contain little snow outside of glacierized surfaces. Therefore, we did not mask off-ice snow cover. Ultimately, in our analysis, moving terrain corresponds to glacierized terrain.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Mont-Blanc case study. Hillshade of Pléiades DEM and land cover for the Mont-Blanc case study. Example glaciers serve to illustrate Section VI-C.</figcaption></figure>

# V. METHODS

## A. Robust statistics and alignment

We use the median instead of the mean as a robust estimator of central tendency, and the normalized median absolute deviation (NMAD) instead of the standard deviation as a measure of statistical dispersion. Both choices are to mitigate the effects of frequent outliers in DEMs [142]. Combining these estimators with the dense sampling of stable terrain also ensures robustness to elevation changes of potentially

unmasked moving terrain. This includes rare events such as landslides or ground subsidence, or events that can occur over a small portion of the analyzed terrain such as volcanic uplift or sediment transport. We coregister DEMs on stable terrain for horizontal and vertical shifts following the aspect-slope relation described in [77] and we correct for possible tilts through least squares optimization of a plane [71].

## B. Heteroscedasticity

We estimate elevation heteroscedasticity by sampling an empirical dispersion of elevation differences $\hat{\sigma}_{dh}$ using the NMAD of binned categories along the terrain slope $\alpha$ [143] and the terrain maximum absolute curvature $c$ (Figs. 4a-c and S4). Maximum absolute curvature is defined as the maximum of the absolute profile curvature and the absolute planform curvature at each location [144]. All terrain attributes are estimated from the Pléiades DEM that contains the least data gaps. When available, the binning can also include an instrument quality factor $q$, such as the quality of stereo-correlation (Supplementary Section I-A, Figs. S5–S7) or interferometric coherence.

We numerically model the empirical dispersion $\hat{\sigma}_{dh}$ as a function $\sigma_{dh}$ of the terrain- and sensor-dependent variables $(\alpha, c, q)$ by multidimensional linear interpolation of the binned data (Fig. S8). The modelling of this variability can also be performed by fitting parametric models, for example an exponential model with the slope or a linear model with the maximum curvature (Fig. S9). These are more robust in the case of small sample sizes of elevation differences.

We standardize the elevation differences $dh$ following Eq. 10, using the modelled dispersion $\sigma_{dh}(\alpha, c, q)$:

$$z_{dh} = \frac{dh}{\sigma_{dh}(\alpha, c, q)} \quad (12)$$

After standardization, we verify that the standard score of the elevation differences matches a normal distribution by quantile-quantile plotting, and by comparison to a normal distribution fit [142], [145] (Fig. S10). The substantial improvement validates our choice of terrain slope and maximum curvature as key variables to describe elevation heteroscedasticity, as those largely explain the departure of random elevation errors from normality.

## C. Spatial correlations

We estimate spatial correlations by sampling an empirical variogram $\hat{\gamma}$ on the standard score $z_{dh}$ using Dowd's estimator [146], [147] (Fig. 5a):

$$2\hat{\gamma}_{z_{dh}}(d) = 2.198 \cdot \text{median}(z_{dh}(x, y) - z_{dh}(x', y'))^2, \quad (13)$$

where $z_{dh}$ is the standard score of elevation differences, and locations $(x, y)$ and $(x', y')$ are separated by a spatial lag $d$.

Dowd's estimator is based on median absolute deviations, and consequently more robust than the Matheron [148] or Cressie-Hawkins [149] estimators classically used (see Supplementary Section II-B based on additional ref. [150]). We
<figure><img src="image_6.png" /><figcaption>Fig. 4. Heteroscedasticity inference from stable terrain as function of slope and curvature for the Mont-Blanc case study. a-b, Violin plots of elevation differences on stable and moving terrain by bins of (a) slope and (b) maximum curvature. Dispersion inferred from stable terrain is showed by a thick line with color matching other panels. Note the logarithmic scales of histograms. c, Heatmap of stable terrain dispersion for slope and maximum curvature. Bins with a relative dispersion difference between stable and moving terrain greater than 30% (dark gray and black dots) contain less than 1% of samples. d, Inferred spatial distribution of vertical precision for all terrain, with inset that matches Fig. 1a.</figcaption></figure>

verify the increased robustness of Dowd’s estimator for the Mont-Blanc case study (Figs. S11 and S12).

To improve the variogram estimation, we introduce a pair-wise subsampling method based on iterative subsetting of pair-wise combinations between a disk and multiple rings centered on a random point (see Supplementary Section II-C). As variograms were historically sampled from point measurements [112], traditional sampling methods are less computationally efficient on large grids. Most critically, they are inefficient at sampling pairwise distances evenly across spatial scales, which is substantially improved by our method to estimate more reliably both short-range and long-range correlations (Fig. S13). Finally, we derive empirical variograms for 100 independent realizations with the same binning. We estimate our final empirical variogram by the mean of all realizations at each spatial lag with, as an empirical uncertainty, the standard error of the mean.

To derive a spatially continuous representation of the variogram, we calibrate an analytical model $\gamma_{z_{dh}}$ with the empirical variogram $\hat{\gamma}_{z_{dh}}$. We fit a sum of $k$ variogram models $V(s_k, r_k, d)$, optimizing their partial sills $s_k$ (i.e. correlated variance) and ranges $r_k$ (i.e. correlation length) simultaneously by weighted least squares, using as weights the squared inverse of the empirical uncertainties previously detailed (Fig. 5a):

$$ \gamma_{z_{dh}}(d) = \sum_{k} V_k(s_k, r_k, d). \quad (14) $$

For the Mont-Blanc study, we find no significant improvement in least-squares residuals when fitting more than three models, which are capable of capturing one short-range and two long-range correlations (Fig. S14, Table S2). The two long-range correlations match the along- and cross-track lengths of low-amplitude undulations in the elevation differences (Fig. 5b). We thus use three models to avoid the possible overfitting of a larger number of summed models. Generally, $k$ should be chosen to reflect the number of distinct

ranges in the patterns of DEM noise. For instance, ASTER undulations are characterized by two wavelengths of 1–2 km and 5–10 km in the along-track direction, and a cross-track distance of 60 km (Fig. S1a), which better fits three distinct long-range models [47] for a total of four ranges.

We also identify a low sensitivity to different variogram model types (Fig. S15, Tables S3 and S4), which shows that adequately modelling the multi-range nature of the spatial correlations is more important than refining that of their spatial form (Fig. 5b). For the Mont-Blanc case study, we reached the smallest least-squares residuals using a gaussian model $G(s, r, x - x')$ at short ranges, and spherical models $S(s, r, x - x')$ at long ranges [151] and used those henceforth:

$$ G(s, r, d) = s \left(1 - e^{-\left(\frac{2d}{r}\right)^2}\right), \quad \text{and} \quad (15) $$

$$ S(s, r, d) = \begin{cases} s \left( \frac{3}{2} \frac{d}{r} - \frac{1}{2} \left( \frac{d}{r} \right)^3 \right) & \text{if } 0 < d < r, \\ s & \text{if } d \ge r. \end{cases} \quad (16) $$

### D. Uncertainty propagation

1) **Simulation methods for elevation derivatives:** For derivatives of elevation with a complex spatial gradient, such as terrain slope and aspect later analyzed, we use simulation methods. We find similar results using Fourier randomization [124], [152] and unconditional Gaussian simulation [117], [153], and thus only use the former in the following. For 1,000 realizations, we simulate a random correlated error field of the standard score $z_{dh}$ based on the modelled spatial correlation $\gamma_{z_{dh}}$ in Eq. 14. We then de-standardize $z_{dh}$ using Eq. 12, and add the resulting elevation error field to the studied DEM. For each of these DEM realizations with an added error field, we then compute the terrain attribute of interest (e.g., terrain slope or aspect), for which we can study the distribution of errors.
<figure><img src="image_7.png" /><figcaption>Fig. 5. Spatial correlation inference from stable terrain for the Mont-Blanc case study. a, Spatial variogram of standardized elevation differences on stable and moving terrain. The empirical variogram is based on Dowd's estimator [146] and modelled by either a short-range spherical model or the sum of short- and long-range spherical models. b, Excerpt of the standardized elevation difference map, which highlights the correlated signals at both short-range (30 meters) and long-range (3.9 km in the along-track direction, 11.2 km in the cross-track direction). c, Standardized elevation uncertainty with increasing circular averaging area validated by empirical Monte Carlo sampling. Note the logit scale of the Y-axis.</figcaption></figure>

2) Theoretical approximation methods for spatial derivatives: For spatial derivatives such as the average $\overline{dh}$ of elevation changes $dh$ in an area $A$, we derive an exact analytical solution of the uncertainty in the spatial average $\overline{\sigma_{dh}}$:

$$ \sigma_{\overline{dh}}^2 = \frac{1}{N^2} \sum_{i=1}^{N} \sum_{j=1}^{N} \rho_{ij} \sigma_{dh_i} \sigma_{dh_j}, \quad (17) $$

where $N$ denotes the number of samples $i$ falling in the area $A$, $\sigma_{dh_i}$ is the vertical precision of pixel $i$, and $\rho_{ij} = (1 - \gamma_{z_{dh}}(d))$ is the spatial correlation between pixel $i$ and pixel $j$ based on their distance $d$.

In practice, Eq. 17 raises the issue of scaling exponentially with the number of samples, possibly resulting in trillions of calculations. To remedy this, we propose an approximation for spatially contiguous areas, inspired by the approach of [65] that computes a single aerial integral by approximating the area $A$ by a disk of the same area. Here, for each pixel $k$ of a random subset of $K$ pixels within the $N$ pixels, we compute the single aerial integral of the variogram numerically. We then approximate the variogram integral by the average of these subset aerial integrations (see Supplementary Section II-E):

$$ \sigma_{\overline{dh}}^2 \approx \overline{\sigma_{dh}}^2 |_{A} \frac{1}{N} \frac{1}{K} \sum_{k=1}^{K} \sum_{i=1}^{N} (1 - \gamma_{z_{dh}}(x_k - x_i)), \quad (18) $$

where $\overline{\sigma_{dh}}^2|_A$ is the average variance of the elevation differences of pixels $i$ in the area $A$:

$$ \overline{\sigma_{dh}^2}|_A = \frac{1}{N} \sum_i \sigma_{dh_i}^2. \quad (19) $$

We show that our method improves the accuracy of the theoretical approximation of [65] by accounting for more complex area shapes than disks while maintaining computational efficiency (Fig. S16). Additionally, these formulations can be linked to a number of effective samples, which describes the number of samples among the $N$ pixels in area $A$ that are statistically independent based on the spatial correlations modelled by $\gamma_{z_{dh}}$ (see Supplementary Section II-D).

Once uncertainties have been integrated from a spatial support (e.g., pixels) to a larger spatially contiguous ensemble (e.g., glaciers), they can be propagated again to a larger ensemble (e.g., all glaciers in a region) following Krige's relation of transitivity [112], [154]. For this, Eq. 17 can be applied for each pair of spatially contiguous ensembles $i$ and $j$ of area $A_i$ with the same variogram $\gamma_{z_{dh}}$ composed of the $k$ summed models $V_k(s_k, r_k, d)$:

$$ \sigma_{\overline{dh}} = \frac{1}{(\sum_i A_i)^2} \sum_k \sum_i \sum_j \left( \sigma_{\overline{dh}_{k,i}} \sigma_{\overline{dh}_{k,j}} - V_k(\sigma_{\overline{dh}_{k,i}}, \sigma_{\overline{dh}_{k,j}}, r_k, d_{i-j}) \right) A_i A_j, \quad (20) $$

where $d_{i-j}$ is the distance between the centroids of ensemble $i$ and $j$, and $\sigma_{\overline{dh}_{k,i}}$ is the spatially integrated uncertainty
of ensemble *i* associated to the variogram model $V_k$, partial sill $s_k$ and range $r_k$ with pixel pairs *n* and *m* (see Eq. 17):

$$ \sigma_{\overline{dh}_{k,i}} = \frac{1}{N^2} \sum_{n=1}^{N} \sum_{m=1}^{N} (s_k - V_k(s_k, r_k, x_n - x_m)) \sigma_{dh_n} \sigma_{dh_m}. \quad (21) $$

Furthermore, we use a Monte Carlo spatial sampling method to validate our uncertainties of spatially averaged elevations, thus indirectly verifying the robustness of our modelled spatial correlations (Fig. 5c). We randomly sample up to 10,000 circular patches of area *A* without replacement. We compute the mean $\overline{dh}$ inside circular patches, keeping only those with more than 80% valid elevation differences *dh* to mitigate the effects of missing data. We use the NMAD of 10,000 realizations to empirically estimate the uncertainty of the spatially averaged *dh* of area *A*, and repeat this procedure for varying area sizes *A* (Fig. 5c). This method substitutes repeated correlated simulation of Fourier randomization or Gaussian simulation by a repeated spatial sampling, relying on the assumption of spatial homogeneity of variance on stable terrain (Section III-C). As it requires a large number of independent patches to produce a robust estimate, the area size *A* for which it can estimate an uncertainty is limited to sizes much smaller than that of the spatial domain. It is also highly dependent on the availability of stable terrain. Therefore, we use it only for validation purposes.

## VI. RESULTS AND DISCUSSION

In Section VI-A below, we discuss the use of stable terrain as an error proxy based on the methods applied to the Mont-Blanc case study. In Sections VI-B and VI-C, we then analyze the impacts of heteroscedasticity and spatial correlations when propagating elevation variance into uncertainties of pixel-scale elevation derivatives such as terrain slope, or spatial derivatives such as glacier volume changes. In those two sections, we provide examples based on the Mont-Blanc case study and determine the impact of our methods for a set of assumptions on the variance properties during uncertainty propagation:

* either homoscedastic elevation (constant variance, shortened "homosc.") or heteroscedastic elevation (variable variance, "heterosc."); and

* either no spatial correlation (shortened "no corr."), or only short-range correlations ("short-range"), or both short- and long-range correlations ("long-range").

In this exercise, the most realistic case refers to the one that accounts for potential elevation heteroscedasticity and potential short- and long-range correlations. Uncertainties are reported as a symmetric confidence interval of $1\sigma$ (68% confidence level) or $2\sigma$ (95%), specified in each case.

### A. Validation of stable terrain as an error proxy

We test the validity of using stable terrain as a proxy of elevation errors for moving terrain on the nearly simultaneous DEMs of the Mont-Blanc case study. We find that elevations on moving terrain exhibit the same heteroscedasticity with

slope and curvature than those on stable terrain, with less than 1% of binned samples that differ by more than 30% (Fig. 4a-c). We additionally verify that this elevation heteroscedasticity is continuous between neighbouring bins when using robust estimators, thereby consolidating our assumption of spatial homogeneity (Section III-C, Eq. 5). By extending this assumption to the case of moving terrain, we infer a complete map of vertical precision (Fig. 4d).

We find similar spatial correlations of errors between stable and moving terrain (Fig. 5a). Values of partial sills and ranges of the variogram models that describe these correlations are within the same orders of magnitude (Table II), despite greater differences at long ranges due to the limited pairwise samples available on moving terrain. Using our Monte Carlo sampling method, we validate the increased robustness of using multiple correlation ranges to estimate uncertainties across spatial scales (Fig. 5c). Our results indicate that using a short-range model alone underestimates elevation uncertainties by several orders of magnitude for areas larger than 0.1 km².

<table><caption>TABLE II
ESTIMATED VARIOLGRAM MODEL PARAMETERS FOR THE SPATIAL CORRELATIONS OF ELEVATION ERRORS IN THE MONT-BLANC CASE STUDY. STABLE AND MOVING TERRAIN IS DISTINGUISHED (COLUMNS). GAUSSIAN COMPONENTS ARE LISTED FOR THE SHORT-RANGE MODEL AND SPHERICAL COMPONENTS FOR THE LONG-RANGE ONES, AS IN FIG. 5. PARTIAL SILLS ARE EXPRESSED AS A PERCENTAGE OF THE TOTAL VARIANCE.</caption><thead><tr><th>Model parameters</th><th>Stable</th><th>Moving</th></tr></thead><tbody><tr><td>Sill of short-range model</td><td>93%</td><td>95%</td></tr><tr><td>Range of short-range model</td><td>30 m</td><td>38 m</td></tr><tr><td>Sill of 1<sup>st</sup> long-range model</td><td>2%</td><td>4%</td></tr><tr><td>Range of 1<sup>st</sup> long-range model</td><td>3,900 m</td><td>2,400 m</td></tr><tr><td>Sill of 2<sup>nd</sup> long-range model</td><td>5%</td><td>1%</td></tr><tr><td>Range of 2<sup>nd</sup> long-range model</td><td>11,200 m</td><td>10,800 m</td></tr></tbody></table>

For elevation heteroscedasticity, our results highlight the importance of elevation standardization to ensure an adequate scaling when inferring on another type of terrain (e.g., from steep, stable terrain to flat, moving terrain). Yet, our analysis only exemplifies snow- and ice-covered terrain with high-resolution stereophotogrammetric DEMs. The physical properties of the observed terrain in relation to the utilized sensor might in some cases invalidate our assumption of spatial homogeneity. For instance, we found that our standardization did not mitigate the larger errors of elevation over forested areas (Fig. S3a). In such a case, an upfront investigation of specific elevation biases is required. After these biases are corrected, a refined modelling of elevation heteroscedasticity based on sensor-dependent variables can help to reach a good description of the properties of spatial homogeneity. We indeed found a strong relationship with the quality of stereo-correlation for the case study of the Northern Patagonian Icefield (Figs. S6 and S7). The rougher resolution (15 m) and spectral range (8 bits) of the ASTER stereo images, compared to those of SPOT-6 and Pléiades (metric resolution and 12-bits), leads to a significant variability in elevation errors with terrain texture.
For spatial correlations, we highlight the value of standard-
ization to reduce variability for empirical variogram estimation
(Figs. S3b-d and S12). It is especially useful to deconvolve
the long-range correlations with small magnitude to the short-
range ones. Heteroscedasticity may indeed explain the short-
range variogram anisotropy found by previous studies [155].
We nevertheless identify a slight difference in the well-
constrained short correlation range between stable and moving
terrain (30 m vs 38 m, respectively; Table II). This difference
might be due to the rougher interpolation of stereophotogram-
metric block-matching algorithms over bright, lower-texture
glacierized terrain. In some cases, sensor properties or pro-
cessing schemes influence not only the magnitude of spatial
variability but also the scale of correlations. Developing a
statistical framework that continuously includes these effects
might be overly complex for most analyses that, instead, could
adjust estimates of short-range correlation depending on the
type of observed terrain.

We conclude that stable terrain is a valid proxy for error
analysis, provided that elevation heteroscedasticity is taken
into account. However, the quality of statistical inference from
this error proxy depends directly on the number of stable
terrain samples available. For some DEMs, these samples
might be scarce in the proximity of continuous expanses of
moving terrain (e.g., at the margins of ice sheets or large
forests) and thus insufficient to perform robust inference. To
address this, the stable terrain of independent DEMs, possibly
located elsewhere, could be utilized if they are generated
from the same instrument and processing chain. Many DEMs
indeed have consistent error properties between segments
acquired under similar conditions around the world (e.g., [47],
[59], [156]). For instruments with correlated noise of varying
amplitude, such as Pléiades or ASTER, long-range correlations
can be more robustly inferred from a multiple-acquisition
average of variograms.

B. Impact on pixel-scale derivatives of elevation: example
with terrain slope and aspect

We illustrate the propagation of elevation uncertainty to the
slope and aspect in a 4 km² area around the Mont-Blanc sum-
mit (Fig. 6a). We select this area due to its wide range of slopes
and aspects, and its small extent facilitating computationally
expensive simulations. To avoid the circularity of the aspect
variable when assessing uncertainty, we divide it into northness
(i.e. cosine of the aspect) and eastness (i.e. sine of the aspect)
which denote, respectively, the north-south and east-west tilt
of the slope.

We propagate uncertainties in the Pléiades DEM by sim-
ulating random elevation error fields (see Section V-D) for
every set of assumptions (Fig. S17). For this example, we
assume that SPOT-6 and Pléiades have random errors of
similar amplitude, and estimate the random errors of the
Pléiades DEM following Eq. 8. We generally note a strong
deviation from normality and asymmetry in the simulated
uncertainty distribution of terrain attributes (Fig. S18). While
this asymmetry requires specific considerations for in-depth
terrain analysis, we here provide a simplified picture by

estimating a symmetric 1-σ uncertainty derived from the half-
difference between the 16<sup>th</sup> and 84<sup>th</sup> percentile of the simulated
slope, northness or eastness of each pixel.

Our analysis reveals that elevation heteroscedasticity plays
a major role in the spatial distribution of uncertainties in
slope and aspect. In particular, it exacerbates errors in steep
and rough terrain. Spatial correlations moderately affect un-
certainties by slightly reducing their amplitude (Fig. 6b-c).
We interpret the latter to be due to an increase in the spatial
coherence of terrain derivatives when the elevation errors are
spatially correlated. Since topographical attributes are derived
over a 3x3 pixel window, the closer the short-range spatial
correlations are to a 3-pixel length, the larger the impact on
the amplitude change (Fig. S19).

By aggregating uncertainties into slope categories, we show
that uncertainties in flat terrain are overestimated when as-
assuming homoscedasticity and no spatial correlation, while
those in steep terrain are underestimated by up to a factor

<figure><img src="image_8.png" /><figcaption>Fig. 6. Uncertainty propagation to terrain slope and aspect at the Mont-Blanc summit. a, Hillshade and terrain attributes based on the Pléiades DEM. b, Slope and aspect uncertainty estimated by the half-difference between the 16^th^ and 84^th^ percentiles of 1,000 simulated terrain attributes at each pixel. c, Distributions of slope and aspect uncertainties by category of terrain slope for each set of assumptions, with boxes denoting the interquartile range and whiskers extending to the entire distribution.</figcaption></figure>
of 10 (Fig. 6c). Slope uncertainties decrease near slopes of 90 degrees, likely because elevation errors tilt the terrain in different orientations while generally maintaining a steep slope, which translates into aspect uncertainties. We reach similar conclusions when aggregating uncertainties by maximum absolute curvature categories, our second variable that describes elevation heteroscedasticity (Fig. S20).

### C. Impact on spatial derivatives of elevation: example with glacier volume changes

<figure><img src="image_9.png" /><figcaption>Fig. 7. Uncertainty propagation to glacier mean elevation changes at the Mont-Blanc massif. a-b, Distributions of uncertainty of glacier mean elevation change by category of area and average terrain slope for each set of assumption, with boxes denoting the interquartile range and whiskers extending to the entire distribution. c, Empirical evaluation of uncertainty ranges for mean glacier elevation changes in the Mont-Blanc case study. Correct uncertainty estimates should cross the vertical zero line in 95% of the cases.</figcaption></figure>

We consider 84 glaciers in the Mont-Blanc massif that have at least 85% of their area covered by valid elevation differences. We analyze the mean elevation changes within the outline of each glacier, which can be converted to volume changes after multiplication by the glacier area, and propagate uncertainties for each set of assumptions.

We find that spatial correlations strongly hamper the decrease in uncertainty with increasing glacier area (Fig. 7a). Long-range correlations are the main contributor to uncertainty for large areas, mirroring the validation of Fig. 5c. While

our case study has long-range correlations of only 7% of the variance, uncertainties of mean elevation changes for glaciers larger than 10 km² are underestimated by a factor of about 25 when based solely on short-range correlation. This is striking, and even more so when realizing that the underestimation is nearly by a factor of 150 when totally omitting spatial correlations. This dramatic increase is explained by the fact that long-range correlations essentially correspond to local biases.

Heteroscedasticity has a moderate influence on the uncertainty of each glacier, impacting its amplitude by a factor of 1 to 3. The uncertainty of glaciers located in flat areas is overestimated when using a homoscedastic assumption due to the larger average variance over rougher, stable terrain. On the contrary, the uncertainty of the steepest glaciers is underestimated (Fig. 7b). Using the empirical comparison provided by the nearly simultaneous volume changes (Fig. 7c), we show that the uncertainties for the mean elevation change are most realistic when accounting for long-range spatial correlation. In such a case, 89% of the ranges intersect zero (the true volume change) at the 2σ level (i.e. 95% confidence), in contrast to only 30% for short ranges and 7% for no correlation. Yet, our uncertainties are slightly too low.

We identify the cause of this underestimation as the omission of a longer-range correlation close to the size of the DEM and thereby difficult to constrain. This longer-range correlation arises from the fact that along-track undulations are fully correlated in the cross-track direction with 20 km swath. Directional variography could help characterize such correlations, but would lead to a more difficult uncertainty propagation, with exacerbated complexity when combining several DEMs. Instead, we maintain an omnidirectional variogram to describe correlations, but assess a conservative estimate based on artificial undulations (Fig. S21). This results in the replacement of the 11.2 km correlation with a 20 km one (DEM swath width) and a partial sill twice larger. We then find that 93% of the uncertainties for glacier larger than 0.2 km² intersect zero at the 95% confidence level, confirming the increased robustness with these considerations. Only 87% do so for smaller glaciers, however. This discrepancy might be explained by unaccounted heteroscedasticity from landform-projected shadows that particularly affects small glaciers in steep and north-facing slopes.

<table><caption>TABLE III
SPATIAL UNCERTAINTY PROPAGATION AMONG GLACIERS IN THE MONT-BLANC MASSIF. MEAN ELEVATION CHANGE UNCERTAINTIES $\sigma_{dh}$ ARE PROPAGATED USING EQ. 20.</caption><thead><tr><th>Uncertainty of spatial average</th><th>Griaz and Bourgeat</th><th>Bossons and Taconnaz</th><th>All glaciers</th></tr></thead><tbody><tr><td>σ<sub>dh</sub> no corr. (m)</td><td>0.027</td><td>0.006</td><td>0.002</td></tr><tr><td>σ<sub>dh</sub> short-range corr. (m)</td><td>0.027</td><td>0.006</td><td>0.002</td></tr><tr><td>σ<sub>dh</sub> long-range corr. (m)</td><td>0.061</td><td>0.049</td><td>0.024</td></tr><tr><td>Glacier area (km²)</td><td>0.9</td><td>16.3</td><td>131.1</td></tr></tbody></table>

When uncertainties of volume change of several glaciers are propagated into that of the massif, correlations also come into play. We illustrate the propagation at different spatial scales
by considering several glacier groups: one group with the two
small neighbouring glaciers of Griaz and Bourgeat, another
group with the two large neighbouring glaciers of Bossons
and Taconnaz, and a third group including all 84 glaciers
(Fig. 3). We find identical uncertainties when considering
no correlation, or only short-range correlations (Table III).
This reflects the fact that all glaciers are separated by at
least 30 m, i.e. a distance larger than that of our short-range
correlation (Table II). Long-range spatial correlations have a
large impact on the total uncertainty, however, with a tenfold
underestimation of the uncertainty for all glaciers in the
massif when omitting them. Increased uncertainties from long-
range correlations mostly affect large neighbouring glaciers,
as shown for Bossons and Taconnaz, but also affect smaller,
disconnected glaciers such as Griaz and Bourgeat. The latter
is true as long as the glaciers are within the correlation range
of 11.2 km.

VII. CONCLUSIONS

In this study, we reviewed the literature on the accuracy and precision of DEMs. On the basis of the raised considerations regarding variable vertical precision and correlated noises, we proposed a non-stationary spatial framework for DEM uncertainty analysis. This framework allows to perform inference on a single difference between a DEM and independent elevation data on stable terrain, and to distinguish elevation biases from elevation variance. We developed robust methods to estimate and model both elevation heteroscedasticity and spatial correlation of elevation errors. We then validated that stable terrain is a reliable error proxy for other terrain types using pairs of DEMs derived from nearly simultaneous acquisitions for the Mont-Blanc massif and the Northern Patagonian Icefield.

We illustrated the impact of our methods when propagating
uncertainties to pixel-scale and spatial derivatives of elevation.
For the pixel-scale terrain slope, uncertainties are underesti-
mated by up to a factor 10 in rough and steep topography
when omitting elevation heteroscedasticity. For glacier volume
changes, the uncertainty of the volume change of a glacier
of 10 km² is underestimated by a factor of 25 when omitting
correlations with ranges of 3.9 and 11.2 km, despite their
small cumulative magnitude of only 7% of the variance. This
underestimation of long-range spatial correlation affects many
studies relying on instruments plagued by noise, such as the
widely used DEMs from SRTM and ASTER.

We provide an implementation of our methods in the
Python package xDEM [157], which includes, in particular,
DEM alignment, correction, and uncertainty analysis. Spatial
statistics have long been used for uncertainty analysis, yet
often suffered from a lack of accessibility [126]. The wider
application of such analysis was still deemed as "unrealized"
a decade ago [133], possibly also due to the scarcity of
open-source and documented tools for spatial statistics. By
providing our methodological tools within the frame of a
package embedded in high-level programming languages that
efficiently pairs with remote sensing analysis, we hope to
foster a consistent, reproducible and accessible uncertainty
analysis of DEMs.

We highlight the genericity of our spatial framework for
uncertainty analysis and of our estimation methods for dense
and outlier-prone grid data. Our framework holds the potential
to be extended to other geospatial data. Gridded surface
displacement, for instance, profit from the same error proxy of
stable terrain and are increasingly used in a variety of appli-
cations. To describe the precision of such spatially structured
data, we advocate for the use of additional metrics. These
metrics should describe potential heteroscedasticity and spatial
correlation of errors, reported, for example, in a tabular man-
ner — parameters of variogram models; discrete categories
of heteroscedasticity. Ultimately, the adoption of such new
metrics is critical to progress towards a realistic description
of error structure in geospatial data, and a robust propagation
of uncertainties in Earth system science assessments.

DATA AND CODE AVAILABILITY STATEMENT

The code developed for the processing and analysis of all data, and
the generation of figures and tables in this article is available at https:
//github.com/rhugonnet/dem_error_study, with routines implemented
in the DEM analysis Python package xDEM available at https:
//github.com/GlacioHack/xdem with supporting documentation at
https://xdem.readthedocs.io, and the spatial statistics Python package
SciKit-GStat available at https://github.com/mmaelicke/scikit-gstat
with supporting documentation at https://scikit-gstat.readthedocs.io.

AUTHOR CONTRIBUTIONS

R.H., F.B. and E.B. designed the study. F.B. performed an early
analysis of spatial correlation for the Mont-Blanc case study. R.H.
developed and tested non-stationary spatial statistics methods for
DEMs with inputs from N.E., A.D., F.B. and E.B. R.H. implemented
the methods in the Python package xDEM with main inputs from
A.D. and E.M. All authors interpreted the results. R.H. performed
the literature review and led the writing of the paper, and all other
co-authors contributed.

ACKNOWLEDGMENTS

We thank Dorde Masovic for realizing Fig. 2 with vector graphics.
We thank Alexandre Wadoux for insightful exchanges on spatial
statistics and Simon Gascoin for discussions on along-track undu-
lations in Pléiades DEMs. We also thank Eyjólfur Magnússon and
Joaquín Belart for providing insights into Gaussian simulations for
DEMs. We thank Jonathan Bamber, César Deschamps-Berger and
Matthias Huss for comments on early versions of the manuscript.
SPOT6/7 data were obtained from GEOSUD (ANR-10-EQPX-20,
programme ‘Investissements d’Avenir’). INRAE is a member of
Labex OSUG. R.H. acknowledges a fellowship from the University
of Toulouse. R.H., A.D., E.M and D.F. acknowledge funding from the
Swiss National Science Foundation, grant number 184634. E.B. ac-
knowledges support from the French Space Agency (CNES) through
ISIS and TOSCA programmes.

REFERENCES

[1] H. Taud, J.-F. Parrot, and R. Alvarez, "DEM generation by contour line
dilation," *Comput. Geosci.*, vol. 25, no. 7, pp. 775-783, Aug. 1999.

[2] Q. Weng, "Quantifying uncertainty of digital elevation models derived from topographic maps," in *Advances in Spatial Data Handling*. Springer Berlin Heidelberg, 2002, pp. 403–418.

[3] P. A. Rosen, S. Hensley, I. R. Joughin, F. K. Li, S. N. Madsen, E. R-
driguez, and R. M. Goldstein, "Synthetic aperture radar interferometry,"
*Proc. IEEE*, vol. 88, no. 3, pp. 333–382, Mar. 2000.

[4] R. Bürgmann, P. A. Rosen, and E. J. Fielding, “Synthetic aperture radar interferometry to measure earth’s surface topography and its deformation,” *Annu. Rev. Earth Planet. Sci.*, vol. 28, no. 1, pp. 169–209, May 2000.
[5] A. S. Walker, "Analogue, analytical and digital photogrammetric workstations: Practical investigations of performance," *The Photogrammetric Record*, 1995.

[6] E. M. Mikhail, J. S. Bethel, and J. C. McGlone, *Introduction to Modern Photogrammetry*. Nashville, TN: John Wiley & Sons, Mar. 2001, vol. 19.

[7] E. P. Baltsavias, "Airborne laser scanning: basic relations and formulas," *ISPRS J. Photogramm. Remote Sens.*, vol. 54, no. 2-3, pp. 199-214, Jul. 1999.

[8] R. O. Dubayah and J. B. Drake, "Lidar remote sensing for forestry," *J. For.*, vol. 98, no. 6, pp. 44-46, Jun. 2000.

[9] S. Bhushan, D. Shean, O. Alexandrov, and S. Henderson, "Automated digital elevation model (DEM) generation from very-high-resolution planet SkySat triplet stereo and video imagery," *ISPRS J. Photogramm. Remote Sens.*, vol. 173, pp. 151-165, Mar. 2021.

[10] J. Biggs and T. J. Wright, "How satellite InSAR has grown from opportunistic science to routine monitoring over the last decade," *Nat. Commun.*, vol. 11, no. 1, p. 3863, Aug. 2020.

[11] J. Rosen, "Shifting ground," *Science*, vol. 371, no. 6532, pp. 876-880, Feb. 2021.

[12] N. Galiatsatos, D. N. M. Donoghue, and G. Philip, "High resolution elevation data derived from stereoscopic CORONA imagery with minimal ground control," *Photogrammetric Engineering & Remote Sensing*, vol. 73, no. 9, pp. 1093-1106, 2007.

[13] A. Dehecq, A. S. Gardner, O. Alexandrov, S. McMichael, R. Hugonnet, D. Shean, and M. Marty, "Automated processing of declassified KH-9 hexagon satellite images for global elevation change analysis since the 1970s," *Front Earth Sci.*, vol. 8, p. 516, 2020.

[14] E. C. Geyman, W. J J van Pelt, A. C. Maloof, H. F. Aas, and J. Kohler, "Historical glacier change on svalbard predicts doubling of mass loss by 2100," *Nature*, vol. 601, no. 7893, pp. 374-379, Jan. 2022.

[15] G. LeFavour and D. Alsdorf, "Water slope and discharge in the amazon river estimated using the shuttle radar topography mission digital elevation model," *Geophys. Res. Lett.*, vol. 32, no. 17, Sep. 2005.

[16] J. Oksanen and T. Sarjakoski, "Error propagation analysis of DEM-based drainage basin delineation," *Int. J. Remote Sens.*, vol. 26, no. 14, pp. 3085-3102, Jul. 2005.

[17] L. Hawker, P. Bates, J. Neal, and J. Rougier, "Perspectives on digital elevation model (DEM) simulation for flood modeling in the absence of a High-Accuracy open access global DEM," *Front Earth Sci.*, vol. 6, p. 233, 2018.

[18] B. Székely and D. Karátson, "DEM-based morphometry as a tool for reconstructing primary volcanic landforms: examples from the börzszöny mountains, hungary," *Geomorphology*, vol. 63, no. 1, pp. 25-37, Nov. 2004.

[19] T. Bolch, U. Kamp, J. Olsenholler, and M. Oluić, "Using ASTER and SRTM DEMs for studying geomorphology and glaciation in high mountain areas," in *New Strategies for European Remote Sensing*, M. Oluić, Ed. Rotterdam: Millpress, 2005, pp. 119-127.

[20] B. Žabota, B. Repe, and M. Kobal, "Influence of digital elevation model resolution on rockfall modelling," *Geomorphology*, vol. 328, pp. 183-195, Mar. 2019.

[21] A. Mueting, B. Bookhagen, and M. R. Strecker, "Identification of debris-flow channels using high-resolution topographic data: A case study in the quebrada del toro, NW argentina," *J. Geophys. Res. Earth Surf.*, vol. 126, no. 12, Dec. 2021.

[22] A. Ganas, G. Papadopoulos, and S. B. Pavlides, "The 7 september 1999 athens 5.9 ms earthquake: Remote sensing and digital elevation model inputs towards identifying the seismic fault," *Int. J. Remote Sens.*, vol. 22, no. 1, pp. 191-196, Jan. 2001.

[23] S. Vassilopoulou, L. Hurni, V. Dietrich, E. Baltsavias, M. Pateraki, E. Lagios, and I. Parcharidis, "Orthophoto generation using IKONOS imagery and high-resolution DEM: a case study on volcanic hazard monitoring of nisyros island (greece)," *ISPRS J. Photogramm. Remote Sens.*, vol. 57, no. 1, pp. 24-38, Nov. 2002.

[24] P. Grosse, B. van Wyk de Vries, P. A. Euillades, M. Kervyn, and I. A. Petrinovic, "Systematic morphometric characterization of volcanic edifices using digital elevation models," *Geomorphology*, vol. 136, no. 1, pp. 114-131, Jan. 2012.

[25] K. K. Singh and A. Singh, "Detection of 2011 sikkim earthquake-induced landslides using neuro-fuzzy classifier and digital elevation model," *Nat. Hazards*, vol. 83, no. 2, pp. 1027-1044, Sep. 2016.

[26] Y. Bühler, S. Kumar, J. Veitinger, M. Christen, A. Stoffel, and Others, "Automated identification of potential snow avalanche release areas based on digital elevation models," *Nat. Hazards Earth Syst. Sci.*, vol. 13, no. 5, pp. 1321-1335, 2013.

[27] A. E. Akay, H. Oğuz, I. R. Karas, and K. Aruga, "Using LiDAR technology in forestry activities," *Environ. Monit. Assess.*, vol. 151, no. 1-4, pp. 117-125, Apr. 2009.

[28] A. Räsänen, M. Kuitunen, E. Tomppo, and A. Lensu, "Coupling high-resolution satellite imagery with ALS-based canopy height model and digital elevation model in object-based boreal forest habitat type classification," *ISPRS J. Photogramm. Remote Sens.*, vol. 94, pp. 169-182, Aug. 2014.

[29] J. Griffin, H. Latief, W. Kongko, S. Harig, N. Horspool, R. Hanung, A. Rojali, N. Maher, A. Fuchs, J. Hossen, S. Upi, D. Edi, N. Rakowsky, and P. Cummins, "An evaluation of onshore digital elevation models for modeling tsunami inundation zones," *Front Earth Sci.*, vol. 3, p. 32, 2015.

[30] S. Kulp and B. H. Strauss, "Global DEM errors underpredict coastal vulnerability to sea level rise and flooding," *Front Earth Sci.*, vol. 4, p. 36, 2016.

[31] A. Hooijer and R. Vernimmen, "Global LiDAR land elevation data reveal greatest sea-level rise vulnerability in the tropics," *Nat. Commun.*, vol. 12, no. 1, p. 3592, Jun. 2021.

[32] M. K. Barker, E. Mazarico, G. A. Neumann, M. T. Zuber, J. Haruyama, and D. E. Smith, "A new lunar digital elevation model from the lunar orbiter laser altimeter and SELENE terrain camera," *Icarus*, vol. 273, pp. 346-355, Jul. 2016.

[33] Y. Wang and B. Wu, "Active machine learning approach for crater detection from planetary imagery and digital elevation models," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 8, pp. 5777-5789, Aug. 2019.

[34] Y.-C. Hsieh, Y.-C. Chan, and J.-C. Hu, "Digital elevation model differencing and error estimation from multiple sources: A case study from the meiyuan shan landslide in taiwan," *Remote Sensing*, vol. 8, no. 3, p. 199, Feb. 2016.

[35] P. Lacroix, "Landslides triggered by the gorkha earthquake in the langtang valley, volumes and initiation processes the 2015 gorkha, nepal, earthquake and himalayan studies: First results 4. seismology," *Earth Planets Space*, vol. 68, no. 1, pp. 1-10, 2016.

[36] D. H. Shugar, M. Jacquemart, D. Shean, S. Bhushan, K. Upadhyay, A. Sattar, W. Schwanghart, S. McBride, M. V. W. de Vries, M. Mergili, A. Emmer, C. Deschamps-Berger, M. McDonnell, R. Bhambri, S. Allen, E. Berthier, J. L. Carrivick, J. J. Clague, M. Dokukin, S. A. Dunning, H. Frey, S. Gascoin, U. K. Haritashya, C. Huggel, A. Kääb, J. S. Kargel, J. L. Kavanaugh, P. Lacroix, D. Petley, S. Rupper, M. F. Azam, S. J. Cook, A. P. Dimri, M. Eriksson, D. Farinotti, J. Fiddes, K. R. Gnyawali, S. Harrison, M. Jha, M. Koppes, A. Kumar, S. Leinss, U. Majeed, S. Mal, A. Muhuri, J. Noetzli, F. Paul, I. Rashid, K. Sain, J. Steiner, F. Ugalde, C. S. Watson, and M. J. Westoby, "A massive rock and ice avalanche caused the 2021 disaster at chamoli, indian himalaya," *Science*, Jun. 2021.

[37] M. Nolan, C. Larsen, and M. Sturm, "Mapping snow depth from manned aircraft on landscape scales at centimeter resolution using structure-from-motion photogrammetry," *The Cryosphere*, vol. 9, no. 4, pp. 1445-1463, Aug. 2015.

[38] R. Marti, S. Gascoin, E. Berthier, M. d. Pinel, T. Houet, and D. Laffly, "Mapping snow depth in open alpine terrain from stereo satellite imagery," *The Cryosphere*, vol. 10, no. 4, pp. 1361-1380, 2016.

[39] C. Deschamps-Berger, S. Gascoin, E. Berthier, J. Deems, E. Gutmann, A. Dehecq, D. Shean, and M. Dumont, "Snow depth mapping from stereo satellite imagery in mountainous terrain: evaluation using airborne laser-scanning data," *The Cryosphere*, vol. 14, no. 9, pp. 2925-2940, Sep. 2020.

[40] M. Bagnardi, P. J. González, and A. Hooper, "High-resolution digital elevation model from tri-stereo pleiades-1 satellite imagery for lava flow volume estimates at fogo volcano," *Geophys. Res. Lett.*, vol. 43, no. 12, pp. 6267-6275, Jun. 2016.

[41] G. B. M. Pedersen, J. M. C. Belart, E. Magnússon, O. K. Vilmundardóttir, F. Kizel, F. S. Sigurmundsson, G. Gísladóttir, and J. A. Benediktsson, "Hekla volcano, iceland, in the 20th century: Lava volumes, production rates, and effusion rates," *Geophys. Res. Lett.*, vol. 45, no. 4, pp. 1805-1813, Feb. 2018.

[42] C. Ticehurst, S. Phinn, and A. Held, "Using multitemporal digital elevation model data for detecting canopy gaps in tropical forests due to cyclone damage: An initial assessment," *Austral Ecol.*, vol. 32, no. 1, pp. 59-69, Feb. 2007.

[43] L. Magruder, H. Leigh, and A. Neuenschwander, "Evaluation of terrain and canopy height products in central african tropical forests," *Int. J. Remote Sens.*, vol. 37, no. 22, pp. 5365-5387, Nov. 2016.

[44] K. Gdulová, J. Marešová, V. Barták, M. Szostak, J. Červenka, and V. Moudrý, "Use of TanDEM-X and SRTM-C data for detection of
deforestation caused by bark beetle in central european mountains," *Remote Sensing*, vol. 13, no. 15, p. 3042, Aug. 2021.

[45] D. Felikson, T. C. Bartholomaus, G. A. Catania, N. J. Korsgaard, K. H. Kjær, M. Morlighem, B. Noël, M. van den Broeke, L. A. Stearns, E. L. Shroyer, D. A. Sutherland, and J. D. Nash, "Inland thinning on the greenland ice sheet controlled by outlet glacier geometry," *Nat. Geosci.*, vol. 10, no. 5, pp. 366-369, May 2017.

[46] D. E. Shean, I. R. Joughin, P. Dutrieux, B. E. Smith, and E. Berthier, "Ice shelf basal melt rates from a high-resolution digital elevation model (DEM) record for pine island glacier, antarctica," *The Cryosphere*, vol. 13, no. 10, pp. 2633-2656, 2019.

[47] R. Hugonnet, R. McNabb, E. Berthier, B. Menounos, C. Nuth, L. Girod, D. Farinotti, M. Huss, I. Dussaillant, F. Brun, and A. Kääb, "Accelerated global glacier mass loss in the early twenty-first century," *Nature*, vol. 592, no. 7856, pp. 726-731, Apr. 2021.

[48] J. Wood, "The geomorphological characterisation of digital elevation models," Ph.D. dissertation, University of Leicester (United Kingdom), Ann Arbor, United States, 1996.

[49] B. H. Carlisle, "Digital elevation model quality and uncertainty in DEM-based spatial modelling," Ph.D. dissertation, University of Greenwich, London, UK, 2002.

[50] J. Oksanen, "Digital elevation model error in terrain analysis," Ph.D. dissertation, University of Helsinki, Faculty of Science, Nov. 2006.

[51] C. Papasaika-Hanusch, "Fusion of digital elevation models," Ph.D. dissertation, ETH Zurich, 2012.

[52] D. F. Maune, *Digital elevation model technologies and applications: the DEM users manual*. Asprs Publications, 2007.

[53] P. F. Fisher and N. J. Tate, "Causes and consequences of error in digital elevation models," *Progress in Physical Geography: Earth and Environment*, vol. 30, no. 4, pp. 467-489, Aug. 2006.

[54] S. P. Wechsler, "Uncertainties associated with digital elevation models for hydrologic applications: a review," *Hydrol. Earth Syst. Sci.*, vol. 11, no. 4, pp. 1481-1500, Aug. 2007.

[55] S. E. Lakshmi and K. Yarrakula, "Review and critical analysis on digital elevation models," *Geofizika*, vol. 35, no. 2, pp. 129-157, 2018.

[56] L. Polidori and M. El Hage, "Digital elevation model quality assessment methods: A critical review," *Remote Sensing*, vol. 12, no. 21, p. 3522, Oct. 2020.

[57] F. Hebeler and R. S. Purves, "The influence of elevation uncertainty on derivation of topographic indices," *Geomorphology*, vol. 111, no. 1, pp. 4-16, Oct. 2009.

[58] T. G. Farr, P. A. Rosen, E. Caro, R. Crippen, R. Duren, S. Hensley, M. Kobrick, M. Paller, E. Rodriguez, L. Roth, D. Seal, S. Shaffer, J. Shimada, J. Umland, M. Werner, M. Oskin, D. Burbank, and D. Alsdorf, "The shuttle radar topography mission," *Rev. Geophys.*, vol. 45, no. 2, p. 1485, May 2007.

[59] P. Rizzoli, M. Martone, C. Gonzalez, C. Wecklich, D. B. Tridon, B. Bräutigam, M. Bachmann, D. Schulze, T. Fritz, M. Huber, B. Wesel, G. Krieger, M. Zink, and A. Moreira, "Generation and performance assessment of the global TanDEM-X digital elevation model," *ISPRS J. Photogramm. Remote Sens.*, vol. 132, pp. 119-139, 2017.

[60] G. Amatulli, D. McInerney, T. Sethi, P. Strobl, and S. Domisch, "Geomorpho90m, empirical evaluation and accuracy assessment of global high-resolution geomorphometric layers," *Sci Data*, vol. 7, no. 1, p. 162, May 2020.

[61] E. Uuemaa, S. Ahi, B. Montibeller, M. Muru, and A. Kmoch, "Vertical accuracy of freely available global digital elevation models (ASTER, AW3D30, MERIT, TanDEM-X, SRTM, and NASADEM)," *Remote Sensing*, vol. 12, no. 21, p. 3482, Oct. 2020.

[62] L. A. Eberhard, P. Sirguey, A. Miller, M. Marty, K. Schindler, A. Stofel, and Y. Bühler, "Intercomparison of photogrammetric platforms for spatially continuous snow depth mapping," *The Cryosphere*, vol. 15, no. 1, pp. 69-94, Jan. 2021.

[63] L. Magruder, A. Neuenschwander, and B. Klotz, "Digital terrain model elevation corrections using space-based imagery and ICESat-2 laser altimetry," *Remote Sens. Environ.*, vol. 264, p. 112621, Oct. 2021.

[64] B. H. Carlisle, "Modelling the spatial distribution of DEM error," *Trans. GIS*, vol. 9, no. 4, pp. 521-540, Oct. 2005.

[65] C. Rolstad, T. Haug, and B. Denby, "Spatially integrated geodetic glacier mass balance and its uncertainty based on geostatistical analysis: Application to the western svartsen ice cap, norway," *J. Glaciol.*, vol. 55, no. 192, pp. 666-680, 2009.

[66] L. Girod, C. Nuth, A. Kääb, R. McNabb, and O. Galland, "MMASTER: Improved ASTER DEMs for elevation change monitoring," *Remote Sensing*, vol. 9, no. 7, p. 704, Jul. 2017.

[67] ISO (International Organisation for Standardisation), *ISO 5725-1: 1994: Accuracy (Trueness and Precision) of Measurement Methods and Results-Part 1: General Principles and Definitions*. Geneva: ISO, 1994.

[68] L. Guan, H. Pan, S. Zou, J. Hu, X. Zhu, and P. Zhou, "The impact of horizontal errors on the accuracy of freely available digital elevation models (DEMs)," *Int. J. Remote Sens.*, vol. 41, no. 19, pp. 7383-7399, Oct. 2020.

[69] S. Mukherjee, P. K. Joshi, S. Mukherjee, A. Ghosh, R. D. Garg, and A. Mukhopadhyay, "Evaluation of vertical accuracy of open source digital elevation model (DEM)," *Int. J. Appl. Earth Obs. Geoinf.*, vol. 21, pp. 205-217, Apr. 2013.

[70] S. A. Kulp and B. H. Strauss, "New elevation data triple estimates of global vulnerability to sea-level rise and coastal flooding," *Nat. Commun.*, vol. 10, no. 1, p. 4844, Oct. 2019.

[71] A. Gruber, B. Wessel, M. Huber, and A. Roth, "Operational TanDEM-X DEM calibration and first validation results," *ISPRS J. Photogramm. Remote Sens.*, vol. 73, pp. 39-49, Sep. 2012.

[72] A. Dehecq, R. Millan, E. Berthier, N. Gourmelen, E. Trouvé, and V. Vionnet, "Elevation changes inferred from TanDEM-X data over the Mont-Blanc area: Impact of the X-Band interferometric bias," *IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing*, vol. 9, no. 8, pp. 3870-3882, Aug. 2016.

[73] P. J. Besl and N. D. McKay, "A method for registration of 3-D shapes," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 14, no. 2, pp. 239-256, Feb. 1992.

[74] A. Myronenko and X. Song, "Point set registration: coherent point drift," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 32, no. 12, pp. 2262-2275, Dec. 2010.

[75] Q.-Y. Zhou, J. Park, and V. Koltun, "Fast global registration," in *Computer Vision – ECCV 2016*. Springer International Publishing, 2016, pp. 766-782.

[76] T. Zhang and M. Cen, "Robust DEM co-registration method for terrain changes assessment using least trimmed squares estimator," *Adv. Space Res.*, vol. 41, no. 11, pp. 1827-1835, Jan. 2008.

[77] C. Nuth and Kääb, "Co-registration and bias corrections of satellite elevation data sets for quantifying glacier thickness change," *The Cryosphere*, vol. 5, no. 1, pp. 271-290, 2011.

[78] M.-J. Noh and I. M. Howat, "Automated coregistration of repeat digital elevation models for surface elevation change measurement using geometric constraints," *IEEE Trans. Geosci. Remote Sens.*, vol. 52, no. 4, pp. 2247-2260, Apr. 2014.

[79] J. B. Abshire, X. Sun, H. Riris, J. M. Sirota, J. F. McGarry, S. Palm, D. Yi, and P. Liiva, "Geoscience laser altimeter system (GLAS) on the ICESat mission: On-orbit measurement performance," *Geophys. Res. Lett.*, vol. 32, no. 21, 2005.

[80] W. Abdalati, H. J. Zwally, R. Bindschadler, B. Csatho, S. L. Farrell, H. A. Fricker, D. Harding, R. Kwok, M. Lefsky, T. Markus, A. Marshak, T. Neumann, S. Palm, B. Schutz, B. Smith, J. Spinhirne, and C. Webb, "The ICESat-2 laser altimetry mission," *Proc. IEEE*, vol. 98, no. 5, pp. 735-751, May 2010.

[81] J. Vaze, J. Teng, and G. Spencer, "Impact of DEM accuracy and resolution on topographic indices," *Environmental Modelling & Software*, vol. 25, no. 10, pp. 1086-1098, Oct. 2010.

[82] J. Gardelle, E. Berthier, and Y. Arnaud, "Impact of resolution and radar penetration on glacier elevation changes computed from DEM differencing," *J. Glaciol.*, vol. 58, no. 208, pp. 419-422, 2012.

[83] J. Dall, S. N. Madsen, K. Keller, and R. Forsberg, "Topography and penetration of the greenland ice sheet measured with airborne SAR interferometry," *Geophys. Res. Lett.*, vol. 28, no. 9, pp. 1703-1706, May 2001.

[84] J. Li, Z.-W. Li, J. Hu, L.-X. Wu, X. Li, L. Guo, Z. Liu, Z.-L. Miao, W. Wang, and J.-L. Chen, "Investigating the bias of TanDEM-X digital elevation models of glaciers on the tibetan plateau: impacting factors and potential effects on geodetic mass-balance measurements," *J. Glaciol.*, vol. 67, no. 264, pp. 613-626, Aug. 2021.

[85] F. Ackermann, "Digital image correlation: Performance and potential application in photogrammetry," *Photogramm. Rec.*, vol. 11, no. 64, pp. 429-439, Oct. 1984.

[86] R. M. Goldstein, H. A. Zebker, and C. L. Werner, "Satellite radar interferometry: Two-dimensional phase unwrapping," *Radio Sci.*, vol. 23, no. 4, pp. 713-720, Jul. 1988.

[87] H. Lee and J. G. Liu, "Analysis of topographic decorrelation in SAR interferometry using ratio coherence imagery," *IEEE Trans. Geosci. Remote Sens.*, vol. 39, no. 2, pp. 223-232, Feb. 2001.
[88] A. Wehr and U. Lohr, "Airborne laser scanning—an introduction and overview," *ISPRS J. Photogramm. Remote Sens.*, vol. 54, no. 2-3, pp. 68–82, Jul. 1999.

[89] F. J. Aguilar and J. P. Mills, "Accuracy assessment of lidar-derived digital elevation models," *Photogramm. Rec.*, vol. 23, no. 122, pp. 148–169, Jun. 2008.

[90] F. J. Aguilar, J. P. Mills, J. Delgado, M. A. Aguilar, J. G. Negreiros, and J. L. Pérez, "Modelling vertical error in LiDAR-derived digital elevation models," *ISPRS J. Photogramm. Remote Sens.*, vol. 65, no. 1, pp. 103–110, Jan. 2010.

[91] P. V. Bolstad and T. Stowe, "An evaluation of DEM accuracy: elevation, slope, and aspect," *Photogrammetric Engineering & Remote Sensing*, vol. 60, no. 11, pp. 1327–1332, 1994.

[92] Y. Gyasi-Agyei, G. Willgoose, and F. P. De Troch, "Effects of vertical resolution and map scale of digital elevation models on geomorphological parameters used in hydrology," *Hydrol. Process.*, vol. 9, no. 3-4, pp. 363–382, Apr. 1995.

[93] K. W. Holmes, O. A. Chadwick, and P. C. Kyriakidis, "Error in a USGS 30-meter digital elevation model and its impact on terrain modeling," *J. Hydrol.*, vol. 233, no. 1, pp. 154–173, Jun. 2000.

[94] T. Toutin, "Three-dimensional topographic mapping with ASTER stereo data in rugged topography," *IEEE Trans. Geosci. Remote Sens.*, vol. 40, no. 10, pp. 2241–2247, 2002.

[95] J. M. Wheaton, J. Brasington, S. E. Darby, and D. A. Sear, "Accounting for uncertainty in DEMs from repeat topographic surveys: improved sediment budgets," *Earth Surf. Processes Landforms*, 2009.

[96] D. J. Milan, G. L. Heritage, A. R. G. Large, and I. C. Fuller, "Filtering spatial error from DEMs: Implications for morphological change estimation," *Geomorphology*, vol. 125, no. 1, pp. 160–171, Jan. 2011.

[97] A. Bouillon, M. Bernard, P. Gigord, A. Orsoni, V. Rudowski, and A. Baudoin, "SPOT 5 HRS geometric performances: Using block adjustment as a key issue to improve quality of DEM generation," *ISPRS J. Photogramm. Remote Sens.*, vol. 60, no. 3, pp. 134–146, May 2006.

[98] I. Dussaillant, E. Berthier, and F. Brun, "Geodetic mass balance of the northern patagonian icefield from 2000 to 2012 using two independent methods," *Front Earth Sci.*, vol. 6, no. February, p. 8, 2018.

[99] R. Perko, H. Raggam, and P. M. Roth, "Mapping with Pléiades—End-to-End workflow," *Remote Sensing*, vol. 11, no. 17, p. 2052, Sep. 2019.

[100] G. Falorni, "Analysis and characterization of the vertical accuracy of digital elevation models from the shuttle radar topography mission," *J. Geophys. Res.*, vol. 110, no. F2, 2005.

[101] M. Simard, M. Neumann, and S. Buckley, "Validation of the new SRTM digital elevation model (NASADEM) with ICESAT/GLAS over the united states," in *2016 IEEE International Geoscience and Remote Sensing Symposium (IGARSS)*. ieexplore.ieee.org, Jul. 2016, pp. 3227–3229.

[102] J. M. C. Belart, E. Berthier, E. Magnússon, L. S. Anderson, F. Pálsson, T. Thorsteinsson, I. M. Howat, G. Aðalgeirsdóttir, T. Jóhannesson, and A. H. Jarosch, "Winter mass balance of drangajökull ice cap (NW Iceland) derived from satellite sub-meter stereo images," *The Cryosphere*, vol. 11, no. 3, pp. 1501–1517, Jun. 2017.

[103] L. Girod, N. I. Nielsen, F. Couderette, C. Nuth, and A. Kääb, "Precise DEM extraction from svalbard using 1936 high oblique imagery," *Geosci. Instrum. Methods Data Syst.*, vol. 7, no. 4, pp. 277–288, Oct. 2018.

[104] W. Gheyle, J. Bourgeois, R. Goossens, and K. Jacobsen, "Scan problems in digital CORONA satellite images from USGS archives," *Photogrammetric Engineering & Remote Sensing*, vol. 77, no. 12, pp. 1257–1264, 2011.

[105] J. Li, Z.-W. Li, J.-J. Zhu, X. Li, B. Xu, Q.-J. Wang, C.-L. Huang, and J. Hu, "Early 21st century glacier thickness changes in the central tien shan," *Remote Sens. Environ.*, vol. 192, pp. 12–29, Apr. 2017.

[106] M. S. Hamid and M. Safy, "InSAR image denoising filter for accurate DEM generation," in *2020 12th International Conference on Electrical Engineering (ICEENG)*, Jul. 2020, pp. 306–310.

[107] L. Bopche and P. P. Rege, "Use of noise reduction filters on stereo images for improving the accuracy and quality of the digital elevation model," *JARS*, vol. 15, no. 1, p. 014508, Feb. 2021.

[108] T. Ai and J. Li, "A DEM generalization by minor valley branch detection and grid filling," *ISPRS J. Photogramm. Remote Sens.*, vol. 65, no. 2, pp. 198–207, Mar. 2010.

[109] T. Smith, A. Rheinwalt, and B. Bookhagen, "Determining the optimal grid resolution for topographic analysis on an airborne lidar dataset," *Earth Surf. Dynam.*, vol. 7, no. 2, pp. 475–489, May 2019.

[110] J. D. Wood and P. F. Fisher, "Assessing interpolation accuracy in elevation models," *IEEE Comput. Graph. Appl.*, vol. 13, no. 2, pp. 48–56, Mar. 1993.

[111] J. Gao, "Resolution and accuracy of terrain representation by grid DEMs at a micro-scale," *Int. J. Geogr. Inf. Sci.*, vol. 11, no. 2, pp. 199–212, Mar. 1997.

[112] A. G. Journel and C. J. Huijbregts, *Mining geostatistics*. Academic press London, 1978, vol. 600.

[113] N. A. C. Cressie, *Statistics for spatial data*. New York: Wiley, 1993, vol. 4.

[114] C. R. Ehlschlaeger and A. Shortridge, "Modeling elevation uncertainty in geographical analyses," in *Proceedings of the International Symposium on Spatial Data Handling, Delft, Netherlands B*, vol. 9, 1996, pp. 15–9B.

[115] P. C. Kyriakidis, A. M. Shortridge, and M. F. Goodchild, "Geostatistics for conflation and accuracy assessment of digital elevation models," *Int. J. Geogr. Inf. Sci.*, vol. 13, no. 7, pp. 677–707, Oct. 1999.

[116] S. P. Wechsler, "Digital elevation model (DEM) uncertainty: evaluation and effect on topographic parameters," in *ESRI User Conference. cnc.purdue.edu, 1999*, pp. 1081–1090.

[117] P. Goovaerts, *Geostatistics for Natural Resources Evaluation*. Oxford University Press, 1997.

[118] G. B. M. Heuvelink, *Error propagation in environmental modelling with GIS*. CRC press, 1998.

[119] A. Shortridge, "Characterizing uncertainty in digital elevation models," in *Spatial Uncertainty in Ecology: Implications for Remote Sensing and GIS Applications*, C. T. Hunsaker, M. F. Goodchild, M. A. Friedl, and T. J. Case, Eds. New York, NY: Springer New York, 2001, pp. 238–257.

[120] G. Wang, G. Z. Gertner, S. Fang, and A. B. Anderson, "A methodology for spatial uncertainty analysis of remote sensing and GIS products," *Photogrammetric Engineering & Remote Sensing*, vol. 71, no. 12, pp. 1423–1432, 2005.

[121] G. B. M. Heuvelink, P. A. Burrough, and A. Stein, "Propagation of errors in spatial modelling with GIS," *International Journal of Geographical Information Systems*, vol. 3, no. 4, pp. 303–322, Oct. 1989.

[122] P. F. Fisher, "First experiments in viewshed uncertainty: the accuracy of the viewshed area," *Photogramm. Eng. Remote Sens.*, vol. 57, no. 10, pp. 1321–1327, 1991.

[123] S. Openshaw, M. Charlton, and S. Carver, "Error propagation: a monte carlo simulation," in *Handling geographical information: methodology and potential applications*. Longman Harlow, Essex, 1991, pp. 78–101.

[124] F. Heße, V. Prykhodko, S. Schlüter, and S. Attinger, "Generating random fields with a truncated power-law variogram: A comparison of several numerical methods," *Environmental Modelling & Software*, vol. 55, pp. 32–48, May 2014.

[125] F. Canters, W. D. Genst, and H. Dufourmont, "Assessing effects of input uncertainty in structural landscape classification," *Int. J. Geogr. Inf. Sci.*, vol. 16, no. 2, pp. 129–149, Mar. 2002.

[126] S. P. Wechsler and C. N. Kroll, "Quantifying DEM uncertainty and its effect on topographic parameters," *Photogrammetric Engineering & Remote Sensing*, vol. 72, no. 9, pp. 1081–1090, 2006.

[127] E. Magnússon, J. Muñoz-Cobo Belart, F. Pálsson, H. Ágústsson, and P. Crochet, "Geodetic mass balance record with rigorous uncertainty estimates deduced from aerial photographs and lidar data - case study from drangajökull ice cap, NW Iceland," *The Cryosphere*, vol. 10, no. 1, pp. 159–177, Jan. 2016.

[128] A. Bachmann and B. Allgöwer, "Uncertainty propagation in wildland fire behaviour modelling," *Int. J. Geogr. Inf. Sci.*, vol. 16, no. 2, pp. 115–127, Mar. 2002.

[129] M. Albani, B. Klinkenberg, D. W. Andison, and J. P. Kimmins, "The choice of window size in approximating topographic surfaces from digital elevation models," *Int. J. Geogr. Inf. Sci.*, vol. 18, no. 6, pp. 577–593, Sep. 2004.

[130] C. Li, Q. Wang, W. Shi, and S. Zhao, "Uncertainty modelling and analysis of volume calculations based on a regular grid digital elevation model (DEM)," *Comput. Geosci.*, vol. 114, pp. 117–129, May 2018.

[131] C. Li, S. Zhao, Q. Wang, and W. Shi, "Uncertainty modeling and analysis of surface area calculation based on a regular grid digital elevation model (DEM)," *Int. J. Geogr. Inf. Sci.*, vol. 32, no. 9, pp. 1837–1859, Sep. 2018.

[132] S. W. Anderson, "Uncertainty in quantitative analyses of topographic change: error propagation and the role of thresholding," *Earth Surf. Processes Landforms*, vol. 44, no. 5, pp. 1015–1033, 2019.
[133] G. B. M. Heuvelink, "Analysing uncertainty propagation in GIS: Why is it not that simple?" in *Uncertainty in Remote Sensing and GIS*. Chichester, UK: John Wiley & Sons, Ltd, Jul. 2006, pp. 155-165.

[134] R. A. Beyer, O. Alexandrov, and S. McMichael, "The ames stereo pipeline: NASA's open source software for deriving and processing terrain data," *Earth and Space Science*, vol. 5, no. 9, pp. 537-548, Sep. 2018.

[135] NASA/METI/AIST/Japan Spacesystems, and U.S./Japan ASTER Science Team, "ASTER level 1A data set - reconstructed, unprocessed instrument data," 2001.

[136] E. Rupnik, M. Daakir, and M. Pierrot Deseilligny, "MicMac - a free, open-source solution for photogrammetry," *Open Geospatial Data, Software and Standards*, vol. 2, no. 1, p. 14, 2017.

[137] NASA JPL, "NASADEM SRTM-only height and height precision mosaic global 1 arc second," 2020.

[138] M.-J. Noh and I. M. Howat, "The surface extraction from TIN based search-space minimization (SETSM) algorithm," *ISPRS J. Photogramm. Remote Sens.*, vol. 129, pp. 55-76, Jul. 2017.

[139] C. Porter, P. Morin, I. Howat, M.-J. Noh, B. Bates, K. Peterman, S. Keesey, M. Schlenk, J. Gardiner, K. Tomko, M. Willis, C. Kelleher, M. Cloutier, E. Husby, S. Foga, H. Nakamura, M. Platson, M. Wethington, Jr, C. Williamson, G. Bauer, J. Enos, G. Arnold, W. Kramer, P. Becker, A. Doshi, C. D'Souza, P. Cummens, F. Laurier, and M. Bojesen, "ArcticDEM," Nov. 2018.

[140] RGI Consortium, "Randolph glacier inventory - a dataset of global glacier outlines: Version 6.0: Technical report, global land ice measurements from space, colorado, USA," 2017.

[141] ESA, "Land cover CCI product user guide version 2. tech. rep," 2017.

[142] J. Höhle and M. Höhle, "Accuracy assessment of digital elevation models by means of robust statistical methods," *ISPRS J. Photogramm. Remote Sens.*, vol. 64, no. 4, pp. 398-406, Jul. 2009.

[143] B. K. P. Horn, "Hill shading and the reflectance map," *Proc. IEEE*, vol. 69, no. 1, pp. 14-47, Jan. 1981.

[144] L. W. Zevenbergen and C. R. Thorne, "Quantitative analysis of land surface topography," *Earth Surf. Processes Landforms*, vol. 12, no. 1, pp. 47-56, Jan. 1987.

[145] M. B. Wilk and R. Gnanadesikan, "Probability plotting methods for the analysis of data," *Biometrika*, vol. 55, no. 1, pp. 1-17, Mar. 1968.

[146] P. A. Dowd, "The variogram and kriging: Robust and resistant estimators," in *Geostatistics for Natural Resources Characterization*, G. Verly, M. David, A. G. Journel, and A. Marechal, Eds. *Dordrecht: Springer Netherlands*, 1984, pp. 91-106.

[147] M. Mälicke and H. D. Schneider, "Scikit-GStat 0.2.6: A scipy flavored geostatistical analysis toolbox written in python," Nov. 2019.

[148] G. Matheron, *Les variables régionalisées et leur estimation: une application de la théorie de fonctions aléatoires aux sciences de la nature*. Masson et CIE, 1965, vol. 4597.

[149] N. Cressie and D. M. Hawkins, "Robust estimation of the variogram: I," *Journal of the International Association for Mathematical Geology*, vol. 12, no. 2, pp. 115-125, Apr. 1980.

[150] M. G. Genton, "Highly robust variogram estimation," *Math. Geol.*, vol. 30, no. 2, pp. 213-221, Feb. 1998.

[151] T. M. Burgess and R. Webster, "Optimal interpolation and isarithmic mapping of soil properties," *J. Soil Sci.*, vol. 31, no. 2, pp. 315-331, Jun. 1980.

[152] S. Müller and L. Schüler, "GeoStat-Framework/GSTools: v1.3.3 'pure pink'," 2021.

[153] E. J. Pebesma and C. G. Wesseling, "gstat: A program for geostatistical modelling, prediction and simulation," *Comput. Geosci.*, vol. 24, no. 1, pp. 17-31, 1998.

[154] R. Webster and M. A. Oliver, *Geostatistics for Environmental Scientists*, 2nd ed. Wiley, Nov. 2007.

[155] H. Liu and K. C. Jezek, "Investigating DEM error patterns by directional variograms and fourier analysis," *Geogr. Anal.*, vol. 31, no. 3, pp. 249-266, Sep. 2010.

[156] E. Rodríguez, C. S. Morris, and J. E. Belz, "A global assessment of the SRTM performance," *Photogrammetric Engineering & Remote Sensing*, vol. 72, no. 3, pp. 249-260, 2006.

[157] xdem contributors, "xdem," May 2021.

<figure><img src="image_10.png" /></figure>

**Romain Hugonnet** received his MSc. degree in 2017 from Ecole Centrale de Lille, France. Since 2018, he has been working as a Ph.D. student for the Laboratory for Space Studies in Geophysics and Oceanography, Toulouse, France and the Laboratory of Hydraulics, Hydrology and Glaciology, ETH Zürich, Switzerland. His research aims at monitoring the Earth system by remote sensing, with a focus on glaciology and digital elevation models. His fields of expertise are satellite imagery and spatiotemporal analyses, including uncertainty quantification.

<figure><img src="image_11.png" /></figure>

**Fanny Brun** received her Ph.D. degree in 2018 from the University Grenoble Alpes, France. Since 2020, she has been working for IRD in the Institute of Environmental Geosciences, Grenoble, France. Her research aims at observing glacier changes from remote sensing and field techniques, with a focus in Asia. She currently serves as Scientific Editor for the Journal of Glaciology.

<figure><img src="image_12.png" /></figure>

**Etienne Berthier** received his Ph.D. degree in 2005 from the University of Toulouse, France. Since 2007, he has been working for CNRS in the Laboratory for Space Studies in Geophysics and Oceanography, Toulouse, France. His research aims at observing and understanding the effect of climate change on mountain glaciers and ice caps and quantifying their contribution to sea level rise. His field of expertise is satellite remote sensing, including the retrieval of ice topography from optical images to measure ice elevation changes. He currently serves as the Scientific Editor for *The Cryosphere and La Météorologie* (in French). He is member of the European Geophysical Union and the International Glaciological Society.

<figure><img src="image_13.png" /></figure>

**Amaury Dehecq** received his Ph.D. degree in remote sensing in 2015 from the University Grenoble Alpes, France. From 2019 to 2021, he worked as a postdoctoral researcher at ETH Zurich's Laboratory of Hydraulics, Hydrology and Glaciology (VAW). Since 2021, he has been working for the French National Research Institute for Sustainable Development (IRD) and based at the Institute for Geosciences and Environmental research (IGE), Grenoble, France. His research focuses on estimating 20th century mountain glacier changes from remote sensing observations, using a combination of modern satellite, declassified spy satellite, and historical airborne/terrestrial images.

<figure><img src="image_14.png" /></figure>

**Erik Schytt Mannerfelt** received his MSc degree in 2020 fromUiT — The Arctic University of Norway. Between 2020 and 2022, he has worked at the Laboratory of Hydraulics, Hydrology and Glaciology, ETH Zurich, Switzerland, as a research assistant. His research focuses on fieldwork-, modelling-, and photogrammetry-based methods to constrain late Little Ice Age to contemporary glacier change and dynamics.
<figure><img src="image_15.png" /></figure>

**Nicolas Eckert** received his Ph.D. degree in 2007 from AgroparisTech. Since 2008 he is researcher at INRAE Grenoble, co-head of the mountain risk team since 2018. He is also task officer in charges of Environmental risks for the AQUA department at INRAE, and for the ALLENCVI research federation. He serves as Associate Chief Editor for Journal of Glaciology and as Scientific Editor for Cold Regions Science and Technology. His research is at the crossroads between geosciences and statistical modelling, with applications to mountain risks,

mountain climatology and glaciers. Additional interests in the socio-historical component of risk makes him involved in interdisciplinary research addressing all dimensions of mountain risks.

<figure><img src="image_16.png" /></figure>

**Daniel Farinotti** received his doctoral degree in 2010 from the Swiss Federal Institute of Technology in Zurich (ETH Zurich). Since 2016 he leads the Professorship of Glaciology at ETH Zurich's Laboratory of Hydraulics, Hydrology and Glaciology (VAW), a position jointly affiliated to the Swiss Federal Institute for Forest, Snow and Landscape Research WSL. His research focuses on the evolution of glaciers and the implications for water resources, notably including the estimation of glacier ice thickness from surface characteristics, the long-

term modelling of glacier mass budgets, the estimation of the runoff contributions from glacierized catchments, or the implications for water resource management in high-mountain environments.