# TICOI: an operational Python package to generate regular glacier velocity time series

Laurane Charrier, Amaury Dehecq, Lei Guo, Fanny Brun, Romain Millan, Nathan Lioret, Luke Copland, Nathan Maier, Christine Dow, Paul Halas

► To cite this version:

Laurane Charrier, Amaury Dehecq, Lei Guo, Fanny Brun, Romain Millan, et al. TICOI: an opera-tional Python package to generate regular glacier velocity time series. The Cryosphere, 2025, 19 (10), pp.4555-4583. 10.5194/tc-19-4555-2025. hal-05376950

HAL Id: hal-05376950

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

Submitted on 21 Nov 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# TICOI: an operational Python package to generate regular glacier velocity time series

Laurane Charrier<fnref n="1" />, Amaury Dehecq<fnref n="1" />, Lei Guo<fnref n="1" />,<fnref n="2" />, Fanny Brun<fnref n="1" />, Romain Millan<fnref n="1" />, Nathan Lioret<fnref n="1" />, Luke Copland<fnref n="3" />, Nathan Maier<fnref n="4" />, Christine Dow<fnref n="5" />, and Paul Halas<fnref n="6" />

**Correspondence:** Laurane Charrier (laurane.charrier@univ-grenoble-alpes.fr)

Received: 1 November 2024 – Discussion started: 20 January 2025

Revised: 18 June 2025 – Accepted: 22 July 2025 – Published: 15 October 2025

**Abstract.** Glacier velocity is a crucial observation as it controls the mass redistribution and future evolution of the geometry of a glacier. While glacier annual velocities are now available in open data worldwide, sub-annual velocity time series are still highly uncertain and available at heterogeneous temporal resolutions. This hinders our ability to understand flow processes such as basal sliding and surges, as well as the integration of these observations into numerical models. The latest could help to better constrain future projections of sea level rise. We introduce an open source and operational Python package called TICOI (Temporal Inversion using linear Combinations of Observations, and Interpolation). TICOI fuses multi-temporal and multi-sensor image-pair velocities produced by different processing chains, using the temporal closure principle. In this article, we provide extensive examples of TICOI applications on the ITS_LIVE dataset and in-house velocity products, to generate monthly velocity time series. The results are evaluated against GNSS data collected on three glaciers with different dynamics in Yukon and western Greenland, including a surging glacier. Comparison with GNSS observations demonstrates a reduction in error by up to 50% in comparison with the raw image-pair velocities and other post-processing methods. This increase in performance comes from the development of methodological strategies to enhance TICOI's robustness to temporal decorrelation and abrupt non-linear changes. In addition, TICOI can retrieve monthly velocity using annual

image-pair velocities only, when there is sufficient temporal redundancy. This package opens the door to the harmonization of various datasets, enabling the creation of standardized sub-annual velocity products.

## 1 Introduction

Glacier velocity monitoring is key to understand how glaciers are changing in a warming climate. Global glacier thinning (The GlaMBIE Team, 2025; Hugonnet et al., 2021) is expected to lead to large changes in ice flow, which, in turn, affect glacier mass redistribution and future geometry (Dehecq et al., 2019; Morlighem et al., 2011). Documenting ice velocity and changes in ice velocity is thus important to understand the future evolution of glaciers. Regional to global mapping of glacier surface velocity from remote sensing helped to constrain estimates of glacier thickness (Millan et al., 2022) and ice fluxes into the ocean (Kochtitzky et al., 2022; Mouginot et al., 2019; Gardner et al., 2018), which led to reductions on present and future glacier contribution to sea level rise.

While these applications rely on multi-annual velocity averages, many glacier processes should be studied at sub-annual scales, in particular to study calving events (Provost et al., 2024; Riel et al., 2021), glacier response to lake drainage (Maier et al., 2023; Main et al., 2023; Wendleder

<footnote n="1">¹Institut des Géosciences de l'Environnement, Université Grenoble-Alpes, Grenoble, France</footnote>

<footnote n="2">²School of Geo-science and Info-physics, Central South University, Changsha, China</footnote>

<footnote n="3">³Department of Geography, Environment and Geomatics, University of Ottawa, Ottawa, Canada</footnote>

<footnote n="4">⁴Los Alamos National Laboratory, Los Alamos, New Mexico</footnote>

<footnote n="5">⁵Department of Geography and Environmental Management, University of Waterloo, Waterloo, Canada</footnote>

<footnote n="6">⁶independent researcher, formerly at: Bjerknes Centre for Climate Research, University of Bergen, Bergen, Norway</footnote>
et al., 2024), surface runoff (Wendleder et al., 2024), changes in the efficiency of subglacial hydrological networks (Nanni et al., 2023; Maier et al., 2023) or glacier surges (Beaud et al., 2022; Quincey et al., 2015; Copland et al., 2011). Sub-annual ice velocity is crucial for the monitoring of natural hazard such as glacier detachment (Kääb et al., 2021; Gilbert et al., 2018), glacier lake outburst flood (Rashid et al., 2020; Round et al., 2017) and potentially volcanic eruption, when glaciers are located nearby (Martin et al., 2025). Additionally, several recent methodological and modeling developments stress the need for precise and temporally resolved velocity products to infer basal conditions (Goldberg et al., 2015; Jay-Allemand et al., 2011) or near-future projection (Choi et al., 2023), for example using transient inverse methods (Choi et al., 2023; Goldberg et al., 2015). This could help to better constrain future sea-level rise.

So called image-pair velocities are derived by calculating the displacement of features between two images. Recent improvements in satellite image quality and resolution, with the launch of Landsat-8 in 2013, Sentinel-1 in 2014 and Sentinel-2 in 2015, among others, have made it possible to derive image-pair velocities, with an enhanced signal-to-noise ratio at a relative high frequency (5 to 16 d at the best). Despite the recent increase in available velocity datasets, post-processed time series of sub-annual glacier surface velocities, sampled on fixed time intervals, are not yet available at a global scale. Only individual image-pair velocities have been released, such as those from the ITS_LIVE project, which provides velocities derived from correlating Landsat-4,5,6,7,8, and Sentinel-1,2 image pairs, separated by temporal baselines of 5 to 546 d. These products remain noisy and sparse, especially over mountain glaciers. Moreover, image-pair velocities are difficult to interpret because they contain velocities measured using image pairs from different satellites, with different temporal baselines. For these reasons, many research teams processed their own image-pair velocities instead of using existing datasets (Provost et al., 2024; Nanni et al., 2023; Halas et al., 2023; Wallis et al., 2023; Yang et al., 2022; Beaud et al., 2022; Derkacheva et al., 2020), to reach a better signal-to-noise ratio and interpretability. However, this requires high computational and storage resources, and often suffers from a lack of reproducibility. For the same reason, image-pair velocities with variable dates are difficult to include into models, whereas seasonal velocity could aid in retrieving basal conditions (Derkacheva et al., 2021) or ice rheology properties (Bolibar et al., 2023). This calls for a standardized framework dedicated to the processing of available image-pair velocity products to produce consistent sub-annual time series of glacier surface velocity.

Several methods have been proposed to produce velocity time series: cubic spline or LOWESS regression (Derkacheva et al., 2020), regression using a dictionary of B-splines (Riel et al., 2021), sinusoidal regression (Greene et al., 2020), Bayesian recursive smoother (Wallis et al., 2023), iterative weighted monthly averaging (Van Wyk de Vries and Wick-

ert, 2021), or temporal closure (Altena et al., 2019; Charrier et al., 2022b). However, most of these methods use only a subset of the available data, such as velocities quantified using short temporal baselines (< 100 d) (Nanni et al., 2023; Wallis et al., 2023; Riel et al., 2021; Derkacheva et al., 2020), or use strong assumptions about glacier behavior, such as sinusoidal variations in seasonal motion (Greene et al., 2020). Moreover, even though efforts have been undertaken to compare some of the image correlation algorithms (Zheng et al., 2023; Jawak et al., 2018; Heid and Kääb, 2012), no study has tried to fuse datasets produced by different research teams to obtain sub-annual time series. Finally, few of these methods are open source.

The objective of this study is to propose an open source and operational package able to fuse image-pair velocities computed using different satellite images, with different temporal baselines, and possibly using different processing chains, in order to obtain regular (i.e., sampled at regular time steps) velocity time series with an associated uncertainty and relative quality indicators. To do this, we rely on a method based on the temporal closure of the displacement network, also known as a Small BAseline Subset (SBAS)-like approach, which originated from the Interferometric Synthetic Aperture Radar (InSAR) community (Doin et al., 2011; Berardino et al., 2002). These approaches have previously been adapted and applied to glaciers to retrieve 2D or 3D velocity time series (Provost et al., 2024; Charrier et al., 2022b, a; Samsonov et al., 2021; Guo et al., 2020), but have not been applied in an operational framework because: (1) they remain sensitive to temporal decorrelation (i.e., very low velocity values which can be measured when strong surface changes occur, or near the margins of the moving object. This phenomenon is especially visible when the temporal baseline is large and/or when the velocity is low.); (2) they often include a regularization term which assumes acceleration to be zero in time (which is not true for surging glaciers or glaciers with a strong seasonality pattern), or a pre-defined mathematical model (i.e. a periodical model) which requires a priori knowledge on the glacier dynamic; (3) they have never been evaluated against Global Navigation Satellite System (GNSS) data; and (4) computational costs are high.

To overcome these issues, we present a new method called Temporal Inversion using linear Combinations of Observations, and Interpolation (TICOI) to derive regular glacier velocity time series. First, we detail our methodological strategy to increase the robustness of previous developments against temporal decorrelation and abrupt non-linear changes. We also propose three criteria to evaluate the quality of our results, including an error propagation, that is based on a robust theoretical framework. Then, we validate this new approach against seven GNSS (Global Navigation Satellite System) stations located on three different glaciers, including one with an active surge phase, and two others showing seasonal variations. We carry out a sensitivity analysis, and illustrate our results along glacier center flowlines. Then, we
show that TICOI is able to retrieve sub-annual velocity time series from annual image-pair velocities only. To finish, we discuss the interest of fusing datasets originating from different processing chains, and the potential application of TICOI at a regional and global scale.

## 2 Method

In this section, we describe the TICOI workflow (Fig. 1). The first part builds on previous developments (Charrier et al., 2022b, a), the second part presents methodological strategies to improve the computational performance and robustness of the method to temporal decorrelation and abrupt non-linear changes in surface velocities. The aim is to provide a robust post-processing package that does not require a priori knowledge of the ice flow behavior.

### 2.1 Previous developments

#### 2.1.1 Temporal closure's principle

The temporal inversion is based on the temporal closure of the displacement measurement network (Berardino et al., 2002). Temporal closure links n measured displacements in Y to p estimated displacements in X (Fig. 1), by a system of linear equations. This system of linear equations can also be written as AX = Y, with A the design matrix linking X to Y. To understand the structure of A, let's take an example with three displacements represented in Fig. 1. Assuming that the displacement is cumulative in time, it can be written that:

