Hostname: page-component-848d4c4894-nr4z6 Total loading time: 0 Render date: 2024-05-24T23:53:01.399Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  19 January 2023

Alexi Morin
Affiliation:
Department of Earth Sciences, Simon Fraser University, Burnaby, British Columbia, Canada
Gwenn E. Flowers*
Affiliation:
Department of Earth Sciences, Simon Fraser University, Burnaby, British Columbia, Canada
Andrew Nolan
Affiliation:
Department of Earth Sciences, Simon Fraser University, Burnaby, British Columbia, Canada
Douglas Brinkerhoff
Affiliation:
Department of Computer Science, University of Montana, Missoula, Montana, USA
Etienne Berthier
Affiliation:
LEGOS, Université de Toulouse, CNES, CNRS, IRD, UPS, 31400 Toulouse, France
*
Author for correspondence: Gwenn E. Flowers, E-mail: gflowers@sfu.ca
Rights & Permissions [Opens in a new window]

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.

Type
Letter
Creative Commons
Creative Common License - CCCreative Common License - BYCreative Common License - NCCreative Common License - SA
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.
Copyright
Copyright © The Author(s), 2023. Published by Cambridge University Press on behalf of The International Glaciological Society

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.

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 outlines are shown, along with ice-thickness measurement locations (orange). Flowline (black with km markers) constructed following Kienholz and others (Reference Kienholz, Rich, Arendt and Hock2014) with tools available in the Open Global Glacier Model (Maussion and others, Reference Maussion2019) and manually adjusted to improve agreement with velocity field. Background image from Sentinel-2, 29 August 2018 (UTM Zone 7N).

2. Bayesian inference

We seek to characterize the posterior distribution of model variables

(1)$${\cal P}( {\bf m}\vert {\bf d}) \propto {\cal P}( {\bf d}\vert {\bf m}) \, {\cal P}( {\bf m}) ,\; $$

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

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

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

(3)$${\partial H\over \partial t}( {\bf x}) = {\partial S\over \partial t}( {\bf x}) = - \nabla \cdot \left[\bar{U}( {\bf x}) \, H( {\bf x}) \, {\bf N}( {\bf x}) \right] + \dot{b}( {\bf x}) ,\; $$

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):

(4)$$w( r) \, {\partial S\over \partial t}( r) = - {\partial\over \partial r} \left[w( r) \, \bar{U}( r) \, H( r) \right] + w( r) \, \dot{b}( r) ,\; $$

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

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

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

