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

HAL
open science

# Exploiting high-slip flow regimes to improve inference of glacier bed topography

Alexi Morin, Gwenn Flowers, Andrew Nolan, Douglas Brinkerhoff, Etienne Berthier

► To cite this version:

Alexi Morin, Gwenn Flowers, Andrew Nolan, Douglas Brinkerhoff, Etienne Berthier. Exploiting high-slip flow regimes to improve inference of glacier bed topography. Journal of Glaciology, 2023, 69 (275), pp.658-664. 10.1017/jog.2022.121. hal-04240456

HAL Id: hal-04240456

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

Submitted on 20 Oct 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

Letter

*Present address: Institut National de la Recherche Scientifique, Quebec City, Quebec, Canada.

**Cite this article:** Morin A, Flowers GE, Nolan A, Brinkerhoff D, Berthier E (2023). Exploiting high-slip flow regimes to improve inference of glacier bed topography. *Journal of Glaciology* 69(275), 658–664. https://doi.org/10.1017/jog.2022.121

Received: 16 August 2022
Revised: 11 November 2022
Accepted: 9 December 2022
First published online: 19 January 2023

**Keywords:**
Glacier modelling; glacier surges; ice dynamics; ice thickness measurements; mountain glaciers

**Author for correspondence:**
Gwenn E. Flowers, E-mail: gflowers@sfu.ca

© The Author(s), 2023. Published by Cambridge University Press on behalf of The International Glaciological Society. This is an Open Access article, distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike licence (http://creativecommons.org/licenses/by-nc-sa/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the same Creative Commons licence is used to distribute the re-used or adapted article and the original article is properly cited. The written permission of Cambridge University Press must be obtained prior to any commercial use.

# Exploiting high-slip flow regimes to improve inference of glacier bed topography

Alexi Morin¹,* ID, Gwenn E. Flowers¹ ID, Andrew Nolan¹ ID, Douglas Brinkerhoff² and Etienne Berthier³ ID

¹Department of Earth Sciences, Simon Fraser University, Burnaby, British Columbia, Canada; ²Department of Computer Science, University of Montana, Missoula, Montana, USA and ³LEGOS, Université de Toulouse, CNES, CNRS, IRD, UPS, 31400 Toulouse, France

## Abstract

Theory and observation show that glacier-flow regimes characterized by high basal slip enhance the projection of topographic detail to the surface, motivating this investigation into the efficacy of using glacier surges to improve bed estimation. Here we adapt a Bayesian inversion scheme and apply it to real and synthetic data as a proof of concept. Synthetic tests show a reduction in mean RMSE between true and inferred beds by more than half, and an increase in the mean correlation coefficient of ~0.5, when data from slip- versus deformation-dominated regimes are used. Multi-epoch inversions, which partition slip- and deformation-dominated regimes, are shown to outperform inversions that average over these flow regimes thereby squandering information. Tests with real data from a surging glacier in Yukon, Canada, corroborate these results, while highlighting the challenges of limited or inconsistent data. With the growing torrent of satellite-based observations, fast-flow events such as glacier surges offer potential to improve bed estimation for some of the world's most dynamic glaciers.

## 1. Introduction and motivation

Knowledge of ice thickness is pre-requisite to any study of ice dynamics and fundamental to solving glaciological problems related to meltwater supply, sea-level rise and hazard potential (e.g. Gilbert and others, 2018; Azam and others, 2021; Fox-Kemper and others, 2021; Millan and others, 2022). This need has fuelled the demand for estimates of glacier thickness, resulting in the development of numerous models that infer ice thickness or bed topography from glacier-surface observables such as elevation, elevation-change rate, mass balance and velocity (see Farinotti and others, 2017, 2020, for an overview). Distributed estimates of glacier thickness have now been made on a global scale (e.g. Farinotti and others, 2019; Millan and others, 2022), and are seeing increased use in regional and local studies despite largely unquantified site-specific uncertainties (cf. Pelto and others, 2020; Werder and others, 2020).

A vigorous surge of an unnamed glacier (RGI ID: RGI50-01.16198, GLIMS ID: G220578E60873N, Fig. 1) from 2017 to 2018 in Yukon, Canada, prompted interest in its study (Main and others, 2019; Samsonov and others, 2021a), and a corresponding need for estimates of ice thickness/bed topography. Theory (e.g. Gudmundsson, 2003) and observation (e.g. Gudmundsson and others, 2003; De Rydt and others, 2013; Sharp, 2021) demonstrate that high-slip flow regimes result in surface projection of shorter-wavelength bed undulations than are detectable in deformation-dominated regimes. Here we aim to test the hypothesis that inference of glacier bed topography can be improved by exploiting high-slip events such as glacier surges. To do this, we present a simple extension of the inference scheme of Brinkerhoff and others (2016) that permits its application to multiple observational epochs. This Bayesian method casts the problem of inferring the glacier bed as a statistical one, with the end result a probability distribution constrained both by data and a priori assumptions such as smoothness and topographic amplitude. By modelling each constituent quantity in the glaciological continuity equation as a random variable, the resulting distribution also provides robust estimates of uncertainty, a critical advantage over alternative deterministic methods. We apply this scheme to synthetic data as a proof of concept, and then to real data from the site in Figure 1.

## 2. Bayesian inference

We seek to characterize the posterior distribution of model variables

$$P(\mathbf{m}|\mathbf{d}) \propto P(\mathbf{d}|\mathbf{m}) P(\mathbf{m}), \quad (1)$$

where vector **m** contains model variables (e.g. surface velocity $U_s$, mass balance $\dot{b}$, ice thickness **H**), which we hereafter refer to as parameters, and vector **d** contains the data (i.e. any observations of the model variables). Bayes' theorem allows us to recast the posterior distribution $P(\mathbf{m}|\mathbf{d})$, which represents the joint probability distribution over the parameters constrained by the data, as proportional to (α) the product of the likelihood $P(\mathbf{d}|\mathbf{m})$, which measures the probability of the data under a hypothesized parameter value, and prior $P(\mathbf{m})$, which

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 1. Unnamed former tributary to Kluane Glacier in the Traditional Territory of the Kluane and White River First Nations. Pre- (yellow) and post- (purple) surge out- lines are shown, along with ice-thickness measurement locations (orange). Flowline (black with km markers) constructed following Kienholz and others (2014) with tools available in the Open Global Glacier Model (Maussion and others, 2019) and manually adjusted to improve agreement with velocity field. Background image from Sentinel-2, 29 August 2018 (UTM Zone 7N).</figcaption></figure>

encodes initial assumptions about the parameters prior to consid- eration of observations (Tarantola, 2005). We model the likeli- hood by assuming that spatially and temporally distributed observations are independent of one another and normally dis- tributed around a mean value that is given directly by the para- meters (in the case of the bed elevation and the rate of surface elevation change) or can be computed via the solution of the gla- ciological continuity equation (as in the case of $U_s$, see below).

We place Gaussian process priors over the bed elevation and rate of surface elevation change which assume that the value of these parameters at any finite collection of spatial locations is multivariate normal. Such priors require assumptions about spatial variability (Table S1). We adopt a Gaussian covariance function

$$K(r, r') = \sigma^2 \exp\left(-\left(\frac{r-r'}{\ell}\right)^2\right), \quad (2)$$