$$ 
\left\{ 
\begin{aligned}
d_{t_0,t_6} &= \hat{d}_{t_0,t_4} + \hat{d}_{t_4,t_5} + \hat{d}_{t_5,t_6} \\
d_{t_0,t_4} &= \hat{d}_{t_0,t_4} \\
d_{t_4,t_6} &= \hat{d}_{t_4,t_5} + \hat{d}_{t_5,t_6}
\end{aligned}
\right. 
\qquad (1) $$

with $d_{t_i,t_j}$ a measured displacement between dates $t_i$ and $t_j$, and $\hat{d}_{t_k,t_l}$ a displacement estimated between dates $t_k$ and $t_l$. This is equivalent to $AX = Y$, which is in this example:

$$ \begin{bmatrix} 1 & 1 & 1 \\ 1 & 0 & 0 \\ 0 & 1 & 1 \end{bmatrix} \begin{bmatrix} \hat{d}_{t_0,t_4} \\ \hat{d}_{t_4,t_5} \\ \hat{d}_{t_5,t_6} \end{bmatrix} = \begin{bmatrix} d_{t_0,t_6} \\ d_{t_0,t_4} \\ d_{t_4,t_6} \end{bmatrix} \qquad (2) $$

It is important here to highlight the temporal redundancy between the measured displacements $d_{t_0,t_6}$, $d_{t_0,t_4}$ and $d_{t_4,t_6}$. The estimated displacement $\hat{d}_{t_0,t_4}$ can be obtained from these three displacements. This is the main interest of using temporal inversion: reducing the noise by using the redundancy between displacement measurements computed with different temporal baselines. Note that X contains either the East/West displacements or the North/South displacements.

#### 2.1.2 Inversion of the system

Most of the time, the system $AX = Y$ is ill-posed, i.e. $rank(A) < p$ (the number of linearly independent rows of

the matrix A is lower than the number of estimations p), and the system has an infinite number of solutions. To overcome this problem, the system can be solved either with a Singular Value Decomposition (SVD) (Berardino et al., 2002), a Least Square (LS) approach (Bontemps et al., 2018; Samsonov and d'Oreye, 2017; Doin et al., 2011) or a L1-norm solution (Lauknes et al., 2010). The L1-norm is more robust to outliers, but computationally expensive, as it requires computing the absolute of the residuals, which is not a differentiable piece-wise function. The SVD solution is equivalent to the minimum-norm LS solution (i.e. it tends to minimize the norm of X) (Berardino et al., 2002). In order to have a more flexible regularization strategy, we use a Weighted Least Square (WLS) approach. The cost function is:

$$ \operatorname{argmin} \left( \|W(AX - Y)\|^2 + \lambda \|\Gamma(X - X_0)\|^2 \right) \qquad (3) $$

where W is a $n \times p$ matrix standing for the weight given to each value in Y, $\lambda$ is a scaling constant and $\Gamma$ is a $p \times p$ matrix representing the regularization matrix, and $X_0$ is first guess solution detailed in Sect. 2.2.2.

Different regularization matrix $\Gamma$ and weights W can be used. The choice of $\Gamma$ will be discussed in Sect. 2.2.2. As for the weights, W could be equal to the identity (Berardino et al., 2002), in which case the WLS solution is equivalent to an Ordinary Least Square. But this ignores the heteroscedasticity of the measurements (i.e. the fact that they have unequal variances). Therefore, it is common to use a priori knowledge of the data quality, for example the InSAR coherence (Yunjun et al., 2019), the shape of the similarity map used in image correlation (Bontemps et al., 2018), errors computed over stable areas, the cosine of the angle between each displacement vector and the spatio-temporal median, or the modified zscore (Charrier et al., 2022a). However, these metrics do not always accurately represent the errors of the measurements. Another complementary strategy is to use the inverse of the residuals (Bontemps et al., 2018; Liang et al., 2020), which can be expressed as the vector of dimension n:

$$ R = A\hat{X} - Y \qquad (4) $$

In the TICOI workflow, the weights can be a combination of a priori knowledge of the data quality (e.g., image correlation score or velocity in stable areas Gardner et al., 2018) and residuals. To be more robust to outliers, we use an Iterative Re-weighted Least Square (IRLS) approach, i.e. we update the weights iteratively using the residuals from the previous inversion.

In the first iteration of the inversion, if the data quality is known, the diagonal elements of the weights, $W_0$, corresponds to the errors scaled between 0 and 1.

In the second iteration, and all following iterations, each diagonal element of the weighted matrix, at the position $(m, m)$ and iteration $u$ is updated as:

$$ W_{m,m}^{u} = \psi(Z_{m}^{u}, c) \qquad (5) $$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** TICOI workflow. In this example, displacements have been measured using images from satellites 1 and 2, which have a repeat cycle of 2 and 5 d, respectively. Outliers have been removed in the displacement dataset, for example from areas impacted by clouds (shown in yellow on the image subsets). The TICOI workflow is applied pixelwise (in this example on the pixel P). First, the system AX = Y is inverted to obtain an irregular time series (Sect. 2.1.2). Second, the time series is interpolated to obtain a regular time series with a homogeneous and optimal temporal sampling (Sect. 2.1.3). Note that the temporal sampling of the final time series τ can be chosen by the end-user.</figcaption></figure>

In this equation, Z is a standardized residual vector of dimension n computed as:

$$Z = \frac{R}{\text{NMAD}(R)} \qquad (6)$$

with NMAD being the Normalized Median Absolute Deviation of the residuals, equal to 1.483 MAD.

$\psi$ is the Tukey's biweight function, which is a common down-weight function (Liang et al., 2020) robust to large outliers, defined as:

$$\psi(Z_m^u, c) = \begin{cases} [1 - (Z_m^u/c)^2]^2, & |Z_m^u| < c \\ 0, & |Z_m^u| > c \end{cases} \qquad (7)$$

where c is a tuning constant which is usually set to 4.685, producing 95% efficiency for a normal distribution (Huber, 1992).

The iterations stop when (mean (|X̂u - X̂u-1|) < δ or (u > 10)) where X̂u corresponds to the results of a given iteration u and X̂u-1 the results of the previous one. δ is set to 0.1 m.

### 2.1.3 Interpolation of the irregular time series

The inversion results in an irregular displacement time series, i.e. the vector X contains displacement between each measured dates (all dates with an image acquisition minus dates rejected by outlier removal) (Fig. 1). However, we need

regular time series to study glacier dynamics. To study variations of fast moving objects such as glaciers we are interested in velocity time series (Derkacheva et al., 2020; Greene et al., 2020), more than cumulative displacement time series (Lacroix et al., 2019; Doin et al., 2011). When looking at velocities, it is necessary to compare velocities with the same temporal sampling for three main reasons: (1) velocities with different temporal sampling are not comparable because they correspond to the average of the instantaneous velocity over different time intervals (e.g., annual image pair-velocities are close to the annual average of the instantaneous velocity whereas short temporal baseline velocity are close to the instantaneous velocity) (Charrier et al., 2022a); (2) velocities with very short temporal sampling can be very noisy because noise decreases with increasing temporal sampling (Millan et al., 2019); and (3) the dates between which the displacements are inverted can be different from one pixel to another because of outliers removal.

To overcome this problem, some studies have proposed to use fractions of displacements inside the network (Samsonov et al., 2021). However, this assumes the velocity to be constant over the corresponding temporal sampling. This assumption is most of the time inaccurate, particularly over long time periods or in cases of surge behavior or seasonal variations, as demonstrated in Charrier et al. (2022c). Therefore, we propose to: (1) compute the cumulative displacement time series by summation of the irregular time series;
(2) interpolate this cumulative displacement time series (here using a cubic spline); and (3) obtain a regular time series with a given temporal sampling, using a discrete derivative. See Charrier et al. (2022a) for more details.

The resulting time series has a constant temporal sam-
pling, for example 4 d in Fig. 1. The larger the temporal sam-
pling, the smoother the time series, i.e. the signal-to-noise
ratio increases, but the temporal resolution decreases. By an-
alyzing the Root Mean Square Error (RMSE) over stable ar-
eas, we have shown that the RMSE according to the tem-
poral sampling has an asymptotic behavior which converges
after around 30 d for glaciers with medium average veloc-
ity (~100 to 200 m yr⁻¹) (Charrier et al., 2022a, b). Note
that this asymptote could be reached with a smaller temporal
sampling for faster glaciers.

## 2.2 Improved robustness and computational performance

An IRLS with a first order Tikhonov regularization term per-
forms poorly in some extreme cases, such as temporal decor-
relation or abrupt non-linear changes, especially when there
are few image-pair velocities. Below, we show improvements
to the method that overcome these challenges.

2.2.1 Robustness to temporal decorrelation

Using a LS approach to solve $AX = Y$ assumes the errors of $Y$ to be normally distributed. However, this assumption is not always true in a real case scenario (Fig. A1). Robust LS regression, like IRLS using Tukey's bi-weight function, helps to reduce the effect of outliers in case of random errors (Charrier et al., 2022b; Liang et al., 2020) but may be inefficient for systematic errors. For example, when temporal decorrelation occurs, the measured displacement is systematically close to 0, instead of the true glacier velocity, which results in a heavy-tailed distribution of errors with a strong kurtosis (Fig. A1). To overcome this problem, we propose to carry on a first LS with short temporal baselines only (lower than 180 d), to automatically detect temporal decorrelation. We compute the residual between each observation (with short and long temporal baselines), to this first small baseline LS solution. Displacements impacted by temporal decorrelation have large residuals, because they do not satisfy temporal closure.

2.2.2 New regularization term robust to abrupt non-linear changes

A crucial choice is the regularization matrix ($\Gamma$ in Eq. 3). The most common regularization matrix is the first order Tikhonov regularization, which assumes the acceleration to be null in time (Charrier et al., 2022b; Samsonov et al., 2021; Lacroix et al., 2019; Bontemps et al., 2018). For this regularization, $\Gamma$ corresponds to a first order derivative operator. The diagonal elements of $\Gamma$ in Eq. (1) are: $\Gamma_{k,k} = 1/\Delta\tau$

and the element above the diagonal will be $\Gamma_{k,k+1} = -1/\Delta\tau$
with $\Delta\tau$ the temporal sampling of time series. When us-
ing the Tikhonov regularization, $X_0$ is a vector filled with
0 in the cost function shown in Eq. (3). In other words,
this strategy minimizes the difference between adjacent esti-
mated velocities, i.e. the acceleration. However, the assump-
tion of an acceleration close to 0 is not always valid, espe-
cially when abrupt non-linear changes occur (e.g., for surge-
type glaciers, or those with strong spring speedups). Pepe
et al. (2016) and López-Quiroz et al. (2009) proposed to use
a model-based regularization, but it also requires a priori
knowledge of the displacement behavior.

Therefore, we propose a new regularization strategy. We constrain the acceleration of the estimated time series, X, with the acceleration of an initial guess, X₀. In the regularization term ||Γ(X − X₀)||²) from Eq. (3), Γ is identical to the Tikhonov regularization matrix but X₀ is not null anymore. The latest is estimated from the spatio-temporal smoothing of linearly interpolated measurements, to ensure that X₀ and X correspond to the same time intervals. We have tested different filters and have obtained the best results with a spatial window of 3 × 3 pixels and a third-order Savitzky-Golay filter with a temporal window of 90 d (see Appendix B2). We use only velocity measurements with temporal baselines <180 d to avoid long temporal baselines, which will be close to the long-term average of velocity, resulting in an oversmoothed solution (see Fig. A2).

2.2.3 Improved computation time

When velocities are derived from images acquired by dif-
ferent satellites, spanning different temporal baselines (e.g.,
from 5 to 546 d in the ITS_LIVE dataset), the length of *Y* can
be very large (e.g., on the order of 10 000 over Kaskawulsh
Glacier, Yukon, between 2013 and 2022). This will result in
very long computation time (several seconds per pixel). To
mitigate this challenge, we solve the Least Square problem
using LSMR, a conjugate-gradient method for sparse least-
squares problems, that leverages the fact that the matrix *A* is
generally sparse, i.e. contains mainly 0. Additionally, we im-
plement embarrassingly parallel processing at the pixel level.
Besides, the lazy mode from Dask, an open source Python li-
brary for parallel computing, allows the user to directly apply
TICOI on ITS_LIVE data sets stored in the Amazon cloud
(Gardner et al., 2025), without the need to download the data
locally. Dask also allows out-of-memory processing by split-
ting the data in chunks.

2.3 Automatic selection of the regularization coefficient

The choice of the regularization coefficient $\lambda$ in Eq. (3) can be empirical (Lacroix et al., 2019; Bontemps et al., 2018) or based on an L-curve (Samsonov et al., 2021), which aims to compare $||W(AX - Y)||^2$ and $||\Gamma(X - X_0)||^2$ in Eq. (3). However, L-curves do not always converge (Vogel, 1996). As
an alternative solution, we propose to use the Velocity Vector Coherence (Dehecq et al., 2015) defined as:

$$VVC = \text{mean}_{(i,j) \in \omega} \left\| \sum_{t=0}^{N} \frac{V(i, j, t)}{\|V(i, j, t)\|} \right\| \quad (8)$$

with $\omega$ the area over which the VVC is computed. The VVC varies from 0 to 1, with 1 corresponding to a perfect coherence of the direction in time.

When the regularization coefficient increases, the solution tends to be smoother. Therefore, the direction of the velocity vectors tends to be constant in time, and the VVC tends to be close to 1. To find a compromise between smoothing and improved signal-to-noise ratio, we need to find the inflection point of this function, which is approximated by:

$$VVC_{\text{optimal\_coef}} = \max(\text{VVC}) - 0.05 \cdot \left( \max(\text{VVC}) - \min(\text{VVC}) \right) \quad (9)$$

If there is a real change in velocity direction over time, the VVC curve will converge before 1, but the curve will still have an inflection point. For better robustness, it is advisable to calculate the VVC over a relatively large area.

## 2.4 Uncertainty

Three metrics are proposed to evaluate the uncertainties of the estimated velocity time series: (1) the VVC of the time series; (2) the number of image-pair velocities that have contributed to each estimation; and (3) the a posteriori covariance matrix.

The VVC has been defined in Eq. (8). The number of image-pair velocities that have contributed to each estimation is defined as:

$$X_{\text{count}} = A^T W \quad (10)$$

This is computed for North/South and East/West components separately. Then, the number of image-pair velocities of the velocity magnitude is taken as the average number of image-pair velocities of two velocity components.

The a posteriori covariance matrix is defined by assuming the errors to be independent (Gavin, 2023; Liang et al., 2020):

$$\Sigma_{\hat{X}} = N^{-1} A^T W \Sigma_{\hat{Y}} W A N^{-1} \quad (11)$$

with $N = A^T W A + \lambda \Gamma^T \Gamma$, which contains a data fidelity term $A^T W A$, and a regularization term $\lambda \Gamma^T \Gamma$. $\Sigma_{\hat{Y}}$ is the covariance matrix of the image-pair velocities, which contains the square of the errors provided with the raw image-pair velocities, converted in meters. If the errors are independent, this matrix is diagonal. The demonstration of this formula is provided in Appendix A.

The a posteriori covariance matrix is computed for each of the components separately and interpolated using the same strategy described in Sect. 2.1.3. Then, the a posteriori covariance matrix of the velocity magnitude is computed following the propagation of uncertainty:

$$\Sigma_{\hat{v}} = \sqrt{\left(\frac{v_x}{v} \cdot \Sigma_{\hat{v}_x}\right)^2 + \left(\frac{v_y}{v} \cdot \Sigma_{\hat{v}_y}\right)^2} \quad (12)$$

with $v_x$ and $v_y$ the x- and y-velocity component, $v$ the velocity magnitude and $\Sigma_k$ the a posteriori covariance matrix of the variable $k$.

Finally, the confidence intervals are defined for each estimated velocity as: $\pm t_{(1-\alpha/2),n-p} \sqrt{\Sigma_{\hat{v}}}$, with $t_{(1-\alpha/2),n-p}$ the value of the student's t-distribution for a degree of freedom of $n-p$ at a confidence level of 100(1-α) %. Here, we chose α to be equal to 0.05, i.e., we provide a 95 % confidence interval.

# 3 Data and study sites

The package is validated on glaciers with low to medium velocity magnitude (~100 to 200 yr⁻¹) (Fig. 2). We chose these glaciers because they have been monitored with continuous GNSS data over multi-year periods and represent both seasonal and surge type dynamics. The possible sites available for this comparison are rare since few glaciers are monitored this way. In the following section we describe each of the glaciers, their associated GNSS data, and surface velocity measurements derived from remote sensing.

## 3.1 Lowell and Kaskawulsh glaciers in Yukon, Canada

Lowell and Kaskawulsh glaciers are two large valley glaciers located in Kluane National Park, at the eastern edge of the St Elias Mountains, in Yukon, Canada. Lowell glacier, also known as Nàlùdäy in Southern Tutchone, is a ~65 km long surge-type glacier. It is composed of a southern and a northern arm divided by a medial moraine. The northern arm joins the main trunk by a ~200 m high icefall, whereas the southern arm originates from a large high-accumulation basin. The terminus of Lowell is divided in two by a large nunatak (Bevington and Copland, 2014). The last 6 surges have been observed in 1948–1950, 1968–1970, 1983–1984, 1997–1998, 2009–2010 (Bevington and Copland, 2014), and 2021–2022 (Van Wychen et al., 2023). During the 5 surges that occurred between 1948 and 2009, the length of the surge active phase ranged from 0.6 to 2 years and the quiescent phase from 11 to 18 years. Bevington and Copland (2014) note that the surge cycle (quiescent + active phase) seems to have decreased in time, which was supported by the start of the most recent surge in 2021. Surges show a rapid terminus advance starting in summer to early fall (late June to early October), continuing through the winter, and ending in June or July of the following year. Velocity peaks during the surge phase are typi-
cally > 3500 m yr⁻¹ in the terminus region, and the fastest-recorded motion is of 11 000 m yr⁻¹ in the lowest part of the south arm during the 1983–1984 surge (Bevington and Cop-land, 2014).

Kaskawulsh, called Tänshj in Southern Tutchone lan-
guage, is approximately 70 km long, with altitudes rang-
ing from 800 to 2500 m above sea level (a.s.l) and flowing
eastwards. It is divided into three main tributaries (North,
Central, and South Arms) (Flowers et al., 2014; Foy et al.,
2011). It is believed to be temperate, at least across its abla-
tion area, and has been monitored for many decades (Clarke,
2014; Flowers et al., 2014; Foy et al., 2011; Arendt et al.,
2002), starting with the Icefield Ranges Research Project
in the 1960s and the 1970s (Anderton, 1973; Clarke et al.,
1967). Annual velocities are about 10 to 50 m yr⁻¹ near the
terminus, and range between 100 and 200 m yr⁻¹ over the
rest of the centerline (Main et al., 2023; Millan et al., 2022;
Waechter et al., 2015). Its surface velocity remained sta-
ble between 1960 and 2012, except for the lower part of
the ablation area (up to ~10 km upstream of the terminus;
Waechter et al., 2015). Annual average velocities across the
lower glacier increased at a rate of 5.5 m yr⁻² between 2010
and 2015, followed by a decrease between 2016 and 2018
of about 8 m yr⁻², mainly in the northern lobe (Main et al.,
2023). The main cause of these changes can be explained by
the drainage of proglacial Slims lake, located northwest of
the terminus, according to Main et al. (2023).

3.2 Land-terminating margin in western Greenland

The land-terminating region in western Greenland is char-
acterized by an absence of marine outlet glaciers where an-
nual velocities are on the order of 100–200 m yr⁻¹ (Joughin
et al., 2018). The considered study area is located inland of
Issunguata Sermia, a land-terminating outlet glacier, particu-
larly well studied over the last two decades using both satel-
lite and in situ instrumentation. Data from this region serves
as the foundation for a better understanding of hydrology-
dynamic coupling in Greenland, where melt-forced velocity
changes have been observed from daily to decadal timescales
(Davison et al., 2019). Surface lake drainage and intense
melt events drive flow accelerations as high as 10 times
above background across daily to weekly timescales as the
subglacial drainage system is temporarily overwhelmed by
the rapid influx of meltwater (Doyle et al., 2015; Tedstone
et al., 2013). Seasonal velocity cycles (two-three times win-
ter velocities) driven by summer melt production have been
well-documented and attributed to basal pressure changes
modulated by melt supply variability and seasonally evolv-
ing drainage (Maier et al., 2022; van de Wal et al., 2015;
Sole et al., 2013; Bartholomew et al., 2012, 2010). Across
decadal timescales, gradual velocity decreases (20 %) have
been found in response to periods of elevated melt rates (Ted-
stone et al., 2015; Halas et al., 2023; Williams et al., 2020).
These variations are hypothesized to be caused by increases

in summer drainage efficiency, which gradually depressur-
izes the ice-base (Williams et al., 2021; Tedstone et al.,
2015).

3.3 Surface velocity measurements

3.3.1 Datasets

We demonstrate the TICOI method with two datasets of surface flow velocity: (1) the *Institut des Geosciences de l'Environnement* (IGE) dataset (Halas et al., 2023; Millan et al., 2022; Halas et al., 2022; Derkacheva et al., 2020; Millan et al., 2019), available on-demand and (2) the NASA MEaSUREs project Inter-mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) dataset (Gardner et al., 2018, 2022), available on-line. The first one derived displacement using a modified version of the Normalized Cross-Correlation (NCC) algorithm AMPCOR from ROI_PAC (Mouginot et al., 2019; Millan et al., 2019). In the Yukon, it is based on images from Sentinel-2 and Landsat-8 and the correlation window size is 16x16 pixels (Mouginot et al., 2023) and corresponds to the dataset published by Millan et al. (2022). In Greenland, the image-pair velocities are based on Landsat-7, Sentinel-2, Landsat-8 and Sentinel-1 data and corresponds to datasets published by Halas et al. (2023); Derkacheva et al. (2020). Image-pair velocities span temporal baselines from 5 to 100 d and from 330 to 400 d in Yukon, from 5 to 32 d (Derkacheva et al., 2020) and from 330 to 400 d over western Greenland (Halas et al., 2023). The spatial sampling of the velocity maps is 50 m in Yukon and 150 m in Greenland. In the dataset, a previous outlier filtering has been performed: displacements that deviate more than three pixels from the median velocity computed over a spatial window of 9 × 9 pixels are considered to be outliers (Millan et al., 2019; Mouginot et al., 2012). The uncertainties depend on the spatial resolution of the images and the temporal baselines as explained in Millan et al. (2019), by assuming the uncertainty in displacement to be around 1/10 pixels.

The second dataset, ITS_LIVE v2.0 (Gardner et al., 2022), contains velocities measured using the NCC algorithm “autoRIFT” (autonomous repeat image feature tracking) on images from Sentinel-1, Sentinel-2, Landsat-7 and 8. Results from the sparse search guide a dense search (Gardner et al., 2018). The size of the correlation window is increased iteratively according to a threshold on the Normalized Displacement coherence, an indicator of the quality of the correlation (Gardner et al., 2018) (i.e., there are uneven correlation window sizes which could lead to uneven smoothing in space). The temporal baselines range between 5 and 546 d. The spatial sampling of the velocity maps is 240 m resampled at a resolution of 120 m using a cubic spline approach (Lei et al., 2022). Only velocities that agree within 4 times the centered 5 × 5 mean absolute deviation are retained. The uncertainty of each image-pair velocity corresponds to the standard error of velocities relative to the stable or slow moving areas.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Study areas of three selected glaciers: Kaskawulsh and Lowell glaciers in the Yukon; Issunguata Sermia Glacier in western Greenland. Black lines represent the center flow-line of Kaskawulsh and Lowell glaciers according to the Randolph Glacier Inventory version 7 (RGI 7.0 Consortium, 2023). It has been slightly modified for Kaskawulsh to reach Slims lake. Orange lines represent the glacier outlines from RGI v7. Blue dots show time evolving locations of GNSS stations. The three main GNSS locations are named Kask L, M, U for the lower, middle and the upper of Kaskawulsh glacier, and Lowell L, M, U for the lower, middle and upper of Lowell glacier. Blue squares represent areas including all GNSS stations for each site. They are used in Sect. 4. The background of each site is the 2017–2018 average velocity from Millan et al. (2022) overlaid over Sentinel-2 images.</figcaption></figure>

The strength of ITS_LIVE is to be available worldwide in open access from 1980s to 2023, while the IGE dataset, published in Millan et al. (2019), covers only two years. The strength of the IGE dataset is its spatial resolution (50 m against 240 m) which allows velocities of relatively small glaciers to be captured. A more detailed comparison is beyond the scope of this paper.

### 3.3.2 Pre-processing of datasets

In the TICOI method, we optionally apply a filtering to the input data. Two main types of strategies for deleting outliers are implemented: the Modified Zcore (Mzscore), which filters the velocity components 3.5 NMAD away from the median of the entire time period (Maronna et al., 2019), and the Median Angle (MA) that removes image-pair velocities that deviate more than 45° from the median vector (Rabatel et al., 2023; Charrier et al., 2022a). In this study, we use the MA filter.

The IGE and ITS_LIVE datasets are reprojected using a nearest-neighbor interpolation on the same coordinate system:

the polar stereographic North projection (EPSG code 3413) for both Greenland and Yukon. This process requires reprojecting both the grid coordinates and the values of the EW and NS velocity components. These velocity components are defined relative to the orientation of the grid, i.e., the EW (NS) component is the projection along the *x* (respectively *y*) axis. Consequently, it is necessary to reproject the grid and recalculate the velocity vector projections along the new axes. To achieve this, we first compute the coordinates of the start and end points of each velocity vector in the new coordinate system. We then calculate the difference between these coordinates along the axes of the reprojected grid to obtain the new component values.

## 3.4 GNSS data

On Kaskawulsh and Lowell glaciers (Yukon), up to six dual frequency GNSS receivers recorded their position at 15 s intervals for 2 or 3 hours per day in winter, and 24 hours per day in summer (Van Wychen et al., 2023; Waechter et al., 2015), from 2013–2022. The positions were post-processed using
Natural Resources Canada's Precise Point Positioning service (http://webapp.geod.nrcan.gc.ca/geod/tools-outils/ppp.php, last access: 1 June 2023), resulting in an accuracy of ~1–2 cm in horizontal and ~5 cm in vertical. We removed any position derived from a daily record of < 1.5 h, which sometimes occurs in winter due to low battery power. Prior to 2017, the stations were Trimble R7 units mounted on poles drilled into the glacier surface. From 2017 onwards, the stations were replaced with Trimble NetR9 units mounted on tripods “floating” on the glacier surface. They were manually moved up glacier every few years to compensate for the glacier displacement, to ensure that they remained in approximately the same location. This creates artifacts in the position time series, which we remove using the Local Outlier Filter (LOF) (Breunig et al., 2000), which computes the local density deviation of a given data point with respect to its neighbors. We apply it on the gradient of the East/West and North/South position. A daily position is considered to be an artifact if LOF > 5. Then, velocities are derived from the discrete derivative of the position time series, and averaged using a temporal window of 5 d, which corresponds to the minimal repeat cycle of the satellites used. Time spans with less than 80% available daily velocities are removed.

In western Greenland, we derived velocities from 15 s position data collected at a field site located 33 km from the terminus of Issunguata Sermia using five Trimble NetR9 GNSS receivers mounted on a poles frozen into the ice. The position data from each receiver was processed against an off-ice base station using TRACK v1.29 differential kinematic processing software (previously detailed in Maier et al., 2022; Maier et al., 2019). From the GNSS positions, ice velocity is estimated on a daily basis (Halas et al., 2022). While the collection period was from 2014–2017, no data are available from each winter due to power limitations.

Since the Yukon GNSS stations typically move around 100–150 m yr⁻¹, we compare them to remote sensing velocities or TICOI results located at the averaged GNSS location of the corresponding year. The GNSS station in Greenland is supposed to be stable spatially because the period of measurements is shorter and the sampling of the velocity maps coarser (considering a spatial sampling of 150 m, an average velocity of 125 yr⁻¹, the GNSS stations have likely moved by 2.5 pixels between 2014–2017). Then, we average daily GNSS velocities to match the same temporal baselines as remote sensing velocities, or the temporal sampling of TICOI time series, which in this study is 30 d.

# 4 Results

## 4.1 Robustness of TICOI method

### 4.1.1 Robustness to temporal decorrelation

One of the main improvements of TICOI compared to previous methods is its ability to eliminate artifacts caused by temporal decorrelation in the input datasets. The strategy detailed in Sect. 2.2.1 rejects long temporal baseline velocities which are biased. These velocities tend to have large errors when compared to the temporal closure solution derived using only short temporal baseline velocities. This results in large residuals (Eq. 4) relative to the initial solution with short temporal baselines, and as a result, Tukey's biweight function assigns a weight of 0 to these raw image-pair velocities.

To evaluate the performance of this strategy, we applied TICOI to a large area around Kaskawulsh Glacier for the period 2013–2022 considering two implementations: with and without an automatic detection of temporal decorrelation. Then, we computed the averaged velocity magnitude of TICOI results obtained with the two implementations (Fig. 3a). The median difference between TICOI with and without an automatic detection of temporal decorrelation is −0.53 m yr⁻¹. However, the difference can reach up to 100 m yr⁻¹ near glacier margins or in narrow, steep areas, where temporal decorrelation is more likely to occur (Fig. 3a), as in point A (Fig. 3c). To assess the uncertainty of each strategy, we compute the NMAD and the median of velocities over stable areas (i.e. areas outside of glaciers boundaries defined by the RGI v.7) for each date of the time series. These indicators are traditionally used to assess the uncertainty of the observations. The maximum and the median value in time of both the NMAD and the median slightly decrease (Fig. A3) with the use of the automatic detection of temporal decorrelation, highlighting a reduction of uncertainty. Note that the conclusions are similar if we calculate the statistics for the EW and NS components instead of the magnitude.

Notably, TICOI remains robust to temporal decorrelation without this specific strategy, when decorrelated image-pairs are in the minority (Fig. 6). In such cases, the Tukey biweight function is able to effectively filter out the decorrelated image-pairs.

### 4.1.2 Robustness to strong changes in velocity

The traditional regularization penalizes abrupt changes in velocity, and is thus unable to resolve the peak of velocity for a surge event. For instance, for the surge of the lower station of Lowell glacier, the traditional regularization retrieves a continuous increase in velocity from mid-2021 onward, while TICOI captures well a sudden increase in velocity and the rapid glacier slow-down in summer 2022 (Fig. 4). In this
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Illustration of the robustness to temporal decorrelation over Kaskawulsh glacier: (a) represents the difference between the averaged velocity magnitude computed using TICOI with and without an automatic detection of temporal decorrelation. Grey lines correspond to glacier outlines according to RGI v.7 (RGI 7.0 Consortium, 2023). (b) ITS_LIVE image-pair velocities on the point A represented in (a). (c) time series of TICOI without (purple) and with an automatic detection of temporal decorrelation (“TICOI_detect_temp”, red) over the same point A. The coordinates of point A is 139.6458° W, 60.6991° N). The coordinates of the grid correspond to the projection EPSG:4326.</figcaption></figure>