(6)$$U_{\rm s}( r) = {s\over w( r) \, H( r) }\int_0^r w( r') \left(\dot{b}( r') - {\Delta S( r') \over \Delta t}\right)\, {\rm d}r',\; $$

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:

(7)$$\eqalign{U_{\rm s}( r,\; \Delta t_{\rm m}) \cr& \hskip-3.7pc= {s( \Delta t_{\rm m}) \over w( r) \, ( S( r,\; \Delta t_{\rm m}) - B( r) ) }\int_0^r w( r') \left(\dot{b}( r',\; \Delta t_{\rm m}) - {\Delta S( r',\; \Delta t_{\rm m}) \over \Delta t_{\rm m}}\right)\, {\rm d}r',\;}$$

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 106 iterations each, with the first 105 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:

(8)$${\rm RMSE} = \sqrt{{\sum_{i = 1}^n( B_i - \hat {{B_i}}) ^2\over n}},\; $$

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

(9)$$r( {\delta B},\; \delta \hat{B}) = {{\rm Cov}( {\delta B},\; {\delta\hat{B}}) \over s_{\delta B} s_{\delta\hat{{B}}}},\; $$

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 (xz 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 kth 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:

(10)$$B_{\rm synth} = B_{\rm F2019} + \sum_{k = 5}^{15} A_{ k} \sin{\left({2\pi\over k \tilde{H}} \right)}.$$

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.

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.

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).

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.

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.

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).

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 uncertainty. (a) Surface flow speed. (b) Surface-elevation change rate. Profiles taken along flowline in Figure 1.

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.

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.

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.

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.

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 × 106 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.

Footnotes

*

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

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), 795814. doi: 10.5194/tc-13-795-2019CrossRefGoogle Scholar
Azam, MF and 9 others (2021) Glaciohydrology of the Himalaya-Karakoram. Science 373(6557), 3668. doi: 10.1126/science.abf3668CrossRefGoogle ScholarPubMed
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), 31233148. doi: 10.5194/tc-16-3123-2022CrossRefGoogle Scholar
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), 494507. doi: 10.1017/jog.2019.32CrossRefGoogle Scholar
Brinkerhoff, DJ, Aschwanden, A and Truffer, M (2016) Bayesian inference of subglacial topography using mass conservation. Frontiers in Earth Science 4, 115. doi: 10.3389/feart.2016.00008CrossRefGoogle Scholar
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), 668673. doi: 10.1038/ngeo2999CrossRefGoogle Scholar
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), 12471250. doi: 10.5194/gmd-7-1247-2014CrossRefGoogle Scholar
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), 407417. doi: 10.5194/tc-7-407-2013CrossRefGoogle Scholar
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), 949970. doi: 10.5194/tc-11-949-2017CrossRefGoogle Scholar
Farinotti, D and 6 others (2019) A consensus estimate for the ice thickness distribution of all glaciers on Earth. Nature Geoscience 12(3), 168173. doi: 10.1038/s41561-019-0300-3CrossRefGoogle Scholar
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.571923Google Scholar
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.Google Scholar
Frappé, TP and Clarke, GKC (2007) Slow surge of Trapridge Glacier, Yukon Territory, Canada. Journal of Geophysical Research 112, F03S32. doi: 10.1029/2006JF000607CrossRefGoogle Scholar
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), 168172. doi: 10.1017/jog.2019.2CrossRefGoogle Scholar
Gagliardini, O and 14 others (2013) Capabilities and performance of Elmer/Ice, a new-generation ice sheet model. Geoscientific Model Development 6(4), 12991318. doi: 10.5194/gmd-6-1299-2013CrossRefGoogle Scholar
Gardner, AS and 6 others (2018) Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years. The Cryosphere 12(2), 521547. doi: 10.5194/tc-12-521-2018CrossRefGoogle Scholar
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.Google Scholar
Gelman, A, Carlin, JB, Stern, HS and Rubin, DB (1995) Bayesian Data Analysis. New York: Chapman and Hall/CRC.CrossRefGoogle Scholar
Gilbert, A and 8 others (2018) Mechanisms leading to the 2016 giant twin glacier collapses, Aru Range, Tibet. The Cryosphere 12(9), 28832900. doi: 10.5194/tc-12-2883-2018CrossRefGoogle Scholar
Gudmundsson, GH (2003) Transmission of basal variability to a glacier surface. Journal of Geophysical Research: Solid Earth 108(B5), 2253. doi: 10.1029/2002JB002107CrossRefGoogle Scholar
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, 9196. doi: 10.3189/172756403781816248CrossRefGoogle Scholar
Hugonnet, R and 9 others (2021) Accelerated global glacier mass loss in the early twenty-first century. Nature 592(7856), 726731. doi: 10.1038/s41586-021-03436-zCrossRefGoogle ScholarPubMed
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), 503519. doi: 10.5194/tc-8-503-2014CrossRefGoogle Scholar
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), 49014907. doi: 10.5194/tc-15-4901-2021CrossRefGoogle Scholar
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.Google Scholar
Maussion, F and 9 others (2019) The open global glacier model (OGGM) v1.1. Geoscientific Model Development 12(3), 909931. doi: 10.5194/gmd-12-909-2019CrossRefGoogle Scholar
Meier, MF and Post, A (1969) What are glacier surges? Canadian Journal of Earth Sciences 6(4), 807817. doi: 10.1139/e69-081CrossRefGoogle Scholar
Millan, R, Mouginot, J, Rabatel, A and Morlighem, M (2022) Ice velocity and thickness of the world's glaciers. Nature Geoscience 15(2), 124129. doi: 10.1038/s41561-021-00885-zCrossRefGoogle Scholar
Mingo, L and Flowers, GE (2010) An integrated lightweight ice-penetrating radar system. Journal of Glaciology 56(198), 709714. doi: 10.3189/002214310793146179CrossRefGoogle Scholar
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), 10511063. doi: 10.1017/jog.2020.75CrossRefGoogle Scholar
Pfeffer, WT and 19 others (2014) The Randolph Glacier Inventory: a globally complete inventory of glaciers. Journal of Glaciology 60(221), 537552. doi: 10.3189/2014JoG13J176CrossRefGoogle Scholar
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.geomorph.2019.106957CrossRefGoogle Scholar
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, 119139. doi: 10.1016/j.isprsjprs.2017.08.008CrossRefGoogle Scholar
Robert, CP and Casella, G (1999) The Metropolis–Hastings algorithm. In Monte Carlo Statistical Methods, New York, NY: Springer, pp. 231–283.CrossRefGoogle Scholar
Samsonov, S, Tiampo, 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), 42214239. doi: 10.5194/tc-15-4221-2021CrossRefGoogle Scholar
Samsonov, S, Tiampo, 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.112343CrossRefGoogle Scholar
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/2011JF002203CrossRefGoogle Scholar
Sharp, M (2021) Amplification of Surface Topography during Surges of Tweedsmuir Glacier. Master's thesis, University of Calgary, Calgary.Google Scholar
Shean, DE and 6 others (2016) An automated, open-source pipeline for mass production of digital elevation models (DEMs) from very-high-resolution commercial stereo satellite imagery. ISPRS Journal of Photogrammetry and Remote Sensing 116, 101117. doi: 10.1016/j.isprsjprs.2016.03.012CrossRefGoogle Scholar
Tarantola, A (2005) Inverse Problem Theory and Methods for Model Parameter Estimation. Philadelphia: SIAM.CrossRefGoogle Scholar
Wasserman, L (2010) All of Statistics: A Concise Course in Statistical Inference. New York: Springer Science & Business Media.Google Scholar
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), 137152. doi: 10.1017/jog.2019.93CrossRefGoogle Scholar
Young, EM, Flowers, GE, Berthier, E and Latto, R (2021) An imbalancing act: the delayed dynamic response of the Kaskawulsh Glacier to sustained mass loss. Journal of Glaciology 67(262), 313330. doi: 10.1017/jog.2020.107CrossRefGoogle Scholar
Figure 0

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 outlines 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).

Figure 1

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.

Figure 2

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.

Figure 3

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.

Figure 4

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 uncertainty. (a) Surface flow speed. (b) Surface-elevation change rate. Profiles taken along flowline in Figure 1.

Figure 5

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.

Figure 6

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.

Supplementary material: PDF

Morin et al. supplementary material

Morin et al. supplementary material

Download Morin et al. supplementary material(PDF)
PDF 2.8 MB