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

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.

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

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, 2005). 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 (2016) 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 ₁ − t ₀, 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, 1999), 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⁶ iterations each, with the first 10⁵ 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:

(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, 2010) 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, 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 _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, 2021), and solving the continuity equation. Positive offsets of 1.95–2.1 m a⁻¹ 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 β from 1 to 3.5 × 10⁻⁴ MPa a m⁻¹ 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, 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).

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, 2019). ITS_LIVE velocity products are produced with auto-RIFT offset tracking (Gardner and others, 2018) 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.

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, 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 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, 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 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⁶ 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.