example, few velocity measurements are available in 2022–2023, because the outlier removal step of ITS_LIVE rejected many image-pair velocities (see Sect. 3.3), therefore the solution is strongly impacted by the regularization term, which minimizes the acceleration in the traditional approach. The TICOI regularization (described in Sect. 2.2.2) relaxes this assumption of minimal acceleration by using an initial guess for the acceleration, which makes it possible to capture the peak of the surge even with limited image-pair velocities (Fig. 4).

## 4.2 Validation with GNSS time series for different glacier dynamics

We evaluate TICOI time series with seven GNSS time series on three different glaciers (Table 1). The metrics for the evaluation are the RMSE (in m yr⁻¹) and Kling–Gupta efficiency (KGE) (no unit) (Gupta et al., 2009). The KGE is a goodness-of-fit indicator, widely used to calibrate hydrological models, in order to make sure that they well capture peak

flows, as well as the seasonality of the flow. It is defined as:

$$KGE = 1 - \sqrt{(r - 1)^2 + (\alpha - 1)^2 + (\beta - 1)^2} \quad (13)$$

with r the Pearson's coefficient, α = σe/σo with σe and σo the variance of the estimated and observed time series, respectively, and β = μe/μo with μe and μo the mean of the estimated and observed time series, respectively. α represents the variability of the estimation and β is the bias term. The KGE values range between -∞ and 1. A perfect agreement between two time series would lead to a KGE of 1.