with hyper-parameters $\sigma$ (amplitude) and $\ell$ (correlation length scale) and $r-r'$ the distance between any two points on the flowline (Fig. S1). The same priors are used for real and syn- thetic data (Table S1), and can be interpreted as uncertainty on input data or as a regularization parameter. Priors on individual variables are described in the Supplementary material.

## 2.1. Continuity equation

For constant bed elevation $B(\mathbf{x}) = S(\mathbf{x}) - H(\mathbf{x})$, continuity can be written

$$\frac{\partial H}{\partial t}(\mathbf{x}) = \frac{\partial S}{\partial t}(\mathbf{x}) = -\nabla \cdot [\bar{U}(\mathbf{x}) H(\mathbf{x}) N(\mathbf{x})] + \dot{b}(\mathbf{x}), \quad (3)$$

for ice thickness $H$, time $t$, map-view coordinate $\mathbf{x}$, surface eleva- tion $S$, depth-averaged flow speed $\bar{U}$, flow direction $N$ and surface mass balance $\dot{b}$. Following Brinkerhoff and others (2016) we reduce the 2-D formulation to one, from the map-view coordinate $\mathbf{x}$ to the flowline coordinate $r$, with flow direction $N(\mathbf{x})$ replaced by the flow-unit width between two streamlines $w(r)$:

$$w(r) \frac{\partial S}{\partial t}(r) = -\frac{\partial}{\partial r} [w(r) \bar{U}(r) H(r)] + w(r) \dot{b}(r), \quad (4)$$

and time-average the forward model over the observation period $\Delta t = t_1 - t_0$, such that

$$w(r) \frac{\Delta S}{\Delta t} = -\frac{1}{\Delta t} \int_{t_0}^{t_1} \left( \frac{\partial}{\partial r} [w(r) \bar{U}(r) H(r)] - w(r) \dot{b}(r) \right) dt, \quad (5)$$

where $\Delta S$ is the change in surface elevation. Writing each term as an average over the observation period, replacing the depth-averaged flow speed $\bar{U}(r)$ with the observable surface flow speed $U_s(r) = s \bar{U}(r)$ and solving Eqn (5) for flow speed yields

