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 ₁ − 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, 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⁶ 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, 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.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 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⁻¹ 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⁻⁴ 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.

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.