TICOI leads to a reduction in RMSE from 9%, in the least favorable case, to 69%, in the most favorable case, in comparison with the image-pair velocities. The median improvement is of 52% (corresponding to a reduction from 5 to 47 m yr⁻¹, with a median reduction of 34 m yr⁻¹). The KGE increases up to 1.87, with a median improvement of 0.57 (Table 1). This improvement is less important for the Lowell L and M stations, due to a larger density of TICOI estimation during the surge, when the velocity magnitude is higher, in comparison with the image-pair velocities. Besides, TICOI leads to a reduction of RMSE from 11 to 81% in comparison
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) Image-pair velocities from ITS_LIVE at Lowell L. Dots and bars represent the central date and temporal baseline, respectively, of each image-pair velocity. (b) Example of velocity time series retrieved on the lower part of Lowell glacier (Lowell L GNSS station) with TICOI regularization term based on an initial guess about the acceleration (red crosses) and with the traditional regularization term (blue triangles), i.e., the first order Tikhonov regularization. The orange dots show the GNSS data averaged to match TICOI temporal sampling.</figcaption></figure>

with a rolling median, with a median improvement of 40% (3 to 250, 14 m yr⁻¹ in median). The KGE increase up to 0.44, with a median improvement of 0.27. For the point Lowell U, the rolling median provides better results than TICOI. This could be due to the low signal-to-noise ratio of image-pair velocities over this area. A strategy to improve this result is discussed in Sect. 5.2.

We note that a rolling median applied on all temporal baseline velocities gives slightly better RMSE than a rolling median with small baselines only, for some of the considered points, as illustrated in Fig. A2. However, the solution using all baselines tends to underestimate the largest velocities, because they include annual temporal baselines which are close to the annual average of the signal. These results show the strength of TICOI: taking advantage of all temporal baselines, without over-smoothing the solution.

## 4.3 Sensitivity analysis and automated choice of the hyperparameters

TICOI was developed as a flexible method, with processing options that can be changed by the user. Several options can be modified: the coefficient of regularization, the possibility to set an initial weight, the strategy to delete outliers, the type of spatio-temporal filter, and the solver.

The regularization coefficient has the greatest impact on the TICOI solutions. To illustrate this, we compute the median RMSE and KGE for the velocity magnitude at the six Yukon GNSS stations (Fig. B1a). Both RMSE and KGE improve drastically, by factors of 3 and 8 respectively, when the regularization coefficient increases from 0.1 to the optimal value of 100. This optimal coefficient corresponds to 1.1min(RMSE) and 0.9max(KGE).

Notably, the RMSE increases slightly and the KGE decreases slightly (by about 5%) when the coefficient is further increased to 10 000, indicating the relative stability of the so-

lution even with large regularization coefficients in general scenarios. However, for surge-type glaciers, using a coefficient greater than 1000 increases the risk of over-smoothing the solution (Fig. B2).

The comparison between TICOI results and GNSS data helps to identify an optimal regularization coefficient of 100. In many cases, though, GNSS observations are unavailable to optimize the regularization coefficient. In such cases, we suggest that the optimal coefficient can also be determined as the approximate point of inflection of the VVC curve, defined in Eq. (9). This method similarly yields an optimal value of 100 for the coefficient. Thus, the VVC approach provides a reliable method for selecting an optimal regularization coefficient (Fig. B1b).

We analyzed the sensitivity of the results to various options in Appendix B. Our analysis shows that the choice of weight and solver has no impact on the results for the tested cases. The selection of spatio-temporal filters introduces a standard deviation in RMSE of approximately 2.6 m yr⁻¹ on average between filters (i.e. 8% of the averaged RMSE). For non-surge glaciers, using MA filters reduces the RMSE by up to 2 m yr⁻¹ compared to a strategy that does not filter outliers. This is largely because the Tukey biweight function, which is applied, already acts as an outlier filter.

## 4.4 Uncertainty of the final product

To provide insights into the uncertainty of TICOI velocity time series, we first use the VVC per pixel. This highlights areas where the direction of TICOI results has a poor temporal coherence. In Fig. 5 we can see that the VVC is lower near the terminus, the borders and the upper parts of the glaciers. We also notice low VVC values (around 0.7) just before the confluence between the southern arm and the main trunk of Lowell glacier. This corresponds to the approximate position of the equilibrium-line altitude, therefore leading to strong
<table><caption>**Table 1.** Comparison of the RMSE (in m yr⁻¹, left) and KGE (unitless, right) of the velocity magnitude between GNSS and: (1) TICOI time series; (2) image-pairs with a temporal baseline lower than 180 d, and (3) a moving 30 d median applied on image-pairs with a temporal baseline lower than 180 d. The lowest values are in bold. The surface velocity measurements used are the ITS_LIVE data in the Yukon, and IGE data in Greenland.</caption><thead><tr><th rowspan="2"></th><th colspan="3">RMSE</th><th colspan="3">KGE</th></tr><tr><th>TICOI</th><th>Image-pairs</th><th>Moving median</th><th>TICOI</th><th>Image-pairs</th><th>Moving median</th></tr></thead><tbody><tr><td>Kask L</td><td><strong>21.3</strong></td><td>68.4</td><td>28.3</td><td><strong>0.84</strong></td><td>-0.27</td><td>0.71</td></tr><tr><td>Kask M</td><td><strong>28.4</strong></td><td>62.5</td><td>31.8</td><td><strong>0.78</strong></td><td>-0.14</td><td>0.72</td></tr><tr><td>Kask U</td><td><strong>43.1</strong></td><td>90.5</td><td>55.5</td><td><strong>0.08</strong></td><td>-1.79</td><td>-0.36</td></tr><tr><td>Lowell L</td><td><strong>58.8</strong></td><td>64.5</td><td>308.7</td><td><strong>0.92</strong></td><td>0.86</td><td>0.5</td></tr><tr><td>Lowell M</td><td><strong>44.7</strong></td><td>60.2</td><td>155.7</td><td>0.93</td><td><strong>0.94</strong></td><td>0.66</td></tr><tr><td>Lowell U</td><td>78.7</td><td>119.9</td><td><strong>53.7</strong></td><td>0.49</td><td>0.53</td><td><strong>0.74</strong></td></tr><tr><td>Iss</td><td><strong>12.5</strong></td><td>28.8</td><td>28.2</td><td><strong>0.86</strong></td><td>0.63</td><td>0.59</td></tr></tbody></table>

changes in surface state (between wet snow, dry snow and bare ice) (NASA Earth Observatory, 2018) that are difficult to tackle. Besides, this is the only part of the glacier that faces north. It is therefore more susceptible to changes in shadows and illumination. (Lacroix et al., 2019). Note that real changes in velocity direction are expected over time in areas of variable flow, such as near confluences, glacier edges or terminus. It could cause a low VVC which is not related to measurement error.

Therefore, to estimate the uncertainty of each TICOI retrieved velocity, in time and space, we propagate the covariance matrix as described in Sect. 2.4. We tested this approach on simulated data with correlated noise described in Appendix C, for different percentage of image-pair velocities and different noise levels. In controlled conditions, the 95 % confidence interval includes both estimated and true velocity (i.e. are correct) for more than 95 % of the estimation, in cases of a low percentage of data and a low noise level, where the confidence interval tends to be slightly underestimated (Fig. B4). This synthetic case provides confidence of the validity of the theoretical framework.

In order to test the validity of the uncertainty calculation, we calculate the 95 % confidence intervals for the real dataset of Kask L (Fig. 6). Confidence intervals are larger before the launch of Sentinel-2 in 2015/2016, when the number and quality of data were lower. They are also larger in wintertime for the same reason. However, only 48 % of the confidence intervals include both the estimated and the GNSS velocities (i.e. are correct), which is much below the expected 95 %. On average, over the six GNSS stations, the percentage is 27 %. The confidence intervals fail to include estimated and true velocity, especially when the number of image-pair velocities is low, and the solution is less strongly constrained. This is a limitation of our approach to calculate confidence intervals that is discussed in Sect. 5.2.

## 4.5 Application to different glacier dynamics

In this section, we apply TICOI to pixels sampled regularly along the centerlines of Lowell and Kaskawulsh glaciers (Fig. 7). On Lowell glacier, we observe an upward propagation of the surge in 2021–2022 (Fig. 7a). This surge was first reported by Van Wychen et al. (2023) using Radarsat Constellation Mission data. However, since this data was only available from winter 2022, the start of the surge was poorly defined. From TICOI results, it is clear that there is a strong positive velocity anomaly in June 2021, with maximum anomalies of 400 m yr⁻¹ at around 2 km from the terminus (Fig. 7a). Velocity anomalies range between 100 and 400 m yr⁻¹ from 1 to 7 km upglacier from the terminus. They remain positive up to 27 km from the terminus, at a location which corresponds to the confluence between the northern and southern arm of Lowell. This positive anomaly remains approximately stable until December 2021 (with an average monthly acceleration of 15 m yr⁻²) except in the last kilometers before the terminus, where the average acceleration was of 1500 m yr⁻². The entire glacier starts to accelerate in December 2021, with an acceleration front propagating upglacier. Anomalies > 300 m yr⁻¹ are recorded up to 27 km from the terminus, just one month later, in January 2022. However, this acceleration front seems to propagate at a slower rate in the southern arm of Lowell glacier, even if this observation has to be interpreted cautiously regarding the low VVC values obtained in this area (Fig. 5). In October 2022, anomalies > 300 m yr⁻¹ are recorded up to at least 47 km from the terminus. In addition, we observe slight positive anomalies occurring each spring before the surge. The intensity of these anomalies increases from 2016 until the start of the surge, especially near the terminus, with annual maximum velocity magnitudes raising from 90 m yr⁻¹ between 2014–2016 to 275 m yr⁻¹ in 2020, at a rate of 9, 33, 70 and 70 m yr⁻² in 2017, 2018, 2019 and 2020 (Fig. 8a). The annual maximum over the entire glacier rises later in time (Fig. 8a). The surge event ends near the terminus in winter
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Spatial variability of VVC, overlaid with the glacier outlines in white according to the RGI v.7 (RGI 7.0 Consortium, 2023), over the same areas as in Fig. 2. (a) and (b) represent Lowell and Kaskawulsh glaciers, respectively. The coordinates of the grid correspond to the projection EPSG:4326.</figcaption></figure>

2022, while anomalies remain high in the upper part of the glacier.

On Kaskawulsh glacier, there is a long-term trend towards decreasing velocities (Fig. 7b), particularly over the lower part of the ablation area, which is likely related to the thinning of the glacier (Main et al., 2023; Dehecq et al., 2019). In the 0 to 5 km section upglacier from the terminus this decrease in velocity between 2015 and 2018 is particularly marked, from 170–180 to 120 m yr⁻¹ (Fig. 8b). This is consistent with the results of Main et al. (2023), who found a large change in velocity near the terminus after drainage of the adjacent proglacial Slims Lake in 2016. We also observe a marked seasonal velocity increase every spring (around April), with maximum of yearly anomalies between 30 and 70 m yr⁻¹. The anomalies propagate up glacier, but the interpretation is more difficult above 30 km from the terminus because the signal becomes noisier (Fig. 7b). The high temporal resolution also reveals very abrupt positive anomalies in March at a position situated from 6 to 9 km from the ter-

minus. This is especially the case for the year 2019, but it can be also noticed in 2017, 2020, 2021, and 2022. This phenomenon is also visible in the GNSS time series (Fig. 6).

## 4.6 Estimating monthly velocities from annual velocities

The TICOI method can be used to study sub-annual velocities even if long temporal baseline velocities only are available. This is important because in slow moving areas, velocities quantified with short temporal baseline mainly show low signal-to-noise ratio, and are frequently filtered out in open data datasets such as ITS_LIVE (Millan et al., 2019; Gardner et al., 2018). For these reasons, previous research chose to perform image correlation on images separated by long time intervals only, for example between 330 and 400 d (Hallas et al., 2023), privileging accuracy over temporal resolution. This strategy is adapted when assessing multi-annual or decadal velocity trends, but hinders the seasonal varia-
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) 95% confidence interval for velocities on the point Kask L. The blue dots and bar correspond to the central date and temporal baselines of the image-pair velocities. GNSS data are represented by orange crosses. Light magenta intervals represent the 95% confidence intervals defined in Sect. 2.4. The colors from light red to dark red correspond to the number of image-pair velocities used to constrain the velocity estimations. The upper and lower limits of the y-axis are defined according to the average ± the standard deviation of the image-pair velocities. (b) TICOI velocity magnitude as a function of the GNSS velocity magnitude. Vertical grey bars correspond to the confidence intervals, which should intersect the red line 1:1 (i.e., encompass the true velocity value) if they are not underestimated. Underestimated confidence intervals are displayed in red, correct one are represented in grey.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Spatio-temporal evolution of monthly velocity anomalies with respect to the averaged velocity magnitude over the period (as defined in Dehecq et al., 2019) over the centerline of (a) Lowell glacier, and (b) Kaskawulsh glacier, plotted as distance from the terminus. The centerline is represented in black on Fig. 2. The y-axis correspond to the central date of each monthly velocity.</figcaption></figure>