$$U_s(r) = \frac{s}{w(r) H(r)} \int_0^r w(r') \left( \dot{b}(r') - \frac{\Delta S(r')}{\Delta t} \right) dr', \quad (6)$$

where the flow-unit width $w(r)$ as well as parameter $s$ are deter- mined by the inversion. We adapt the original model above first by replacing ice thickness $H(r,t)$ with $S(r,t) - B(r)$, and then writing Eqn (6) to accommodate multiple observational epochs:

$$U_s(r, \Delta t_m) = \frac{s(\Delta t_m)}{w(r)(S(r, \Delta t_m) - B(r))} \int_0^r w(r') \left( \dot{b}(r', \Delta t_m) - \frac{\Delta S(r', \Delta t_m)}{\Delta t_m} \right) dr', \quad (7)$$

where each term represents an average over an observational epoch $\Delta t_m$. Epochs need not be contiguous in time. These small changes to the model allow us to isolate data from time intervals associated with different flow regimes, and use multiple observational epochs to reconstruct a single bed. Hereafter, we use the term ‘full epoch’ to describe inversions where averages over the full observation period are used, as in the original model (Eqn (6)).

## 2.2. Model implementation and evaluation

The product of the likelihood and prior distributions is only proportional to the posterior distribution, so we approximate the posterior as is standard in Bayesian inference. Here we use the Metropolis-Hastings algorithm (Robert and Casella, 1999), a Markov chain Monte Carlo (MCMC) method, to draw a collec- tion of samples from the posterior distribution (Eqn (1)), which
can then be used to characterize the distribution in the sense of a histogram.

At least three Monte Carlo Markov chains are run for a minimum of $10^6$ iterations each, with the first $10^5$ results discarded ('burn-in'), and only 1 in 10 results retained ('thinning') to avoid auto-correlation (Gelman and others, 1995). Convergence of the Markov chains is assessed using the Gelman–Rubin statistic. We quantify model performance by computing the RMSE (Chai and Draxler, 2014) between the true and posterior bed:

$$ \text{RMSE} = \sqrt{\frac{\sum_{i=1}^{n} (B_i - \hat{B}_i)^2}{n}}, \quad (8) $$

where $B_i$ is the true bed elevation at discrete location $i$ of $n$, and $\hat{B}_i$ is the co-located mean of the inferred bed posterior. We also compute the Pearson correlation coefficient $r = [-1, 1]$ (Wasserman, 2010) defined as

$$ r(\delta B, \delta \hat{B}) = \frac{\text{Cov}(\delta B, \delta \hat{B})}{s_{\delta B s_{\delta \hat{B}}}}, \quad (9) $$

for the true ($\delta B$) and inferred ($\delta \hat{B}$) bed perturbations, and their respective sampled standard deviations $s_{\delta B}$ and $s_{\delta \hat{B}}$. Bed perturbations are defined by subtracting the same reference bed from any given bed profile. The purpose of using perturbations is to isolate and compare short-wavelength variations in bed topography ($\le 10H$) without the influence of long-wavelength bed structure, which is already assessed by the RMSE.

We use a horizontal model grid spacing of 200 m. Hyper-parameter values (Table S1) are chosen to be realistic and representative of the study glacier, and real observational uncertainties. Some trial and error is nevertheless required to determine values that yield plausible results.

In what follows, we perform four inversions on each dataset: one each for a deformation-dominated ('quiescent') glacier flow regime, a slip-dominated ('surge') flow regime, a single time interval that includes both quiescent and surge regimes ('full-epoch inversion') and multiple time intervals that include but partition quiescent and surge regimes ('multi-epoch inversion'). The first three of these use the original model in Eqn (6), while the multi-epoch inversion uses Eqn (7).

# 3. Synthetic case

## 3.1. Data

We use the open-source finite-element model Elmer/Ice (Gagliardini and others, 2013) to generate synthetic glacier profiles, and associated synthetic data, for different glacier beds. We solve the Stokes equations in two dimensions ($x-z$ flowband) for incompressible flow using a Glen–Nye-type constitutive law for temperate ice. Effects of variable flowband width and lateral drag are neglected, and a linear friction law is used as the basal boundary condition (Gagliardini and others, 2013). The numerical mesh has a horizontal spacing of 50 m and 10 vertical layers.

We use smoothed surface and bed profiles from Farinotti and others (2019), hereafter 'F2019', extracted along the centreline of the surging tributary in Figure 1 to define the initial model geometry. A series of synthetic glacier beds is produced by adding sinusoidal perturbations to the F2019 bed $B_{\text{F2019}}$: $B_{\text{synth}_k} = B_{\text{F2019}} + A_k \sin(2\pi/k\hat{H})$, where $A_k = \lambda_k R$ (m) and $\lambda_k$ (m) are the amplitude and wavelength of the $k$th bed perturbation, respectively, $R = 0.01$ is the prescribed amplitude-to-wavelength ratio, $\hat{H} = 100$ m is the characteristic ice thickness and $k = [5, 6, ..., 15]$. We choose the value of $R$ such that perturbations

<figure><img src="image_6.png" /><figcaption>**Fig. 2.** Synthetic input data (solid lines) for composite-bed inversions: deformation-only (quiescent) regime in blue, high-slip (surge) regime in orange and full epoch in purple (time-weighted average of blue and orange curves). Multi-epoch inversion uses data in blue and orange. Shading indicates prescribed std dev. used to represent observational uncertainty. (a) Surface flow speed. (b) Surface-elevation change rate.</figcaption></figure>

are significant in amplitude but do not change the overall shape of the bed (see Supplementary material). The value of $\hat{H}$ is based on the mean ice thickness of the reference geometry. The minimum value of $k$ corresponds to a perturbation wavelength of $k\hat{H} = 500$ m, and is chosen to be slightly larger than the Nyquist wavelength of 400 m determined from the 200 m horizontal grid spacing of the inverse model. Finally, a composite bed $B_{\text{synth}}$ is also defined by adding the sum of all perturbations for $k = [5, 6, ..., 15]$ to $B_{\text{F2019}}$:

$$ B_{\text{synth}} = B_{\text{F2019}} + \sum_{k=5}^{15} A_k \sin \left( \frac{2\pi}{k\hat{H}} \right). \quad (10) $$

To produce self-consistent synthetic datasets that represent both deformation- and slip-dominated glacier flow regimes, we first compute steady-state surface profiles (with no sliding) for each bed that resemble the real glacier in its pre-surge state. We do this experimentally by applying offsets to prescribed annual surface mass balances, which have been modelled for the region (Young and others, 2021), and solving the continuity equation. Positive offsets of $1.95-2.1$ m a$^{-1}$ ice-equivalent were required to produce synthetic glaciers of comparable volume to the real glacier prior to the surge (Supplementary material). We use the resulting glacier geometry and velocities to represent the deformation-dominated (quiescent) regime, to which we associate an arbitrarily prescribed time interval of 10 years, equal to the time interval associated with the slip-dominated flow regime. Although 10 years is shorter than most observed quiescent intervals (e.g. Meier and Post, 1969), it is a plausible interval over which the necessary data for these inversions may be available. A more realistic representation of quiescence would include a non-steady surface profile in which mass accumulates above and diminishes below the dynamic balance line.

To simulate a slip-dominated flow regime, we begin with the steady-state glacier and instantaneously reduce the slip coefficient $\beta$ from 1 to $3.5 \times 10^{-4}$ MPa a m$^{-1}$ over the entire domain for another 10 years. This simulation produces terminus advance comparable to that observed during the 2016–18 surge of the study glacier (Fig. 1). Requirements for numerical stability prevent us from simulating the most extreme flow speeds associated with surges, however, the high-slip regimes we simulate are characterized by flow speeds 5–10 times higher than those associated with quiescence and basal contributions to surface flow speed in excess
of 80% (Fig. 2 and Supplementary material). These flow speeds are sufficient to achieve our aim of producing plausible fields of model variables from a slip-dominated flow regime.

The synthetic data include surface velocities and elevation change rates for quiescent and surge-like flow regimes on each of 12 different glacier beds (see Fig. 2 for the case of the composite bed). For use in the inversions, the synthetic data are time-averaged through trapezoidal integration. We also specify the locations of zero ice thickness at the head and terminus of the glacier, as well as known bed elevations where the bed is exposed prior to glacier advance.

## 3.2. Results

Inversions of synthetic data demonstrate the increase in performance made possible by isolating the high-slip (surge) flow regime (Figs 3, 4). Using data from the high-slip regime alone yields considerably more bed structure and a narrower posterior distribution than using data from the deformation-only (quiescent) regime alone (Fig. 3b cf. Fig. 3a). It also reduces the mean RMSE between true and inferred bed by >60% (Fig. 4a) and increases the mean correlation coefficient between true and inferred bed perturbations by ~0.5 (Fig. 4b). Results of the multi-epoch inversion (Fig. 3d) closely resemble those of the high-slip case, demonstrating that this approach effectively captures the information content of the high-slip regime. In contrast,

<figure><img src="image_7.png" /><figcaption>**Fig. 3.** Posterior distributions of the bed for synthetic glacier. (a) Deformation-only (quiescent) regime. (b) High-slip (surge) regime. (c) Full-epoch inversion. (d) Multi-epoch inversion. Wider confidence intervals on the posteriors in (a) and (c) reflect less information content in the datasets. Black dots represent known bed elevations that are input to the model.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Error metrics for inversions using synthetic data and composite bed (Fig. 3). Histograms generated by computing respective error metrics between each true value of bed elevation or bed perturbation amplitude and all co-located realizations of the corresponding quantities in the posterior distributions. (a) RMSE (Eqn (8)) between posterior distributions of bed elevation and true bed elevation. (b) *r* (Eqn (9)) between posterior distributions of bed perturbations and true bed perturbations.</figcaption></figure>

information is lost in the full-epoch inversion (Fig. 3c), where variables are time-averaged over the full interval of observation which includes quiescence and surging. The performance of the full-epoch inversion is accordingly intermediate between the quiescent and surge-like regimes (Fig. 4). For more realistic prescriptions of the quiescent time interval (decades rather than 10 years) compared to the 10-year surge interval, the performance of the full-epoch inversion declines, making the multi-epoch approach even more advantageous (Supplementary material).

# 4. Real case

## 4.1. Data

Surface-elevation change rates were calculated from DEMs for 13 September 2007, 30 September 2016 and 1 October 2018. The 2007 and 2018 DEMs are derived from SPOT5 and SPOT7 data and were used in a previous study (Young and others, 2021), while the 2016 DEM is derived from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) stereo-images using the Ames Stereo Pipeline (Shean and others, 2016; Brun and others, 2017). The processing of the DEMs follows the workflow presented in Berthier and Brun (2019): briefly, all DEMs are coregistered to the global TanDEM-X DEM (Rizzoli and others, 2017) on stable terrain, masking out glaciers using the Randolph Glacier Inventory (RGI) v6.0 (Pfeffer and others, 2014). We associate elevation change rates from 2007 to 2016 with quiescence, and those from 2016 to 2018 with the surge (Fig. 5).

Surface velocities for the 2007–16 (quiescent) epoch are taken from the publicly available NASA MEaSURES Inter-Mission Time Series of Land Ice Velocity and Elevation (ITS_LIVE) image pairs (Gardner and others, 2019). ITS_LIVE velocity products are produced with auto-RIFT offset tracking (Gardner and others, 2018)
<figure><img src="image_9.png" /><figcaption>**Fig. 5.** Real input data (solid lines) for study glacier: 2007–16 quiescent regime in blue, 2016–18 surge regime in orange and 2007–18 full-epoch in purple (time- weighted average of blue and orange curves). Multi-epoch inversion uses data in blue and orange. Shading indicates 1 std dev. used to represent observational uncer- tainty. (a) Surface flow speed. (b) Surface-elevation change rate. Profiles taken along flowline in Figure 1.</figcaption></figure>

using Landsat 4, 5, 7 and 8 imagery. Velocities for the 2016–18 (surge) epoch are taken from Samsonov and others (2021a) and derived using synthetic aperture radar methods described in Samsonov and others (2021b). A combination of these two datasets are needed, because ITS_LIVE velocities are unreliable during the surge, while the dataset from Samsonov and others (2021a) begins only in 2016. A self-consistent velocity dataset spanning both quiescence and surging is currently under development (personal communication from L. Copland, 2022).

The mean surface mass-balance profile for 2007–18 is identical to that used to initialize the synthetic simulations (Young and others, 2021), but without the applied offset (see Supplementary material). Ground-based ice-thickness measurements were made in several locations after the surge in July 2019 and July 2021 (Fig. 1) using an impulse radar system (Mingo and Flowers, 2010) and resistively loaded dipole antennas of 10 MHz centre frequency. Bed elevations were derived by subtracting the ice-thickness measurements from the October 2018 (post-surge) surface DEM.

## 4.2. Results

Inversions of the real data bear a qualitative similarity to the synthetic results: inversion of data from the quiescent regime alone yields the smoothest bed and widest posterior distribution (Fig. 6a), while inversion of data from the surge regime alone and the multi-epoch inversion again bear a strong resemblance to one another (Figs 6b, 6d). Both produce shorter wavelength undulations and have narrower posteriors than the full-epoch inversion (Fig. 6c). Unlike the synthetic tests however, the bed shape arising from the quiescent data alone (Fig. 6a) is not simply a smoothed version of that arising from the surge regime. Instead, it has a unique structure, with a deeper bed along the first ~6 km of the flowline and shallower bed from ~6–13 km.

While the true bed of the entire study glacier is unknown, we compare the inversion results to sparse ice-thickness measurements (locations in Fig. 1, comparison in Fig. 6). Visual inspection suggests the inversion of data from the surge regime alone performs best (Fig. 6b) and quiescent regime alone worst (Fig. 6a), with the full- and multi-epoch inversions in between. Quantification of the mismatch between measured bed elevations and the posterior bed distributions in Figure 6 clearly illustrates this hierarchy of performance (Fig. 7), albeit on the basis of few observations.

<figure><img src="image_10.png" /><figcaption>**Fig. 6.** Posterior distributions of the bed for inversions using real data. Measurements of bed elevation are shown as filled circles where they intersect the flowline (Fig. 1) and as open circles where they are 95–125 m from the flowline. (a) Quiescent regime, 2007–16. (b) Surge regime, 2016–18. (c) Full-epoch inversion, 2007–18. (d) Multi-epoch inversion, 2007–16 and 2016–18. Black dots represent known bed elevations that are input to the model.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 7.** RMSE between bed posterior distributions for inversions using real data (Fig. 6) and measurements of bed elevation that intersect the flowline (filled circles in Fig. 6). Histograms were generated by computing RMSE between each measurement and all co-located realizations of bed elevation in the posterior distributions.</figcaption></figure>

# 5. Discussion and recommendations

## 5.1. Study shortcomings and limitations

This pilot study has demonstrated potential for data from high-slip glacier flow regimes to reveal bed topographic detail normally
obscured by deformation-dominated flow regimes. It is limited, however, in its application to a flowband rather than the full glacierr, and by the quality, consistency and completeness of the real data. The profile in Figure 1 crosses flow units in attempting to capture tributary and trunk glaciers, while data have been extracted along the profile rather than from a swath surrounding the profile and therefore do not represent laterally averaged values as a flowband model would assume (e.g. Sergienko, 2012). Comparison of inversion results with measurements and with other bed estimates would be more defensible if swath, rather than profile, data were used.

Bed topography reveals itself in high-slip regimes through perturbations to surface velocity and elevation (e.g. Figs 2, 5). The extent to which a 2-D flowband model represents 3-D reality therefore depends on how well the 2-D model represents slip-induced perturbations to surface velocity and elevation. 2-D flowband models have been shown to exaggerate perturbations in the presence of bedrock bumps, and in the case of surface velocity, even produce perturbation structures that differ from those of 3-D models (Sergienko, 2012). The extent to which the slip-induced amplification of these structures differs between 2-D and 3-D models is relevant here: Sergienko (2012) shows that absolute amplification of surface elevation and velocity perturbations with increased slip is greater for 2-D than 3-D models, but relative amplification is sometimes comparable. One exception is the case from no slip to slip, where elevation perturbations decrease slightly in 3-D, but increase in 2-D (Sergienko, 2012). While our use of 2-D flowband models leads to overestimation of the information content in high-slip regimes, and thus an optimistic assessment of the power of exploiting these regimes in inferences of bed topography, observations indeed demonstrate measurable amplification of 3-D surface perturbations following glacier surges (e.g. Sharp, 2021) suggesting our approach has merit.

The tidy and intuitive results of the synthetic tests benefit from internally consistent noise-free data (see Supplementary material for simulations with added noise) with shorter-wavelength spectral content compared to the corresponding real data (spectra not shown). Despite the lower magnitudes of the synthetic surge velocities and elevation changes rates, the information content of the synthetic datasets may still be higher. The synthetic data also benefit from an unambiguous transition between quiescent and surge-like flow regimes that is abrupt in time and occurs uniformly along the flowline. The extent to which this situation mimics reality varies (e.g. Frappé and Clarke, 2007).

The real data further suffer from some internal inconsistency, particularly between surge velocities and elevation change rates near the terminus, making satisfaction of the continuity equation impossible (Supplementary material). Such situations are not unique to this study, but can lead, for example, to unphysical posteriors such as negative surface velocities. We would anticipate improvement with internally consistent datasets, as well as with more precise temporal isolation of the high-slip regime which, in our case, was limited by the number of surface DEMs. The poor results associated with the real quiescent data tarnish even the multi-epoch inversion results. These data do not appear to suffer from the same inconsistency as the surge data (based on posterior distributions of flow speed, elevation change rate and mass balance, see Supplementary material), while sensitivity tests (not shown) point away from the mis-specification of hyperparameters or model priors as the problem. We hypothesize the poor performance arises from low signal-to-noise ratios associated with the exceptionally low quiescent surface velocities (Fig. 5a). Some improvement would be expected if bed-elevation measurements are used in defining priors, as opposed to being reserved for model evaluation as done here.

Finally, the mathematical methods used in this study suffer from numerical constraints when dealing with the high-resolution topography that is desirable as model input. The grid spacing we employed results in spatially correlated samples in the MCMC, leading to non-convergence of the traces for certain nodes in the mesh. While fewer than 10% of all nodes failed to converge in most inversions, the multi-epoch inversion with real data was an exception, with >50% nodes having a Gelman–Rubin statistic >1.1. Simulations with more than three chains and 5 × 10⁶ iterations yielded similar results and convergence statistics.

## 5.2. Outlook and future work

Logical extensions to this pilot study would include application of the method to 2-D plan-view bed inversions, and to surge-type glaciers with consistent and comprehensive data including widespread measurements of bed topography. One could extend the multi-epoch approach from two to any number of observational epochs. Compared to the full-epoch approach, the advantage of the multi-epoch approach is in isolating, rather than averaging over, high-slip regimes, which reveal additional information about bed topography (e.g. Gudmundsson, 2003). In the limit of numerous short epochs, the multi-epoch approach becomes one of solving the time-dependent problem. While using data from the surge-regime only may be the best approach, the multi-epoch approach represents an operationally convenient means of improving bed estimates for a spectrum of dynamic glaciers and has potential for automation. Numerous studies attest to a growing ability to detect and monitor glacier surges from space (e.g. Altena and others, 2019; Fu and others, 2019; Rashid and others, 2020; Hugonnet and others, 2021; Leclercq and others, 2021; Beaud and others, 2022). With this trend towards near-real-time processing of satellite data (e.g. Gardner and others, 2019) from which most of the required variables can be derived, we see an opportunity for ongoing improvement of bed estimation for glaciers that experience high-slip events such as surges.

**Supplementary material.** The supplementary material for this article can be found at https://doi.org/10.1017/jog.2022.121

**Data.** Inverse model code and inputs: https://github.com/aleximorin/TimeDeltaBayesianModel. Synthetic data and metadata: https://github.com/andrewdnolan/Harmonic_Beds. NASA ITS_LIVE surface velocities obtained from https://nsidc.org/apps/itslive/. Surface velocities during surge obtained from at https://doi.org/10.17632/zf67rsgydv.1 (Samsonov and others, 2021a). Surface DEMs: https://zenodo.org/record/6458189.

**Acknowledgements.** This research was enabled in part by West Grid (www.westgrid.ca) and Compute Canada (www.computecanada.ca). Financial support was provided by the Natural Sciences and Engineering Research Council of Canada, Polar Continental Shelf Program and Simon Fraser University. Field data were collected with permission and support from Kluane and White River First Nations, Parks Canada and Yukon Government. Mass-balance model output was provided by Erik Young. SPOT7 data were obtained thanks to public funds received in the framework of GEOSUD, a project (ANR-10-EQPX-20) of the programme 'Investissements d'Avenir' managed by the French National Research Agency. EB acknowledges support from the French Space Agency (CNES) through the TOSCA and DINAMIS projects. We are grateful for the publicly accessible data used in this study and input from S. Samsonov in particular. Thanks to our two reviewers for helping us write a better paper.

**Author contributions.** AM and GF conceived of the project. AM wrote the inversion code and generated the inversion results. AN produced the synthetic data. GF processed and interpreted the radar data. EB produced the surface DEMs. GF, AM and AN wrote the manuscript. DB provided guidance in implementing the inversion scheme and contributed to revisions. All authors reviewed the manuscript.
## References

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

Azam MF and 9 others (2021) Glaciohydrology of the Himalaya-Karakoram. *Science* 373(6557), 3668. doi: 10.1126/science.abf3668

Beaud F, Aati S, Delaney I, Adhikari S and Avouac JP (2022) 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(8), 3123–3148. doi: 10.5194/tc-16-3123-2022

Berthier F and Brun F (2019) Karakoram geodetic glacier mass balances between 2008 and 2016: persistence of the anomaly and influence of a large rock avalanche on Siachen Glacier. *Journal of Glaciology* 65(251), 494–507. doi: 10.1017/jog.2019.32

Brinkerhoff DJ, Aschwanden A and Truffer M (2016) Bayesian inference of subglacial topography using mass conservation. *Frontiers in Earth Science* 4, 1–15. doi: 10.3389/feart.2016.00008

Brun F, Berthier E, Wagnon P, Kääb A and Treichler D (2017) A spatially resolved estimate of High Mountain Asia glacier mass balances from 2000 to 2016. *Nature Geoscience* 10(9), 668–673. doi: 10.1038/ngeo2999

Chai T and Draxler RR (2014) Root mean square error (RMSE) or mean absolute error (MAE)? – arguments against avoiding RMSE in the literature. *Geoscientific Model Development* 7(3), 1247–1250. doi: 10.5194/gmd-7-1247-2014

De Rydt J, Gudmundsson GH, Corr H and Christoffersen P (2013) Surface undulations of Antarctic ice streams tightly controlled by bedrock topography. *The Cryosphere* 7(2), 407–417. doi: 10.5194/tc-7-407-2013

Farinotti D and 9 others (2017) How accurate are estimates of glacier ice thickness? Results from ITMIX, the Ice Thickness Models Intercomparison eXperiment. *The Cryosphere* 11(2), 949–970. doi: 10.5194/tc-11-949-2017

Farinotti D and 6 others (2019) A consensus estimate for the ice thickness distribution of all glaciers on Earth. *Nature Geoscience* 12(3), 168–173. doi: 10.1038/s41561-019-0300-3

Farinotti D and 9 others (2020) Results from the Ice Thickness Models Intercomparison eXperiment Phase 2 (ITMIX2). *Frontiers in Earth Science* 8, 484. doi: 10.3389/feart.2020.571923

Fox-Kemper B and 17 others (2021) Ocean, Cryosphere and Sea Level Change. In: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change.

Frappé TP and Clarke GKC (2007) Slow surge of Trapridge Glacier, Yukon Territory, Canada. *Journal of Geophysical Research* 112, F03S32. doi: 10.1029/2006JF000607

Fu X, Li Z and Zhou J (2019) Characterizing the surge behavior of Alakesayi Glacier in the West Kunlun Shan, Northwestern Tibetan Plateau, from remote-sensing data between 2013 and 2018. *Journal of Glaciology* 65 (249), 168–172. doi: 10.1017/jog.2019.2

Gagliardini O and 14 others (2013) Capabilities and performance of Elmer/Ice, a new-generation ice sheet model. *Geoscientific Model Development* 6(4), 1299–1318. doi: 10.5194/gmd-6-1299-2013

Gardner AS and 6 others (2018) Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years. *The Cryosphere* 12(2), 521–547. doi: 10.5194/tc-12-521-2018

Gardner AS, Fahnestock MA and Scambos TA (2019) ITS_LIVE Regional Glacier and Ice Sheet Surface Velocities. Data archived at National Snow and Ice Data Center, data archived at National Snow and Ice Data Center.

Gelman A, Carlin JB, Stern HS and Rubin DB (1995) Bayesian Data Analysis. New York: Chapman and Hall/CRC.

Gilbert A and 8 others (2018) Mechanisms leading to the 2016 giant twin glacier collapses, Aru Range, Tibet. *The Cryosphere* 12(9), 2883–2900. doi: 10.5194/tc-12-2883-2018

Gudmundsson GH (2003) Transmission of basal variability to a glacier surface. *Journal of Geophysical Research: Solid Earth* 108(B5), 2253. doi: 10.1029/2002JB002107

Gudmundsson GH, Adalgeirsdóttir G and Björnsson H (2003) Observational verification of predicted increase in bedrock-to-surface

amplitude transfer during a glacier surge. *Annals of Glaciology* 36, 91–96. doi: 10.3189/172756403781816248

Hugonnet R and 9 others (2021) Accelerated global glacier mass loss in the early twenty-first century. *Nature* 592(7856), 726–731. doi: 10.1038/s41586-021-03436-z

Kienholz C, Rich J, Arendt A and Hock R (2014) A new method for deriving glacier centerlines applied to glaciers in Alaska and northwest Canada. *The Cryosphere* 8(2), 503–519. doi: 10.5194/tc-8-503-2014

Leclercq PW, Kääb A and Altena B (2021) Brief communication: Detection of glacier surge activity using cloud computing of Sentinel-1 radar data. *The Cryosphere* 15(10), 4901–4907. doi: 10.5194/tc-15-4901-2021

Main B and 6 others (2019) Surge of Little Kluane Glacier in the St. Elias Mountains, Yukon, Canada, from 2017–2018. *AGU Fall Meeting Abstracts*, 31.

Maussion F and 9 others (2019) The open global glacier model (OGGM) v1.1. Geoscientific Model Development 12(3), 909–931. doi: 10.5194/gmd-12-909-2019

Meier MF and Post A (1969) What are glacier surges? *Canadian Journal of Earth Sciences* 6(4), 807–817. doi: 10.1139/e69-081

Millan R, Mouginot J, Rabatel A and Morlighem M (2022) Ice velocity and thickness of the world's glaciers. *Nature Geoscience* 15(2), 124–129. doi: 10.1038/s41561-021-00885-z

Mingo L and Flowers GE (2010) An integrated lightweight ice-penetrating radar system. *Journal of Glaciology* 56(198), 709–714. doi: 10.3189/002214310793146179

Pelto BM, Maussion F, Menounos B, Radić V and Zeuner M (2020) Bias-corrected estimates of glacier thickness in the Columbia River Basin, Canada. *Journal of Glaciology* 66(260), 1051–1063. doi: 10.1017/jog.2020.75

Pfeffer WT and 19 others (2014) The Randolph Glacier Inventory: a globally complete inventory of glaciers. *Journal of Glaciology* 60(221), 537–552. doi: 10.3189/2014JoG13J176

Rashid I, Majeed U, Jan A and Glasser NF (2020) The January 2018 to September 2019 surge of Shisper Glacier, Pakistan, detected from remote sensing observations. *Geomorphology* 351 106957. doi: 10.1016/j.geo-morph.2019.106957

Rizzoli P and 9 others (2017) Generation and performance assessment of the global TanDEM-X digital elevation model. *ISPRS Journal of Photogrammetry and Remote Sensing* 132, 119–139. doi: 10.1016/j.isprsjprs.2017.08.008

Robert CP and Casella G (1999) The Metropolis-Hastings algorithm. In *Monte Carlo Statistical Methods*, New York, NY: Springer, pp. 231–283.

Samsonov S, Tiamo K and Cassotto R (2021a) 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(9), 4221–4239. doi: 10.5194/tc-15-4221-2021

Samsonov S, Tiamo K and Cassotto R (2021b) SAR-derived flow velocity and its link to glacier surface elevation change and mass balance. *Remote Sensing of Environment* 258 112343. doi: 10.1016/j.rse.2021.112343

Sergienko O (2012) The effects of transverse bed topography variations in ice-flow models. *Journal of Geophysical Research: Earth Surface* 117(F3), F03011. doi: 10.1029/2011JF002203

Sharp M (2021) Amplification of Surface Topography during Surges of Tweedsmuir Glacier. Master's thesis, University of Calgary, Calgary.

Shean DE and 6 others (2016) An automated, open-source pipeline for mass production of digital elevation models (DEM's) from very-high-resolution commercial stereo satellite imagery. *ISPRS Journal of Photogrammetry and Remote Sensing* 116, 101–117. doi: 10.1016/j.isprsjprs.2016.03.012

Tarantola A (2005) Inverse Problem Theory and Methods for Model Parameter Estimation. Philadelphia: SIAM.

Wasserman L (2010) All of Statistics: A Concise Course in Statistical Inference. New York: Springer Science & Business Media.

Werder MA, Huss M, Paul F, Dehecq A and Farinotti D (2020) A Bayesian ice thickness estimation model for large-scale applications. *Journal of Glaciology* 66(255), 137–152. doi: 10.1017/jog.2019.93

Young EM, Flowers GE, Berthier E and Lotto R (2021) An imbalancing act: the delayed dynamic response of the Kaskawulsh Glacier to sustained mass loss. *Journal of Glaciology* 67(262), 313–330. doi: 10.1017/jog.2020.107