## 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, Reference Gilbert2018; Azam and others, Reference Azam2021; Fox-Kemper and others, Reference Fox-Kemper2021; Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022). 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, Reference Farinotti2017, Reference Farinotti2020, for an overview). Distributed estimates of glacier thickness have now been made on a global scale (e.g. Farinotti and others, Reference Farinotti2019; Millan and others, Reference Millan, Mouginot, Rabatel and Morlighem2022), and are seeing increased use in regional and local studies despite largely unquantified site-specific uncertainties (cf. Pelto and others, Reference Pelto, Maussion, Menounos, Radić and Zeuner2020; Werder and others, Reference Werder, Huss, Paul, Dehecq and Farinotti2020).

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, Reference Main2019; Samsonov and others, Reference Samsonov, Tiampo and Cassotto2021a), and a corresponding need for estimates of ice thickness/bed topography. Theory (e.g. Gudmundsson, Reference Gudmundsson2003) and observation (e.g. Gudmundsson and others, Reference Gudmundsson, Adalgeirsdóttir and Björnsson2003; De Rydt and others, Reference De Rydt, Gudmundsson, Corr and Christoffersen2013; Sharp, Reference Sharp2021) 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 (Reference Brinkerhoff, Aschwanden and Truffer2016) 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

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 ${\cal P}( {\bf m}\vert {\bf d})$, which represents the joint probability distribution over the parameters constrained by the data, as proportional to (∝) the product of the likelihood ${\cal P}( {\bf d}\vert {\bf m})$, which measures the probability of the data under a hypothesized parameter value, and prior ${\cal P}( {\bf m})$, which encodes initial assumptions about the parameters prior to consideration of observations (Tarantola, Reference Tarantola2005). We model the likelihood by assuming that spatially and temporally distributed observations are independent of one another and normally distributed around a mean value that is given directly by the parameters (in the case of the bed elevation and the rate of surface elevation change) or can be computed via the solution of the glaciological 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

with hyper-parameters *σ* (amplitude) and *ℓ* (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 synthetic 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*(**x**) = *S*(**x**) − *H*(**x**), continuity can be written

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

and time-average the forward model over the observation period Δ*t* = *t* _{1} − *t* _{0}, such that

where Δ*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_{\rm s}( r) = s \, \bar {U}( r)$ and solving Eqn (5) for flow speed yields

where the flow-unit width *w*(*r*) as well as parameter *s* are determined 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:

where each term represents an average over an observational epoch Δ*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, Reference Robert and Casella1999), a Markov chain Monte Carlo (MCMC) method, to draw a collection 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, Reference Gelman, Carlin, Stern and Rubin1995). Convergence of the Markov chains is assessed using the Gelman–Rubin statistic. We quantify model performance by computing the RMSE (Chai and Draxler, Reference Chai and Draxler2014) between the true and posterior bed:

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, Reference Wasserman2010) defined as

for the true (*δB*) and inferred ($\delta \hat {B}$) bed perturbations, and their respective sampled standard deviations *s* _{δ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 ($\lesssim \!10 H$) 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, Reference Gagliardini2013) 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, Reference Gagliardini2013). 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 (Reference Farinotti2019), 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* _{F2019}: $B_{{\rm synth}_k} = B_{\rm F2019} + A_{k} \sin {( {2\pi }/{k \tilde {H}} ) }$, where *A* _{k} = *λ* _{k} *R* (m) and *λ* _{k} (m) are the amplitude and wavelength of the *k*th bed perturbation, respectively, *R* = 0.01 is the prescribed amplitude-to-wavelength ratio, $\tilde {H} = 100$ m is the characteristic ice thickness and *k* = [5, 6, …, 15]. We choose the value of *R* such that perturbations are significant in amplitude but do not change the overall shape of the bed (see Supplementary material). The value of $\tilde {H}$ is based on the mean ice thickness of the reference geometry. The minimum value of *k* corresponds to a perturbation wavelength of $k \tilde {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* _{synth} is also defined by adding the sum of all perturbations for *k* = [5, 6, …, 15] to *B* _{F2019}:

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, Reference Young, Flowers, Berthier and Latto2021), 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, Reference Meier and Post1969), 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 *β* from 1 to 3.5 × 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, 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, Reference Young, Flowers, Berthier and Latto2021), 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, Reference Shean2016; Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017). The processing of the DEMs follows the workflow presented in Berthier and Brun (Reference Berthier and Brun2019): briefly, all DEMs are coregistered to the global TanDEM-X DEM (Rizzoli and others, Reference Rizzoli2017) on stable terrain, masking out glaciers using the Randolph Glacier Inventory (RGI) v6.0 (Pfeffer and others, Reference Pfeffer2014). 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, Reference Gardner, Fahnestock and Scambos2019). ITS_LIVE velocity products are produced with auto-RIFT offset tracking (Gardner and others, Reference Gardner2018) using Landsat 4, 5, 7 and 8 imagery. Velocities for the 2016–18 (surge) epoch are taken from Samsonov and others (Reference Samsonov, Tiampo and Cassotto2021a) and derived using synthetic aperture radar methods described in Samsonov and others (Reference Samsonov, Tiampo and Cassotto2021b). A combination of these two datasets are needed, because ITS_LIVE velocities are unreliable during the surge, while the dataset from Samsonov and others (Reference Samsonov, Tiampo and Cassotto2021a) 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, Reference Young, Flowers, Berthier and Latto2021), 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, Reference Mingo and Flowers2010) 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.

## 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 glacier, 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, Reference Sergienko2012). 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, Reference Sergienko2012). The extent to which the slip-induced amplification of these structures differs between 2-D and 3-D models is relevant here: Sergienko (Reference Sergienko2012) 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, Reference Sergienko2012). While our use of 2-D flowband models leads to over-estimation 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, Reference Sharp2021) 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, Reference Frappé and Clarke2007).

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 hyper-parameters 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^{6} 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, Reference Gudmundsson2003). 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, Reference Altena, Scambos, Fahnestock and Kääb2019; Fu and others, Reference Fu, Li and Zhou2019; Rashid and others, Reference Rashid, Majeed, Jan and Glasser2020; Hugonnet and others, Reference Hugonnet2021; Leclercq and others, Reference Leclercq, Kääb and Altena2021; Beaud and others, Reference Beaud, Aati, Delaney, Adhikari and Avouac2022). With this trend towards near-real-time processing of satellite data (e.g. Gardner and others, Reference Gardner, Fahnestock and Scambos2019) 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, Reference Samsonov, Tiampo and Cassotto2021a). 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.