tions and rapid summer velocity changes (Fig. 9). In the case of Issunguata Sermia glacier, satellite observations show an average velocity of ~125 m yr⁻¹ with slight variations from year to year (Fig. 9). The application of TICOI to this dataset provides monthly velocities that match very well with the GNSS data after 2016, when Sentinel-2 images become available (Sentinel-2A was launched in June 2015), result-

ing in an increasing number of image-pair velocities (red and dark red in Fig. 9). The RMSE between GNSS and TICOI time series is 25.8 m yr⁻¹ between 2014 and 2017, and 17.7 m yr⁻¹ between 2016 and 2017, with stronger errors in winter time when optical images are impacted by night and clouds. Hence, TICOI can retrieve monthly velocities using only image-pair velocities with long temporal baselines. It
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Average yearly maximal velocity magnitude over the entire glacier in blue, and near the terminus in red, for (a) Lowell glacier, (b) Kaskawulsh glacier. Note the logarithm scale in (a).</figcaption></figure>

takes advantage of the temporal closure that relies on redundancy of annual velocities. However, it still requires a sufficient amount of observations (> 500) to obtain a reliable time series (Figs. 9 and C1). This is also illustrated over the Lower part of Kaskawulsh glacier using ITS_LIVE dataset (Fig. C1).

# 5 Discussion

## 5.1 Fusion of velocity measurements from different processing chains

Datasets from several processing chains can be included as input in TICOI, with the datasets reprojected on the same grid as explained in Sect. 3.3.2. Fusing different datasets, like the ones from IGE and ITS_LIVE, can particularly improve the signal-to-noise ratio in the upper parts of glaciers (Table 2). For example, this provides an improvement in RMSE of 11% for the upper part of Lowell glacier, and 32% for the upper part of Kaskawulsh glacier. Over these two areas the surface is snow-covered most of the year, which produces poor results from image correlation algorithms. Fusing velocity results from different processing chains takes advantage of the different sets of correlation parameters used and strategies to delete outliers. However, we do not see improvements in the middle and lower parts of the glaciers when applying this technique; the RMSE even slightly increases (Table 2). Considering the increasing computation time when including additional dataset, this technique may only be suited for very noisy areas or areas with a low number of image-pair velocities.

**Table 2.** Comparison of the RMSE of the velocity magnitude (in m yr⁻¹) between GNSS and TICOI time series, with datasets from (1) ITS_LIVE and (2) IGE & ITS_LIVE over upper and middle GNSS stations of Lowell and Kaskawulsh glaciers. We consider the time-period which is common between the two datasets, i.e. 2016–2022. The lowest values are in bold.

<table><thead><tr><th></th><th>ITS_LIVE</th><th>IGE & ITS_LIVE</th></tr></thead><tbody><tr><td>Lowell U</td><td>46.2</td><td><strong>41.5</strong></td></tr><tr><td>Kask U</td><td>46.7</td><td><strong>35.3</strong></td></tr><tr><td>Kask M</td><td><strong>26.0</strong></td><td>28.5</td></tr><tr><td>Lowell M</td><td><strong>30.1</strong></td><td>33.7</td></tr></tbody></table>

## 5.2 Uncertainty

Despite providing satisfying confidence intervals on synthetic data, our framework to calculate TICOI confidence intervals underestimates uncertainties on real data when compared to GNSS measurements (see Sect. 4.4). Here, we discuss the different limitations in our approach that could explain these discrepancies. First, the error of the input image-pair velocities may be underestimated. These errors are computed either using stable areas, which may not be representative of the glacier texture (Zheng et al., 2023; Altena et al., 2022, 2019) or by assuming a 1/10 pixel error, which corresponds to an error of 1 m for Sentinel-2 and 3 m for Landsat-8. To estimate the underestimation of the errors provided with image-pair velocities, we consider the true error as the difference between ITS_LIVE and GNSS data. The averaged RMSE between the errors provided by the ITS_LIVE dataset and the true errors are 9 m and 7 m for the East/West and North/South components, respectively. The difference between the average true error and ITS_LIVE error is of 1.8 and 2.6 m, highlighting a clear underestimation in ITS_LIVE errors. Moreover, there are no strong temporal correlations between the provided errors and the true errors, with Pearson's coefficient of 0.27 on average for both com-
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Example of monthly velocity time series (in red) retrieved using annual velocities from Halas et al. (2023) (in blue) on the Issunguata Sermia glacier. GNSS data are represented by orange crosses. The colors from light red to dark red correspond to the number of image-pair velocities used to constrain the velocity estimations.</figcaption></figure>

ponents. Second, the underestimation of our confidence in-
tervals could be caused by biases in the image-pair velocities,
for example due to shadows or seasonal illumination changes
(Lacroix et al., 2019). The errors in the ITS_LIVE dataset are
based on the standard error in component velocities relative
to stable surface velocity; they characterize random errors.
Therefore, our confidence intervals only account for random
errors and not systematic biases. Third, the real errors may
have a stronger correlation in time than what has been simu-
lated in Sect. 4.4, for instance due to seasonal source of errors
(shadows, surface changes), since we only take into account
the correlation of errors between displacement with common
acquisition dates. In case of highly correlated errors, the *a priori*
covariance matrix cannot be diagonal anymore, lead-
ing to much higher computational cost. The cost of explicitly
computing the inverse of the error covariance matrix is pro-
portional to *n* (the number of image-pair velocities) if the
matrix is diagonal and to *n*³ in the general case (Ruggiero
et al., 2016).

To address this issue, a potential solution involves scal-
ing the confidence interval by a specific correction factor.
This factor could be a function of the VVC, as the VVC
serves as an effective proxy for the relative uncertainty be-
tween pixels. Indeed, the RMSE decreases linearly with the
VVC, until reaching a vertical asymptote close to 1, which is
the VVC maximum value (Fig. B5a). To estimate this correc-
tion factor, we determine the 95th quantile of the theoretical
95 % confidence intervals divided by the true errors, for each
GNSS stations (Fig. B5b). We then perform a linear regres-
sion on the resulting seven points. Although the RMSE of
the linear regression is around 3 and the sample size is small,

this method provides an empirical correction factor. By mul-
tiplying the confidence intervals by this factor, the confidence
intervals contain both estimated and GNSS data with an av-
erage percentage of 86 % (against 27 % without a correction
factor). Additionally, the reliability of the confidence inter-
vals derived with the correction factor appears to improve
with the number of observations used. For example, when
selecting only TICOI estimations with more than 500 obser-
vations, the confidence intervals encompass both estimated
and GNSS data with an average coverage of 91 %.

Another strategy could be to augment the observation vec-
tor with the first- and second-order spatial derivatives of the
original observations, as described in Ruggiero et al. (2016);
Brankart et al. (2009). However, this requires proper char-
acterization of the spatial and temporal correlation of errors
of surface velocities, which could be the scope of future re-
search.

With the current state of knowledge in velocity errors, we
recommend relying on the VVC and number of contributed
image-pair velocities. The VVC is a quality metric which
characterize random errors, by analazing the temporal coherence
of the direction. The number of contributed image-pair
velocities indicates the robustness of the TICOI estimation:
the solution seems to be correctly constrained by the obser-
vations when this number is higher than 100.

Note, that the TICOI package also offers the possibil-
ity to compute the Normalized Median Absolute Deviation
(NMAD) over stable areas. This is a widely used and ro-
bust metric for characterizing random errors in glacier ve-
locity fields (Dehecq et al., 2015). However, as previously
demonstrated for ITS_LIVE scene-pair velocities, such er-
errors are often underestimated due to the differences in texture between glacier surfaces and stable ground. Moreover, the NMAD over stable ground do not capture the spatial variability in errors because they provide only one value for the entire scene at a given time. For example, the RMSE between TICOI and the GNSS is about 43 m yr⁻¹ on the upper of Kaskawulsh glacier, which is 20 m yr⁻¹ below the maximal NMAD obtained in the area (Fig. A3b). To enhance the flexibility of the package, the application of the along-flow shear strain rate proposed by Zheng et al. (2023) is also implemented. It provides insight into the smoothness of the velocity solution.

## 5.3 Large scale application

Here, we discuss the possibility of applying TICOI at a large scale. We have shown that it can be applied to all kinds of glacier dynamics because it does not include any *a priori* information about the glacier behavior, unlike a wide range of post-processing approaches (Riel et al., 2021; Samsonov et al., 2021; Greene et al., 2020; López-Quiroz et al., 2009). This flexibility also allows for the detection of unexpected events and trends, such as the annual acceleration in March over Kaskawulsh glacier. However, data-driven approaches may encounter limitations when data density is very low. In this case, *a priori* information, if available, may help to constrain the time series. The regularization term can be modified to include a model, similar to López-Quiroz et al. (2009). An example is given in Eq. (C1). By doing so, the inversion solves both the temporal closure and a parametric regression problem. For more flexibility, it may also be possible to use a dictionary of functions (Riel et al., 2021; Hetland et al., 2012).

TICOI can be applied in a nearly automatic way, at a reasonable computational cost. We proposed a general strategy to automatically select the regularization coefficient. Moreover, the computation time of the TICOI processing chain including loading, pre-processing, inversion and saving is about 0.1 s per pixel by using 32 CPUs (on a Intel(R) Xeon(R) Gold 6426Y with a processing rate of 3.3 GHz), for datasets containing 80 000 layers in time (corresponding to the period 2013–2024). This means that the processing over a region of 100 km × 100 km requires around 19 h using 32 CPUs only. Note that the processing time per pixel scales with the size of the data in time, and that the 2013–2024 period has high density of measurements compared to previous years. This computation time remains affordable at the regional scale, and even at the global scale with a large number of CPUs. The computation time could be further reduced, for instance, by taking advantage of GPUs, or by reducing the number of input data with a stricter outlier filter.

With these two points in mind, it appears relevant to apply TICOI at a large scale. First, it reduces the size of the data in comparison with raw image-pair velocities by removing redundant information. On average, for our study sites, the size

of the data is reduced by a factor 100. Second, it produces regular velocity time series, with relevant quality indicators. This would make the integration of sub-annual velocities in numerical models much more affordable.

## 6 Conclusions

To derive sub-annual velocity variations over glaciers, we propose an operational Python package, called TICOI, based on the temporal closure principle. This package fuses multi-temporal and multi-sensor image-pair velocities computed by different processing chains to generate regular velocity time series (i.e., sampled at regular time steps).

We improved previous methods based on two strategies: (1) a regularization term based on a coarse initial estimate, which enhances the resilience of the temporal closure inversion to abrupt non-linear changes; and (2) an iteratively reweighted approach, which automatically detects temporal decorrelation. TICOI is entirely data-driven (i.e., it does not require strong *a priori* information on the glacier dynamic), making it applicable to any glacier dynamics. The validation of TICOI results using GNSS data highlights an improvement in RMSE and KGE of around 50% and 0.4, respectively, compared to both the raw image-pair velocities and a rolling median. Furthermore, it is able to retrieve monthly velocities using annual image-pair velocities only when there is sufficient temporal redundancy in the dataset. This could be especially useful for slow moving areas, where annual image-pair velocities may be of better quality than image-pair velocities with short temporal baselines. Moreover, TICOI can be used to combine datasets from different processing chains. This has the potential to reduce the uncertainty in the upper part of glaciers, such as in the accumulation area, where image-pair velocities are more noisy.

We recommend using three criteria to assess the quality of the retrieved velocity series: (1) the VVC, a proxy of the temporal coherence of the direction; (2) the number of contributing image-pair velocities, and (3) a 95% confidence interval derived from the *a posteriori* covariance matrix. The application of TICOI provides velocity time series with an unprecedented temporal resolution. On Lowell glacier (Yukon, Canada), we are able to observe that summer velocities near the terminus started to increase five years before the surge. On Kaskawulsh glacier (Yukon, Canada), we are able to resolve velocity peaks in March in a very localized part of the lower ablation area.

Finally, the TICOI workflow offers reasonable computation time for application at the regional to global scale (0.1 s per pixel for large dataset with 80 000 layers in time on 32 CPUs). The code is open source and can be applied to any datasets and regions. This paves the way for the integration of a wide range of image-pair velocities and the production of standardized post-processed sub-annual velocity time series.
## Appendix A: Demonstration of the posteriori covariance matrix formula

We demonstrate how the posteriori covariance matrix formula (Eq. 11) has been derived. First, the analytic solution of a least square problem with a Tikhonov regularization is:

$$ \hat{X} = N^{-1} A^T W Y \quad (A1) $$

with $N = A^T W A + \lambda \Gamma^T \Gamma$

The posteriori covariance matrix is:

$$ \Sigma_{\hat{X}} = \left[ \frac{\partial \hat{X}}{\partial Y} \right] \Sigma_{\hat{Y}} \left[ \frac{\partial \hat{X}}{\partial Y} \right]^T \quad (A2) $$

It corresponds to:

$$ \Sigma_{\hat{X}} = N^{-1} A^T W \Sigma_{\hat{Y}} W A N^{-1} \quad (A3) $$

<figure><img src="image_12.png" /><figcaption>**Figure A1.** Distribution of the errors estimated by comparing measured displacements from remote sensing images and GNSS displacements. The left column shows errors in East/West ($D_x$) and North/South ($D_y$) displacements according to the temporal baseline. The right column shows the distribution of this error. Skewness is a measure of the symmetry of a distribution, a value of 0 indicates a symmetric distribution. Kurtosis refers to the degree of “tailedness” of a distribution relative to a normal distribution. Strong kurtosis (> 3) reveals heavy-tailed distribution.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

**Figure A2.** Scatterplot of GNSS velocity magnitude and 30 d rolling median applied to velocity magnitude observations, (a) with a temporal baseline lower than 180 d and (b) with every temporal baseline. The RMSE is better while using all temporal baselines, but there is a clear underestimation for velocities larger than 180 m yr⁻¹.

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

**Figure A3.** Normalized Median Absolute Deviation of the velocity magnitude over stable areas represented in Fig. 3 using TICOI without (a) and with (b) an automatic detection of temporal decorrelation. Median of the velocity magnitude over stable areas using TICOI without (c) and with (d) an automatic detection of temporal decorrelation, which hilight an oversmoothing of the solution.
# Appendix B: Sensitivity analysis and automated choice of the hyperparameters

## B1 Regularization coefficient

The regularization coefficient is the parameter that has the greatest impact on the results. We proposed using the VVC to strike a balance between enhancing the signal-to-noise ratio and avoiding oversmoothing (Fig. B1).

## B2 Spatio-temporal filter

We compared the performance of five types of spatio-temporal filters: Savitzky-Golay, Gaussian, Locally Weighted Scatterplot Smoothing (LOWESS), and median. The choice of filter results in a standard deviation in RMSE of about 2.6 m yr⁻¹ on average (i.e., 8% of the averaged RMSE), except for the GNSS station Lowell L, where the LOWESS filter produces an RMSE approximately 170 m yr⁻¹ higher than the other filters (Table B1).

LOWESS is a non-parametric moving regression that fits a model to the k nearest points (Derkacheva et al., 2020), which tends to over-smooth data during periods with low observation density. For example, the LOWESS solution flattens the surge peak of Lowell L, due to the limited number of observations available during that time (Fig. B3a).

We note that both the LOWESS and median filters can provide slightly better results for non-surge type glaciers, with improvements ranging from 0.8 to 4.5 m yr⁻¹ (i.e., 3% to 10%). However, they can also lead to over-smoothing (Fig. B3b) and LOWESS require 1.5 times more computational time. Therefore, we recommend using the Savitzky-Golay filter, which offers a good balance between computational efficiency and accuracy in general scenarios.

## B3 Solver

We compare four different solvers: the Least Square solver (LS), LSMR, LSMR with an initialisation and LSQR. The RMSE between GNSS and TICOI time-series are really similar among GNSS stations. However, the computational time of the LS is twice larger than for the other solvers. Therefore, we recommend using LSMR, LSMR\_ini or LSQR for a question of computational time.

## B4 Strategy to delete outliers

We compared two strategies for outlier removal: the Modified Z-score (Mzscore), which filters out values more than 3.5 NMAD away from the median of the entire time period, and the Median Angle (MA), which removes observations deviating by more than 45° from the median vector (Charrier et al., 2022a). We strongly advise against using the Mzscore for surge-type glaciers. For non-surge glaciers, the improvement is at most 2 m yr⁻¹ compared to a strategy that does not filter outliers, largely because the Tukey biweight function, which is used, already acts as an outlier filter. Nevertheless, we recommend applying at least the MA filter to reduce the number of observations input into TICOI.

**Table B1.** Comparison of RMSE of the velocity magnitude between GNSS and TICOI time-series in m yr⁻¹ for different spatio-temporal filters (columns) and GNSS stations (rows). The averaged computational time between all the GNSS stations are given in the last row in s. The filter svagol corresponds to the Savitzky-Golay filter. The LOWESS filter corresponds to the statsmodel.nonparametric implementation and the Savitzky-Golay to scipy.signal implementation. The Standard Deviation (SD) at each GNSS for the different filters is given in the last column. The fraction used for the LOWESS filter is 60/N with N the number of observations over the period. The temporal window of the Savitzky-Golay, median and gaussian filters are of 90 d. We note that the better performance obtained by the median and LOWESS filters on Lowell M is mainly due to the absence of GNSS data during the maximum of the surge.

<table><thead><tr><th></th><th>savgol</th><th>gaussian</th><th>median</th><th>LOWESS</th><th>SD</th></tr></thead><tbody><tr><td>Lowell L</td><td><strong>61.4</strong></td><td>76.7</td><td>72.9</td><td>229.6</td><td>70.7</td></tr><tr><td>Lowell M</td><td>43.5</td><td>46.1</td><td>36.5</td><td><strong>31.7</strong></td><td>6.0</td></tr><tr><td>Lowell U</td><td><strong>85.7</strong></td><td>86.0</td><td>88.4</td><td>87.4</td><td>1.0</td></tr><tr><td>Kask L</td><td>22.3</td><td>23.8</td><td>20.1</td><td><strong>20.0</strong></td><td>1.8</td></tr><tr><td>Kask M</td><td>26.0</td><td>27</td><td><strong>25.2</strong></td><td>25.4</td><td>0.8</td></tr><tr><td>Kask U</td><td>42.2</td><td>46.1</td><td><strong>37.7</strong></td><td>39.7</td><td>3.6</td></tr><tr><td>Computational time [s]</td><td>67.8</td><td>67.2</td><td>69.4</td><td>101.6</td><td></td></tr></tbody></table>
**Table B2.** Comparison of RMSE of the velocity magnitude between GNSS and TICOI time-series in m yr⁻¹ for different solvers (columns) and GNSS stations (rows). The average computational time between all the GNSS stations are given in the last row in s. The solver Least Square (LS) corresponds to the function `lstsq` of ` numpy.lintalg`, the solvers `LSMR` and `LSQR` are respectively the functions `lsmr` and `lsqr` of `scipy.lintalg` and `scipysparse`. `LSMR_ini` corresponds to the solver `LSMR` with an initialisation, corresponding to the spatio-temporal filtering observations $X_0$ defined in Appendix 2.2.2.

<table><thead><tr><th></th><th>LSMR</th><th>LSMR_ini</th><th>LSQR</th><th>LS</th></tr></thead><tbody><tr><td>Lowell L</td><td>61.16</td><td>61.39</td><td>61.11</td><td>61.13</td></tr><tr><td>Lowell M</td><td>46.0</td><td>45.99</td><td>46.0</td><td>46.0</td></tr><tr><td>Lowell U</td><td>86.24</td><td>86.24</td><td>86.24</td><td>86.25</td></tr><tr><td>Kask L</td><td>20.9</td><td>20.89</td><td>20.9</td><td>20.91</td></tr><tr><td>Kask M</td><td>25.03</td><td>25.03</td><td>25.03</td><td>25.04</td></tr><tr><td>Kask U</td><td>40.12</td><td>40.15</td><td>40.13</td><td>40.14</td></tr><tr><td>Comptutional time [s]</td><td><strong>40.8</strong></td><td><strong>48.8</strong></td><td><strong>43.0</strong></td><td><strong>85.2</strong></td></tr></tbody></table>

**Table B3.** Comparison of RMSE of the velocity magnitude between GNSS and TICOI time-series in m yr⁻¹ for different strategy to delete outliers (columns) and GNSS stations (rows). The minimal values are displayed in blod.

<table><thead><tr><th></th><th>median_angle</th><th>mz_score</th><th>no_delete_outliers</th></tr></thead><tbody><tr><td>Lowell L</td><td><strong>58.79</strong></td><td>711.44</td><td>59.15</td></tr><tr><td>Lowell M</td><td><strong>44.71</strong></td><td>511.52</td><td>44.8</td></tr><tr><td>Lowell U</td><td><strong>78.72</strong></td><td>96.74</td><td>81.72</td></tr><tr><td>Kask L</td><td><strong>21.27</strong></td><td>23.2</td><td>21.93</td></tr><tr><td>Kask M</td><td>28.37</td><td><strong>25.46</strong></td><td>27.36</td></tr><tr><td>Kask U</td><td>43.15</td><td><strong>41.25</strong></td><td>43.2</td></tr><tr><td>Iss</td><td><strong>12.53</strong></td><td>15.1</td><td>15.4</td></tr></tbody></table>

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

**Figure B1.** (a) Evolution of the RMSE and KGE of TICOI results with respect to the six Yukon GNSS stations according to the regularization coefficient. (b) Evolution of the VVC computed over the six red squares shown in Fig. 2 as a function of the regularization coefficient. The optimal coefficient corresponds to max(VVC) - 0.05·(max(VVC) - min(VVC)). The optimal coefficient found is 100 using both approaches (a) and (b), which confirms the value of using VVC for selecting the regularization coefficient.
<figure><img src="image_16.png" /><figcaption>**Figure B2.** Evolution of the RMSE and KGE of TICOI results with respect to GNSS data for station Lowell L, where a surge occurs. When the coefficient increases the acceleration of the TICOI estimations tend to be close the initial guess of acceleration, which in this case slightly oversmooths the peak of the surge. This is why the RMSE and KGE reach a plateau after a coefficient of around 5000.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure B3.** (a) Visual comparison of the spatio-temporal filters for the point Lowell M. (b) Visual comparison of the spatio-temporal filters for the point Lowell L.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure B4.** Percentage of the estimated 95 % confidence intervals that include both the estimated and the true displacement values (i.e. the valid confidence intervals), using simulated data described in Appendix C. The displacement noise factors correspond to a factor of the simulated instantaneous velocity amplitude. A factor of 1, 0.1 and 0.02 are equivalent to a standard deviation of 4.6, 0.4 and 0.09 m, respectively. The percentage of data corresponds to the percentage of simulated image-pair displacements in comparison with the total number of possible image-pair displacements. The value of each configuration of displacement noise factor and percentage of image-pair displacements relies on 50 experiments.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure B5.** (a) Scatter plot of VVC and RMSE over the seven stations analyzed in this study, from 2013–2022. The RMSE corresponds to the comparison between TICOI results and GNSS data, while the VVC is an indicator of the temporal coherence of the directions of the TICOI results. (b) Correction factor that should be applied to the theoretical 95 % confidence intervals so that they include both the TICOI estimate and the GNSS data.</figcaption></figure>
## Appendix C: Simulated data

### C1 Synthetic instantaneous velocity and position time-series

The synthetic instantaneous velocity is taken as: $v(t) = a + b \sin(\frac{2\pi}{T}t) + c \cos(\frac{2\pi}{T}t)$ with $T = 365.25$ as in Greene et al. (2020). To make sure that the coefficients $a, b$ and $c$ represent well the data, instead of an arbitrary choice, these coefficients are estimated by an IRLS in-version by adding a regularization term corresponding to a displacement model with the coefficients $a, b$ and $c$ as parameters as performed in López-Quiroz et al. (2009). The corresponding system of equations is given in Eq. (C1). The system is solved for Sentinel-2 data on the point represented in blue in Charrier et al. (2022c) Fig. S1. The coefficients are found to be: $a = -0.49$, $b = -0.0788$ and $c = 0.018$.

$$ \begin{bmatrix} 1 & 1 & 0 & \dots & 0 & 0 & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & 1 & 0 & 0 & 0 & 0 & 0 \\ \lambda & 0 & 0 & \dots & 0 & -[\tau_{1+\Delta\tau} - \tau_1] & -\frac{T}{2\pi}[\sin(\frac{2\pi}{T}\tau_{1+\Delta\tau}) - \sin(\frac{2\pi}{T}\tau_1)] & -\frac{T}{2\pi}[\cos(\frac{2\pi}{T}\tau_{1+\Delta\tau}) - \cos(\frac{2\pi}{T}\tau_1)] & 1 & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \dots & \lambda & -[\tau_{1+n\Delta\tau} - \tau_{1+(n-1)\Delta\tau}] & -\frac{T}{2\pi}[\sin(\frac{2\pi}{T}\tau_{1+n\Delta\tau}) - \sin(\frac{2\pi}{T}\tau_{1+(n-1)\Delta\tau})] & -\frac{T}{2\pi}[\cos(\frac{2\pi}{T}\tau_{1+n\Delta\tau}) - \cos(\frac{2\pi}{T}\tau_{1+(n-1)\Delta\tau})] & 1 & 1 \end{bmatrix} $$

(C1)

where $T$ is the period of the sinusoidal signal. $a, b, c, d$ are the coefficients of the model.

The position time-series is defined as the integral of the instantaneous velocity.

### C2 Selection of acquisition dates

Then, we randomly select the acquisition dates in a list of dates ranging from the 1 January 2015 and 31 December 2020, every 5 d. By doing so, some dates between the 1 January 2015 and 31 December 2020 will not be included in the simulated dataset. It represents the effect of clouds: the pixels covered by clouds will be systematically rejected.

### C3 Noise

For each acquisition date, we add a Gaussian noise to the position value. This accounts for the fact that the noise depends mainly on the image texture (clouds, snow, crevasses, etc.). Therefore, the noise of each displacement is the sum of the noises of each date of the pair.
## C4 Image-pair velocity

We randomly select a temporal baseline between 5 and 400. Then, we compute image-pair velocity by taking the difference between the position at the second date of acquisition and the first date of acquisition.

<figure><img src="image_20.png" /><figcaption>**Figure C1.** Example of monthly velocity time series (in red) retrieved using annual velocities from ITS_LIVE image pair velocities, with a temporal baseline from 200 to 400 d (in blue) on the Lower part of the Kaskawulsh glacier (Kask L). GNSS data are represented by orange crosses. The colors from light red to dark red correspond to the number of image-pair velocities used to constrain the velocity estimations.</figcaption></figure>

*Code and data availability.* The TICOI package is available on github (https://github.com/ticoi/ticoi) and Zenodo (https://doi.org/10.5281/zenodo.17209282, Charrier et al., 2025). ITS_LIVE data are available on https://its-live.jpl.nasa.gov/, last access: September 2023.

*Author contributions.* LaC, AD, FB, and RM designed the study. LaC, LG, AD, and FB proposed the methodological improvements. LaC, LG and NL worked on the python package with advice from AD. LaC generated the velocity time series and analyzed the results with feedbacks from AD, FB and RM. LuC and CD provided processed GNSS data over Yukon, and helped with analysis of the results. LaC carried out the post-processing of Yukon GNSS data. NM provided processed GNSS data over western Greenland. LaC led the writing, and all authors contributed to it.

*Competing interests.* The contact author has declared that none of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Also, please note that this paper has not received English language copy-editing.

*Acknowledgements.* We acknowledge Antoine Rabatel for his help in designing the study and his relevant feedbacks, Ghislain Picard and Fabien Gillet-Chaulet for their advice concerning the code optimization and the regularization term. We acknowledge the use of ChatGPT and DeepL for assistance in rephrasing and proofreading some sentences in the manuscript.

*Financial support.* Laurane Charrier acknowledges support from the Centre National d'Etudes Spatiales for her postdoctoral fellowship. Amaury Dehecq and Laurane Charrier acknowledge support from the French Programme National de Télédétection Spatiale (PNTS). Lei Guo thanks to the China Scholarship Council (grant no. 202306370154) for his scholarship. Luke Copland and Christine Dow thank the Natural Sciences and Engineering Research Council of Canada, Polar Continental Shelf Program, Canada Foundation for Innovation, Ontario Research Fund, and New Frontiers Research Fund for funding to purchase and operate the Yukon GNSS stations,
and the staff of Kluane Lake Research Station and graduate students from University of Ottawa and University of Waterloo for assistance and support in the field.

*Review statement.* This paper was edited by Johannes J. Fürst and reviewed by Benjamin Wallis and Maximillian Van Wyk de Vries.

## References

Altena, B., Scambos, T., Fahnestock, M., and Kääb, A.: Extracting recent short-term glacier velocity evolution over southern Alaska and the Yukon from a large collection of Landsat data, The Cryosphere, 13, 795–814, https://doi.org/10.5194/tc-13-795-2019, 2019.

Altena, B., Kääb, A., and Wouters, B.: Correlation dispersion as a measure to better estimate uncertainty in remotely sensed glacier displacements, The Cryosphere, 16, 2285–2300, https://doi.org/10.5194/tc-16-2285-2022, 2022.

Anderton, P. W.: Structural glaciology of a glacier confluence, Kaskawulsh glacier, Yukon Territory, Canada, Tech. rep., Research Foundation and the Institute of Polar Studies, The Ohio State, 1973.

Arendt, A. A., Echelmeyer, K. A., Harrison, W. D., Lingle, C. S., and Valentine, V. B.: Rapid wastage of Alaska glaciers and their contribution to rising sea level, Science, 297, 382–386, 2002.

Bartholomew, I., Nienow, P., Mair, D., Hubbard, A., King, M. A., and Sole, A.: Seasonal evolution of subglacial drainage and acceleration in a Greenland outlet glacier, Nat. Geosci., 3, 408–411, 2010.

Bartholomew, I., Nienow, P., Sole, A., Mair, D., Cowton, T., and King, M. A.: Short-term variability in Greenland Ice Sheet motion forced by time-varying meltwater drainage: Implications for the relationship between subglacial drainage system behavior and ice velocity, J. Geophys. Res.-Earth, 117, https://doi.org/10.1029/2011JF002220, 2012.

Beaud, F., Aati, S., Delaney, I., Adhikari, S., and Avouac, J.-P.: Surge dynamics of Shisper Glacier revealed by time-series correlation of optical satellite images and their utility to substantiate a generalized sliding law, The Cryosphere, 16, 3123–3148, https://doi.org/10.5194/tc-16-3123-2022, 2022.

Berardino, P., Fornaro, G., Lanari, R., and Sansosti, E.: A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms, IEEE Transactions on geoscience and remote sensing, 40, 2375–2383, https://doi.org/10.1109/TGRS.2002.803792, 2002.

Bevington, A. and Copland, L.: Characteristics of the last five surges of Lowell Glacier, Yukon, Canada, since 1948, J. Glaciol., 60, 113–123, 2014.

Bolibar, J., Sapienza, F., Maussion, F., Lguensat, R., Wouters, B., and Pérez, F.: Universal differential equations for glacier ice flow modelling, Geosci. Model Dev., 16, 6671–6687, https://doi.org/10.5194/gmd-16-6671-2023, 2023.

Bontemps, N., Lacroix, P., and Doin, M.-P.: Inversion of deformation fields time-series from optical images, and application to the long term kinematics of slow-moving landslides in Peru, Remote Sens. Environ., 210, 144–158, 2018.

Brankart, J.-M., Ubelmann, C., Testut, C.-E., Cosme, E., Brasseur, P., and Verron, J.: Efficient parameterization of the observation error covariance matrix for square root or ensemble Kalman filters: application to ocean altimetry, Mon. Weather Rev., 137, 1908–1927, 2009.

Breunig, M. M., Kriegel, H.-P., Ng, R. T., and Sander, J.: LOF: identifying density-based local outliers, in: Proceedings of the 2000 ACM SIGMOD international conference on Management of data, 93–104, https://doi.org/10.1145/342009.335388, 2000.

Charrier, L., Yan, Y., Colin Koeniguer, E., Mouginot, J., Millan, R., and Trouvé, E.: FUSION OF MULTI-TEMPORAL AND MULTI-SENSOR ICE VELOCITY OBSERVATIONS, ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, V-3-2022, 311–318, https://doi.org/10.5194/isPRS-annals-V-3-2022-311-2022, 2022a.

Charrier, L., Yan, Y., Koeniguer, E. C., Leinss, S., and Trouve, E.: Extraction of velocity time series with an optimal temporal sampling from displacement observation networks, IEEE Transactions on Geoscience and Remote Sensing, 60, 1–10, https://doi.org/10.1109/TGRS.2021.3128289, 2022b.

Charrier, L., Yan, Y., Trouvé, E., Koeniguer, E. C., Mouginot, J., and Millan, R.: Fusion of Multitemporal Multi-sensor Velocities Using Temporal Closure of Fractions of Displacements, IEEE Geosci. Remote Sens. Lett., 19, 1–5, https://doi.org/10.1109/LGRS.2022.3227413, 2022c.

Charrier, L., Guo, L., and Dehecq, A.: TICOI Software (v.0.1.1). Zenodo [code], https://doi.org/10.5281/zenodo.17209295, 2025.

Choi, Y., Seroussi, H., Morlighem, M., Schlegel, N.-J., and Gardnei, A.: Impact of time-dependent data assimilation on ice flow model initialization and projections: a case study of Kjer Glacier, Greenland, The Cryosphere, 17, 5499–5517, https://doi.org/10.5194/tc-17-5499-2023, 2023.

Clarke, G. K.: A short and somewhat personal history of Yukon glacier studies in the twentieth century, Arctic, 1–21, 2014.

Clarke, G. K. C., Bushnell, V., and Ragle, R.: Geophysical measurements on the Kaskawulsh and Hubbard Glaciers, Arctic Institute of North America, 1967.

Copland, L., Sylvestre, T., Bishop, M. P., Schroder, J. F., Seong, Y. B., Owen, L. A., Bush, A., and Kamp, U.: Expanded and recently increased glacier surging in the Karakoram, Arct. Antarct. Alp. Res., 43, 503–516, 2011.

Davison, B. J., Sole, A. J., Livingstone, S. J., Cowton, T. R., and Nienow, P. W.: The influence of hydrology on the dynamics of land-terminating sectors of the Greenland Ice Sheet, Front. Earth Sci., p. 10, https://doi.org/10.3389/feart.2019.00010, 2019.

Dehecq, A., Gourmelen, N., and Trouvé, E.: Deriving large-scale glacier velocities from a complete satellite archive: Application to the Pamir–Karakoram–Himalaya, Remote Sens. Environ., 162, 55–66, 2015.

Dehecq, A., Gourmelen, N., Gardner, A. S., Brun, F., Goldberg, D., Nienow, P. W., Berthier, E., Vincent, C., Wagnon, P., and Trouvé, E.: Twenty-first century glacier slowdown driven by mass loss in High Mountain Asia, Nat. Geosci., 12, 22–27, 2019.

Derkacheva, A., Mouginot, J., Millan, R., Maier, N., and Gillet-Chaulet, F.: Data Reduction Using Statistical and Regression Approaches for Ice Velocity Derived by Landsat-8, Sentinel-1 and Sentinel-2, Remote Sens., 12, 1935, https://doi.org/10.3390/rs12121935, 2020.
Derkacheva, A., Gillet-Chaulet, F., Mouginot, J., Jager, E., Maier, N., and Cook, S.: Seasonal evolution of basal environment conditions of Russell sector, West Greenland, inverted from satellite observation of surface flow, The Cryosphere, 15, 5675–5704, https://doi.org/10.5194/tc-15-5675-2021, 2021.

Doin, M.-P., Lodge, F., Guillaaso, S., Jolivet, R., Lasserre, C., Ducret, G., Grandin, R., Pathier, E., and Pinel, V.: Presentation of the Small Baselin NSBAS Processing Chain on a Case Example: The Etan Deformation Monitoring from 2003 to 2010 Using Envisat Data, in: Fringe Symposium, HAL-02185213, 2011.

Doyle, S. H., Hubbard, A., Van De Wal, R. S., Box, J. E., Van As, D., Scharrer, K., Meierbachtol, T. W., Smeets, P. C., Harper, J. T., Johansson, E., et al.: Amplified melt and flow of the Greenland ice sheet driven by late-summer cyclonic rainfall, Nat. Geosci., 8, 647–653, 2015.

Flowers, G. E., Copland, L., and Schoof, C. G.: Contemporary glacier processes and global change: recent observations from Kaskawulsh Glacier and the Donjek Range, St. Elias Mountains, Arctic, 22–34, https://doi.org/10.14430/arctic4356, 2014.

Foy, N., Copland, L., Zdanowicz, C., Demuth, M., and Hopkinson, C.: Recent volume and area changes of Kaskawulsh Glacier, Yukon, Canada, J. Glaciol., 57, 515–525, 2011.

Gardner, A. S., Moholdt, G., Scambos, T., Fahnstock, M., Ligtenberg, S., van den Broeke, M., and Nilsson, J.: Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years, The Cryosphere, 12, 521–547, https://doi.org/10.5194/tc-12-521-2018, 2018.

Gardner, A. S., Fahnstock, M., and Scambos, T.: MEaSURES ITS_LIVE Landsat Image-Pair Glacier and Ice Sheet Surface Velocities, Version 1, https://doi.org/10.5067/IMR9D3PEI28U, 2022.

Gardner, A. S., Greene, C. A., Kennedy, J. H., Fahnstock, M. A., Liukis, M., López, L. A., Lei, Y., Scambos, T. A., and Dehecq, A.: ITS_LIVE global glacier velocity data in near-real time, The Cryosphere, 19, 3517–3533, https://doi.org/10.5194/tc-19-3517-2025, 2025.

Gavin, H. P.: Fitting Models to Data: Generalized Linear Least Squares and Error Analysis, https://apisemanticscholar.org/ CorpusID:145046314 (last access: 1 November 2024), 2023.

Gilbert, A., Leinss, S., Kargel, J., Kääb, A., Gascoin, S., Leonard, G., Berthier, E., Karki, A., and Yao, T.: Mechanisms leading to the 2016 giant twin glacier collapses, Aru Range, Tibet, The Cryosphere, 12, 2883–2900, https://doi.org/10.5194/tc-12-2883-2018, 2018.

Goldberg, D. N., Heimbach, P., Joughin, I., and Smith, B.: Committed retreat of Smith, Pope, and Kohler Glaciers over the next 30 years inferred by transient model calibration, The Cryosphere, 9, 2429–2446, https://doi.org/10.5194/tc-9-2429-2015, 2015.

Greene, C. A., Gardner, A. S., and Andrews, L. C.: Detecting seasonal ice dynamics in satellite images, The Cryosphere, 14, 4365–4378, https://doi.org/10.5194/tc-14-4365-2020, 2020.

Guo, L., Li, J., Li, Z.-w., Wu, L.-x., Li, X., Hu, J., Li, H.-l., Li, H.-y., Miao, Z.-l., and Li, Z.-Q.: The Surge of the Hispar Glacier, Central Karakoram: SAR 3-D Flow Velocity Time Series and Thickness Changes, J. Geophys. Res.-Sol. Ea., 125, e2019JB018945, https://doi.org/10.1029/2019JB018945, 2020.

Gupta, H. V., Kling, H., Yilmaz, K. K., and Martinez, G. F.: Decomposition of the mean squared error and NSE performance criteria:

Implications for improving hydrological modelling, J. Hydrol., 377, 80–91, 2009.

Halas, P., Mouginot, J., de Fleurian, B., and Langebroek, P. M.: Southwest Greenland Ice Sheet Yearly Ice Velocities dataset from 1984 to 2020, Version 1, Zenodo [data set], https://doi.org/10.5281/zenodo.7418361, 2022.

Halas, P., Mouginot, J., de Fleurian, B., and Langebroek, P. M.: Impact of seasonal fluctuations of ice velocity on decadal trends observed in Southwest Greenland, Remote Sens. Environ., 285, 113419, https://doi.org/10.1016/j.rse.2022.113419, 2023.

Heid, T. and Kääb, A.: Evaluation of existing image matching methods for deriving glacier surface displacements globally from optical satellite imagery, Remote Sens. Environ., 118, 339–355, https://doi.org/10.1016/j.rse.2011.11.024, 2012.

Hetland, E. A., Musé, P., Simons, M., Lin, Y. N., Agram, P. S., and DiCaprio, C. J.: Multiscale InSAR Time Series (MInTS) analysis of surface deformation, J. Geophys. Res.-Sol. E., 117, https://doi.org/10.1029/2011JB008731, 2012.

Huber, P. J.: Robust estimation of a location parameter, in: Breakthroughs in statistics, 492–518, Springer, https://doi.org/10.1007/978-1-4612-4380-9_35, 1992.

Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss, M., Dussaillant, I., Brun, F., et al.: Accelerated global glacier mass loss in the early twenty-first century, Nature, 592, 726–731, 2021.

Jawak, S. D., Kumar, S., Luis, A. J., Bartanwala, M., Tummala, S., and Pandey, A. C.: Evaluation of geospatial tools for generating accurate glacier velocity maps from optical remote sensing data, Multidisciplinary Digital Publishing Institute Proceedings, 2, 341, https://doi.org/10.3390/e crs-2-05154, 2018.

Jay-Allemand, M., Gillet-Chaulet, F., Gagliardini, O., and Nodet, M.: Investigating changes in basal conditions of Variegated Glacier prior to and during its 1982–1983 surge, The Cryosphere, 5, 659–672, https://doi.org/10.5194/tc-5-659-2011, 2011.

Joughin, I., Smith, B. E., and Howat, I. M.: A complete map of Greenland ice velocity derived from satellite data collected over 20 years, J. Glaciol., 64, 1–11, 2018.

Kääb, A., Jacquemart, M., Gilbert, A., Leinss, S., Girod, L., Huggel, C., Falaschi, D., Ugalde, F., Petrakov, D., Chernomorets, S., Dokukin, M., Paul, F., Gascoin, S., Berthier, E., and Kargel, J. S.: Sudden large-volume detachments of low-angle mountain glaciers – more frequent than thought?, The Cryosphere, 15, 1751–1785, https://doi.org/10.5194/tc-15-1751-2021, 2021.

Kochtitzky, W., Copland, L., Van Wychen, W., Hugonnet, R., Hock, R., Dowdeswell, J. A., Benham, T., Strozzi, T., Glazovsky, A., Lavrentiev, I., Rounce, D. R., Millan, R., Cook, A., Dalton, A., Jiskoot, H., Cooley, J., Jania, J., and Navarro, F.: The unquantified mass loss of Northern Hemisphere marine-terminating glaciers from 2000–2020, Nature Communications, 13, https://doi.org/10.1038/s41467-022-33231-x, 2022.

Lacroix, P., Araujo, G., Hollingsworth, J., and Taipe, E.: Self-Entrainment Motion of a Slow-Moving Landslide Inferred From Landsat-8 Time Series, J. Geophys. Res.-Earth, 124, 1201–1216, https://doi.org/10.1029/2018JF004920, 2019.

Lauknes, T. R., Zebker, H. A., and Larsen, Y.: InSAR deformation time series using an $L_{(1)}$-norm small-baseline approach, IEEE transactions on geoscience and remote sensing, 49, 536–546, https://doi.org/10.1109/TGRS.2010.2051951, 2010.
Lei, Y., Gardner, A. S., and Agram, P.: Processing methodology for the ITS_LIVE Sentinel-1 ice velocity products, Earth Syst. Sci. Data, 14, 5111–5137, https://doi.org/10.5194/essd-14-5111-2022, 2022.

Liang, H., Zhang, L., Ding, X., Lu, Z., Li, X., Hu, J., and Wu, S.: Suppression of Coherence Matrix Bias for Phase Linking and Ambiguity Detection in MTInSAR, IEEE Transactions on Geoscience and Remote Sensing, 59, https://doi.org/10.1109/IGRS.2020.3000991, 2020.

López, L. A., Gardner, A. S., Greene, C. A., Kennedy, J. H., Liukis, M., Fahnestock, M. A., Scambos, T., and Fahnestock, J. R.: ITS_LIVE: A Cloud-Native Approach to Monitoring Glaciers From Space, Comput. Sci. Eng., 25, 49–56, 2023.

López-Quiroz, P., Doin, M.-P., Tupin, F., Briole, P., and Nicolas, J.-M.: Time series analysis of Mexico City subsidence constrained by radar interferometry, J. Appl. Geophys., 69, 1–15, 2009.

Maier, N., Humphrey, N., Harper, J., and Meierbachtol, T.: Sliding dominates slow-flowing margin regions, Greenland Ice Sheet, Sci. Adv., 5, https://doi.org/10.1126/sciadv.aaw5406, 2019.

Maier, N., Humphrey, N., Meierbachtol, T., and Harper, J.: Deformation motion tracks sliding changes through summer, western Greenland, J. Glaciol., 68, 187–196, 2022.

Maier, N., Andersen, J. K., Mouginot, J., Gimbert, F., and Gagliardi, N.: Wintertime supraglacial lake drainage cascade triggers large-scale ice flow response in Greenland, Geophys. Res. Lett., 50, https://doi.org/10.1029/2022GL02251, 2023.

Main, B., Copland, L., Smeda, B., Kochtitzky, W., Samsonov, S., Dudley, J., Skidmore, M., Dow, C., Van Wychen, W., Medrzycka, D., Higgs, E., and Mingo, L.: Terminus change of Kaskawulsh Glacier, Yukon, under a warming climate: retreat, thinning, slowdown and modified proglacial lake geometry, J. Glaciol., 69, 936–952, 2023.

Maronna, R. A., Martin, R. D., Yohai, V. J., and Salibián-Barrera, M.: Robust statistics: Theory and methods (with R) (2nd ed.), John Wiley & Sons, https://doi.org/10.1002/9781119214656, 2019.

Martin, M. D., Barr, I., Edwards, B., Spagnolo, M., Symeonakis, E., Mallalieu, J., Adamson, K., Mullan, D., and de Vries, M. V. W.: Glacier speed-up as a possible precursor to volcanic eruptions at Mount Veniaminof, Alaska, J. Glaciol., 1–18, https://doi.org/10.1017/jog.2024.107, 2025.

Millan, R., Mouginot, J., Rabatel, A., Jeong, S., Cusicanqui, D., Derkacheva, A., and Chekki, M.: Mapping Surface Flow Velocity of Glaciers at Regional Scale Using a Multiple Sensors Approach, Remote Sens., 11, 2498, https://doi.org/10.3390/rs11212498, 2019.

Millan, R., Mouginot, J., Rabatel, A., and Morlighem, M.: Ice velocity and thickness of the world's glaciers, Nat. Geosci., 15, 124–129, https://doi.org/10.1038/s41561-021-00885-z, 2022.

Morlighem, M., Rignot, E., Seroussi, H., Larour, E., Ben Dhia, H., and Aubry, D.: A mass conservation approach for mapping glacier ice thickness, Geophys. Res. Lett., 38, https://doi.org/10.1029/2011GL048659, 2011.

Mouginot, J., Scheuchl, B., and Rignot, E.: Mapping of ice motion in Antarctica using synthetic-aperture radar data, Remote Sens., 4, 2753–2767, 2012.

Mouginot, J., Rignot, E., Bjørk, A. A., Van den Broeke, M., Milan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.:

Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, P. Natl. Acad. Sci. USA, 116, 9239–9244, 2019.

Mouginot, J., Rabatel, A., Ducasse, E., and Millan, R.: Optimization of cross correlation algorithm for annual mapping of alpine glacier flow velocities; application to Sentinel-2, IEEE Trans. Geosci. Remote Sens., 61, 1–12, 2023.

Nanni, U., Scherler, D., Ayoub, F., Millan, R., Herman, F., and Avouac, J.-P.: Climatic control on seasonal variations in mountain glacier surface velocity, The Cryosphere, 17, 1567–1583, https://doi.org/10.5194/tc-17-1567-2023, 2023.

NASA Earth Observatory: Snow Swamp on Lowell Glacier, https://earthobservatory.nasa.gov/images/92699/snow-swamp-on-lowell-glacier (last access: 20 September 2025), 2018.

Pepe, A., Solaro, G., Calo, F., and Dema, C.: A minimum acceleration approach for the retrieval of multiplatform InSAR deformation time series, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 9, 3883–3898, 2016.

Provost, F., Zigone, D., Le Meur, E., Malet, J.-P., and Hbert, C.: Surface dynamics and history of the calving cycle of Astrolabe Glacier (Adélie Coast, Antarctica) derived from satellite imagery, The Cryosphere, 18, 3067–3079, https://doi.org/10.5194/tc-18-3067-2024, 2024.

Quincey, D. J., Glasser, N. F., Cook, S. J., and Luckman, A.: Heterogeneity in Karakoram glacier surges, J. Geophys. Res.-Earth, 120, 1288–1300, https://doi.org/10.1002/2015JF003515, 2015.

Rabatel, A., Ducasse, E., Millan, R., and Mouginot, J.: Satellite-derived annual glacier surface flow velocity products for the European alps, 2015–2021, Data, 8, 66, https://doi.org/10.3390/data8040066, 2023.

Rashid, I., Majeed, U., Jan, A., and Glasser, N. F.: The January 2018 to September 2019 surge of Shisper Glacier, Pakistan, detected from remote sensing observations, Geomorphology, 351, 106957, https://doi.org/10.1016/j geomorph.2019.106957, 2020.

RGI 7.0 Consortium: Randolph Glacier Inventory – A Dataset of Global Glacier Outlines, Version 7.0, https://doi.org/10.5067/f6jmovy5navz, [Data Set], 2023.

Riel, B., Minchew, B., and Joughin, I.: Observing traveling waves in glaciers with remote sensing: new flexible time series methods and application to Sermeq Kujalleq (Jakobshavn Isbræ), Greenland, The Cryosphere, 15, 407–429, https://doi.org/10.5194/tc-15-407-2021, 2021.

Round, V., Leinss, S., Huss, M., Haemmig, C., and Hajnsek, I.: Surge dynamics and lake outbursts of Kyagar Glacier, Karakoram, The Cryosphere, 11, 723–739, https://doi.org/10.5194/tc-11-723-2017, 2017.

Ruggiero, G. A., Cosme, E., Brankart, J.-M., Le Sommer, J., and Ubelmann, C.: An efficient way to account for observation error correlations in the assimilation of data from the future SWOT high-resolution altimeter mission, J. Atmos. Ocean. Technol., 33, 2755–2768, 2016.

Samsonov, S., Tiampo, K., and Cassotto, R.: Measuring the state and temporal evolution of glaciers in Alaska and Yukon using synthetic-aperture-radar-derived (SAR-derived) 3D time series of glacier surface flow, The Cryosphere, 15, 4221–4239, https://doi.org/10.5194/tc-15-4221-2021, 2021.

Samsonov, S. V. and d'Oreye, N.: Multidimensional Small Baseline Subset (MSBAS) for Two-Dimensional Deformation Analysis:
Case Study Mexico City, Can. J. Remote Sens., 43, 318–329, https://doi.org/10.1080/07038992.2017.1344926, 2017.

Sole, A., Nienow, P., Bartholomew, I., Mair, D., Cowton, T., Tedstone, A., and King, M. A.: Winter motion mediates dynamic response of the Greenland Ice Sheet to warmer summers: INTERANNUAL GREENLAND ICE FLOW CHANGES, Geophys. Res. Lett., 40, 3940–3944, https://doi.org/10.1002/grl.50764, 2013.

Tedstone, A. J., Nienow, P. W., Sole, A. J., Mair, D. W., Cowton, T. R., Bartholomew, I. D., and King, M. A.: Greenland ice sheet motion insensitive to exceptional meltwater forcing, P. Natl. Acad. Sci. USA, 110, 19719–19724, 2013.

Tedstone, A. J., Nienow, P. W., Gourmelen, N., Dehecq, A., Goldberg, D., and Hanna, E.: Decadal slowdown of a land-terminating sector of the Greenland Ice Sheet despite warming, Nature, 526, 692–695, 2015.

The GlaMBIE Team: Community estimate of global glacier mass changes from 2000 to 2023, Nature, 639, 382–388, https://doi.org/10.1038/s41586-024-08545-z, 2025.

van de Wal, R. S. W., Smeets, C. J. P. P., Boot, W., Stoffelen, M., van Kampen, R., Doyle, S. H., Wilhelms, F., van den Broeke, M. R., Reijmer, C. H., Oerlemans, J., and Hubbard, A.: Self-regulation of ice flow varies across the ablation area in south-west Greenland, The Cryosphere, 9, 603–611, https://doi.org/10.5194/tc-9-603-2015, 2015.

Van Wychen, W., Bayer, C., Copland, L., Brummell, E., and Dow, C.: Radarsat Constellation Mission Derived Winter Glacier Velocities for the St. Elias Icefield, Yukon/Alaska: 2022 and 2023, Can. J. Remote Sen., 49, 2264395, https://doi.org/10.1080/07038992.2023.2264395, 2023.

Van Wyk de Vries, M. and Wickert, A. D.: Glacier Image Velocimetry: an open-source toolbox for easy and rapid calculation of high-resolution glacier velocity fields, The Cryosphere, 15, 2115–2132, https://doi.org/10.5194/tc-15-2115-2021, 2021.

Vogel, C. R.: Non-convergence of the L-curve regularization parameter selection method, Inverse Problems, 12, 535, https://doi.org/10.1088/0266-5611/12/4/013, 1996.

Waechter, A., Copland, L., and Herdes, E.: Modern glacier velocities across the Icefield Ranges, St Elias Mountains, and variability at selected glaciers from 1959 to 2012, J. Glaciol., 61, 624–634, 2015.

Wallis, B. J., Hogg, A. E., van Wessem, J. M., Davison, B. J., and van den Broeke, M. R.: Widespread seasonal speed-up of west Antarctic Peninsula glaciers from 2014 to 2021, Nat. Geosci., 16, 231–237, 2023.

Wendleder, A., Bramboeck, J., Izzard, J., Erbertseder, T., d'Angelo, P., Schmitt, A., Quincey, D. J., Mayer, C., and Braun, M. H.: Velocity variations and hydrological drainage at Baltoro Glacier, Pakistan, The Cryosphere, 18, 1085–1103, https://doi.org/10.5194/tc-18-1085-2024, 2024.

Williams, J. J., Gourmelen, N., and Nienow, P.: Dynamic response of the Greenland ice sheet to recent cooling, Sci. Rep., 10, 1647, https://doi.org/10.1038/s41598-020-58355-2, 2020.

Williams, J. J., Gourmelen, N., and Nienow, P.: Complex multi-decadal ice dynamical change inland of marine-terminating glaciers on the Greenland Ice Sheet, J. Glaciol., 67, 833–846, 2021.

Yang, R., Hock, R., Kang, S., Guo, W., Shangguan, D., Jiang, Z., and Zhang, Q.: Glacier surface speed variations on the Kenai Peninsula, Alaska, 2014–2019, J. Geophys. Res.-Earth, 127, e2022JF006599, https://doi.org/10.1029/2022JF006599, 2022.

Yunjun, Z., Fattahi, H., and Amelung, F.: Small baseline InSAR time series analysis: Unwrapping error correction and noise reduction, Comput. Geosci., 133, 104331, https://doi.org/10.1016/j.cageo.2019.104331, 2019.

Zheng, W., Bhushan, S., Van Wyk De Vries, M., Kochtitzky, W., Shean, D., Copland, L., Dow, C., Jones-Ivey, R., and Pérez, F.: GLAcier Feature Tracking testkit (GLAFT): a statistically and physically based framework for evaluating glacier velocity products derived from optical satellite image feature tracking, The Cryosphere, 17, 4063–4078, https://doi.org/10.5194/tc-17-4063-2023, 2023.