Hostname: page-component-848d4c4894-v5vhk Total loading time: 0 Render date: 2024-06-14T02:04:29.921Z Has data issue: false hasContentIssue false
  • English
  • Français

Inferring time-dependent calving dynamics at Helheim Glacier

Published online by Cambridge University Press:  10 August 2022

Jacob Downs Open the ORCID record for Jacob Downs [Opens in a new window] ,
Douglas Brinkerhoff  and
Mathieu Morlighem Open the ORCID record for Mathieu Morlighem [Opens in a new window]
Jacob Downs*
Affiliation:
Department of Computer Science, University of Montana, Missoula, MT, USA
Douglas Brinkerhoff
Affiliation:
Department of Computer Science, University of Montana, Missoula, MT, USA
Mathieu Morlighem
Affiliation:
Department of Earth Sciences, Dartmouth College, Hanover, NH 03755, USA Department of Earth System Science, University of California, Irvine, CA 92697, USA
*
Author for correspondence: Jacob Downs, E-mail: jacob.downs@umontana.edu
Rights & Permissions [Opens in a new window]

Abstract

We perform Bayesian inference of the parameters of a time-dependent model of ice flow and calving at Helheim Glacier, East Greenland. We find that, while a time-independent calving parameterization can recover the mean observed terminus position, such a model is unable to recover sub-annual variability, even when forced with seasonally varying climate. To address this, we develop a simple stochastic model relating surface runoff rates and calving threshold. Again inferring model parameters from observations, we find that this parameterization is able to reproduce observations with respect to both mean position and characteristic temporal variability. This result demonstrates the importance of considering potential sub-annual controls on calving rates in numerical models, which may include variable undercutting rates or surface runoff-dependent surface crevasse propagation.

Keywords

Glacier calvingglacier flowglacier modeling
Type
Article
Information
Journal of Glaciology , Volume 69 , Issue 274 , April 2023 , pp. 381 - 396
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
Copyright © The Author(s), 2022. Published by Cambridge University Press

1. Introduction

Mass loss from the Greenland ice sheet (GrIS) can be partitioned between surface mass balance (SMB) and ice dynamics, with a large percentage ($\approx 66\percnt$ annually) being attributed to ice dynamics (Mouginot and others, Reference Mouginot2019). Warming surface temperatures result in direct mass loss via surface melting, while losses due to ice dynamics occur due to discharge of glacier ice into the ocean, where it is lost to calving and submarine melting. Marine-terminating outlet glaciers drain ~88% of the GrIS by area (Rignot and Mouginot, Reference Rignot and Mouginot2012). Therefore, understanding the dynamics of Greenland's marine-terminating glaciers is critical for estimating future mass loss from the GrIS.

A principal challenge in understanding the dynamics of marine-terminating glaciers is that their behavior is governed by a myriad of interrelated physical processes affected by both marine and atmospheric forcings. Changing ocean circulation patterns have led to increased subaqueous melting of tidewater glaciers, which has been linked to glacier thinning and flow acceleration (e.g. Holland and others, Reference Holland, Thomas, de Young, Ribergaard and Lyberth2008; Motyka and others, Reference Motyka2011; Rignot and others, Reference Rignot, Fenty, Menemenlis and Xu2012). In addition to causing direct mass loss, subaqueous melting impacts ice dynamics by undercutting calving fronts, potentially amplifying calving rates (Fried and others, Reference Fried2019; Slater and others, Reference Slater, Benn, Cowton, Bassis and Todd2021).

Frontal ablation from calving and submarine melting is affected not only by the marine environment, including ocean thermal forcing and salinity, but also by surface processes. Meltwater pooling in crevasses may result in fracture propagation and thereby cause iceberg calving (van der Veen, Reference van der Veen1998). Seasonal variations in meltwater input to the glacier bed cause changes in subglacial water pressure and basal sliding speed (e.g. Sundal and others, Reference Sundal2011). While seasonal acceleration is more pronounced on land-terminating outlet glaciers than marine-terminating glaciers (Davison and others, Reference Davison2020), seasonal velocity variations may still have an important impact on calving (Cook and others, Reference Cook2014). More importantly, surface runoff is linked to submarine melt rates via subglacial discharge (Xu and others, Reference Xu, Rignot, Fenty, Menemenlis and Flexas2013). Hence, greater subglacial discharge results in greater subaqueous melting and undercutting.

The multitude of interrelated physical processes acting on marine-terminating glaciers makes modeling their dynamics challenging and obscures the primary causes of mass loss and retreat. Recent modeling studies have made strides in coupling models of ice-dynamics, subglacial hydrology, frontal plume melting and calving (e.g. Cook and others, Reference Cook, Christoffersen and Todd2021). However, there are still significant uncertainties in state of the art models. Questions remain as to which sliding laws and calving laws are most appropriate in numerical models (Amaral and others, Reference Amaral, Bartholomaus and Enderlin2020; Åkesson and others, Reference Åkesson, Morlighem, O'Regan and Jakobsson2021), and work on understanding these processes is an active area of research. Even if a specific model formulation is assumed correct, effort is needed to rigorously constrain unknown model inputs and parameters.

In this work, we aim to understand the primary physical processes controlling seasonal terminus position variability at Helheim Glacier from 2007 to 2019. Helheim Glacier is one of the largest outlet glaciers in Greenland in terms of ice-flux (Mankoff and others, Reference Mankoff2019). Helheim has retreated ~6 km since 2003 (Fig. 1), but this retreat has not been linear. Observations indicate a seasonal signal in terminus position with advance during winter and retreat in summer (Bevan and others, Reference Bevan, Luckman, Khan and Murray2015; Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017). Processes driving this seasonal cycle are uncertain, but modeling indicates that Helheim may be sensitive to factors such as the depth of water in surface crevasses, which could influence calving rates on a seasonal time scale, and subglacial water pressure, which could cause seasonal acceleration (Cook and others, Reference Cook2014). Evidence based on geometry and nearby glaciers also suggests that Helheim may soon undergo substantial retreat (Williams and others, Reference Williams, Gourmelen, Nienow, Bunce and Slater2021).

Fig. 1. Terminus position projected along a central flowline (a) and surface runoff (b) at Helheim glacier between 2003 and 2019. Shaded blue regions indicate periods of terminus advance.

To better understand the physical processes controlling terminus position variability at Helheim, we model its dynamics between 2007 and 2019 using the Ice-sheet and Sea-level System Model (ISSM), a finite element based, thermo-mechanically coupled ice-sheet model (Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012). To identify the most important physical processes affecting terminus position, we first estimate a number of time-independent model parameters and their uncertainties by comparing modeled and observed terminus positions. We find that static parameters do not capture observed seasonal fluctuations in terminus position. To rectify this issue, we perform sensitivity experiments using seasonally varying model forcings. Based on the results of these sensitivity experiments, we develop and test a stochastic model that relates the calving stress threshold in a von-Mises calving law (Morlighem and others, Reference Morlighem2016) to surface runoff at Helheim, such that the resulting instances of modeled terminus position exhibit statistical agreement with observations of front position.

We structure the manuscript in an unconventional manner, such that results of a given experiment are presented before describing subsequent experiments. This presentation is intended to help motivate each experiment given previous findings. In Section 2, we present the numerical model and estimate a number of time-independent model parameters and present the results. Motivated by the time-independent parameter estimates, Section 3 focuses on experiments using seasonal forcings. Finally, building on all previous results, we present a Markov model for estimating time-dependent calving rates in Section 4.

2. Methods

2.1. Ice dynamics model

Ice dynamics at Helheim are modeled using ISSM (Larour and others, Reference Larour, Seroussi, Morlighem and Rignot2012). Ice velocity is modeled using the shallow shelf approximation (MacAyeal, Reference MacAyeal1989). Dirichlet boundary conditions on velocity are imposed on the interior and lateral boundaries of the computational domain using a composite ice velocity map from 1985 to 2018 from the ITS_LIVE project (Gardner and others, Reference Gardner, Fahnestock and Scambos2021a). Lateral boundaries are determined by flowline velocity tracing. We use an anisotropic computational mesh refined based on the ITS_LIVE velocity map with a maximum resolution of ~5 km in the interior and ~200 m in the fast flowing regions of Helheim and in the fjord. Sliding follows a linear, Budd-like, friction law of the form

(1)$$\tau_{\rm b} = \beta^2 N {\bf u}$$

where τb is the basal stress, N is effective pressure, u is basal ice velocity and $\beta _0^2$ is an uncertain traction coefficient (Budd and others, Reference Budd, Keage and Blundy1979). Effective pressure is given by

(2)$$N = P_{\rm i} - P_{\rm w}$$

where P i is ice overburden pressure and P w is water pressure. Radar interferometry at Helheim indicates a relatively small tidal influence on near terminus velocity, suggesting a strong coupling to the bed (Voytenko and others, Reference Voytenko2015). We therefore assume that water pressure is a fixed fraction (85%) of overburden pressure (Wright and others, Reference Wright, Harper, Humphrey and Meierbachtol2016), and that there is not a strong hydraulic connection to the ocean. However, we further examine the role of subglacial water pressure in Section 3. Bed and surface geometries for Helheim are based on BedMachine v3 (Morlighem and others, Reference Morlighem2017). Surface elevation from BedMachine yields a nominal ice-sheet geometry corresponding to the year 2007, which we select as the starting point for all model runs.

We adopt a von Mises calving law in ISSM (Morlighem and others, Reference Morlighem2016) to determine the calving rate at the ice-front position, which is tracked using a levelset method (Bondzio and others, Reference Bondzio2016). The calving rate c is given by

(3)$$c = \Vert{{\bf u}}\Vert {\tilde{\sigma}\over \sigma}.$$

where σ is an uncertain stress threshold and $\tilde {\sigma }$ is the tensile von Mises stress. The latter is defined as

(4)$$\tilde{\sigma} = \sqrt{3} B {\tilde{\dot{\varepsilon}}_{\rm e}}^{1/3}$$

where B is an uncertain, temperature-dependent ice hardness. The effective tensile strain rate $\tilde {\dot {\varepsilon }}_{\rm e}$ is defined by

(5)$$\tilde{\dot{\varepsilon}}_{\rm e}^2 = {1\over 2} \left(\max ( 0,\; \dot{\varepsilon_1}) ^2 + \max ( 0,\; \dot{\varepsilon_2}) ^2 \right)$$

where $\dot {\varepsilon _1}$ and $\dot {\varepsilon _2}$ represent the eigenvalues of the 2D horizontal strain rate vector. At the calving front, a water pressure boundary condition is applied automatically where the bed is below sea level and a zero stress boundary condition is used otherwise.

2.2. Estimation of time-independent parameters

In order to better understand the physical processes controlling terminus position variability at Helheim, we seek probability distributions for model parameters that best explain observed changes in terminus position at Helheim from 2007 to 2019. We first examine several static (time-independent) model parameters. Here, we define prior distributions that we use to sample uncertainty in forcings including SMB, boundary conditions (basal friction coefficient) and model parameters (calving stress threshold, submarine melting and ice temperature). Prior distributions are intended to be very broad in order to define a large space of admissible parameter combinations. Prior distributions are randomly sampled to generate training data for a surrogate model, or emulator, that approximates the output of ISSM as outlined in Section 2.3. Using this surrogate model, refined parameter estimates are obtained by incorporating terminus position observations as described in Section 2.5.

2.2.1. Surface mass balance

We use seasonally varying SMB $\dot {a}_0( x,\; \, y,\; \, i,\; \, j)$ from the Modèle Atmosphérique Régional (MAR) indexed by year i and month j (Tedesco and Fettweis, Reference Tedesco and Fettweis2020). One way to parameterize the uncertainty in SMB would be to introduce a simple scaling parameter that multiples $\dot {a}_0( x,\; \, y,\; \, i,\; \, j)$ field. A limitation of this approach is that it would not account for spatial variability in uncertainty in the SMB field. For example, seasonal variability is greater near the terminus than in the interior. To address this, we parameterize uncertainty in SMB in a way that captures its natural spatial and temporal variability. We introduce an SMB bias parameter $p_{\dot {a}}$ and let

(6)$$\dot{a}( x,\; y,\; i,\; j) = \left\{\matrix{ \vert{\,p_{\dot{a}}\vert} \dot{a}_{\min}( x,\; y,\; j) \hfill & p_{\dot{a}} < -1 \hfill \cr \vert{\,p_{\dot{a}}\vert} \dot{a}_{\min}( x,\; y,\; j) + ( 1 - \vert{\,p_{\dot{a}}\vert}) \dot{a}_0( x,\; y,\; i,\; j) \hfill & -1 \leq p_{\dot{a}} \leq 0 \hfill \cr \vert{\,p_{\dot{a}}\vert} \dot{a}_{\max}( x,\; y,\; j) + ( 1 - \vert{\,p_{\dot{a}}\vert}) \dot{a}_0( x,\; y,\; i,\; j) \hfill & 0 < p_{\dot{a}} \leq 1 \hfill \cr \vert{\,p_{\dot{a}}\vert} \dot{a}_{\max}( x,\; y,\; j) \hfill & p_{\dot{a}} > 1 \hfill \cr \hfill }\right..$$

Here $-1.33 \leq p_{\dot {a}} \leq 1.33$, and a min(x,  y,  j) and a max(x,  y,  j) represent the minimum and maximum SMB values at a given location and month across all years between 2007 and 2019. In effect, $p_{\dot {a}}$ values greater than zero indicate seasonally high SMB, whereas values less than zero indicate seasonally low values. A value of zero recovers the default MAR SMB (Fig. 2). Values of $p_{\dot {a}}$ above 1 or below − 1 allow us to explore sensitivity to SMB by allowing for broader spectrum of seasonal mass-balance values slightly beyond the minimum and maximum seasonal values in MAR.

Fig. 2. Idealized, conceptual example of summer SMB at a given point in time along a flowline for different values of $p_{\dot {a}}$. Values of − 1, 0 and 1 map to $\dot {a}_{\min }$, $\dot {a}_0$ and $\dot {a}_{\max }$.

2.2.2. Basal friction

We estimate a baseline basal traction coefficient $\beta _0^2$ by inverting for surface velocity (Morlighem and others, Reference Morlighem2010). To maximize spatial coverage of velocity observations, we again use the 1985–2018 composite ice velocity field from ITS_LIVE (Gardner and others, Reference Gardner, Fahnestock and Scambos2021a). The basal traction coefficient β 2 is then given by

(7)$$\beta^2 = p_{\beta_0} \beta_0^2.$$

Uncertainty in β 2 is sampled by allowing $p_{\beta _0}$ to vary between 0.66 and 1.33. The parameter $p_{\beta _0}$ therefore uniformly scales the traction field from 66 to 133% of its default value. This range is slightly broader than that tested in Downs and Johnson (Reference Downs and Johnson2022), where a 50% change in the basal traction coefficient was found to cause a peak change 50% in near terminus velocity using a similar model formulation on a marine-terminating glacier. We adopt a slightly larger range than in previous work so as to provide sufficient training data for the surrogate model (Section 2.3).

2.2.3. Calving stress threshold

The calving stress threshold σ used in the von Mises calving law is allowed to vary in the interval 25 × 103 ≤ σ ≤ 2000 × 103 Pa. A lower threshold of 25 × 103 Pa represents a physical extreme, where ice calves nearly immediately when entering the sea. The upper threshold of 2000 × 103 Pa represents the opposite physical extreme in which ice rarely calves. Note that calving in ISSM occurs only for ice in contact with the ocean.

2.2.4. Submarine melting

Since the ice front of Helheim is lightly grounded (Melton and others, Reference Melton2022), we choose to impose a frontal submarine melting or undercutting rate as opposed to modeling melt on the underside of floating ice. Enderlin and Howat (Reference Enderlin and Howat2013) estimated a mean submarine melt rate, $\dot {m}$, of ~0.56 m d−1 at Helheim between 2000 and 2010. Rignot and others (Reference Rignot, Koppes and Velicogna2010) estimated melt rates of between 1 and 4 m d−1 at three grounded tidewater glaciers in West Greenland. Based on this range of estimates for marine-terminating glaciers in Greenland, we constrain the submarine melting rate $\dot {m}$ between 0 and 4 m d−1. As with other parameters, we use a broad range of plausible melt values for the prior distribution as opposed to a specific range estimated for Helheim.

2.2.5. Ice temperature

To estimate a baseline ice temperature field T 0, we solve for vertically averaged steady-state ice temperature. As boundary conditions, we use a 30-year averaged surface temperature field from 1980 to 2010 from Box (Reference Box2013), and a uniform geothermal heat flux of 70 mW m−2 at the basal boundary based on Martos and 5 others (Reference Martos2018). Using a built-in steady-state solver in ISSM, we then solve for the steady-state 3D ice temperature and flow fields. The resulting temperature field is vertically averaged for use in the vertically integrated stress balance model. We explore uncertainty in the resulting temperature approximation by parameterizing ice temperature as follows:

(8)$$T = \min \left(0,\; T_0 + \Delta T \right).$$

Here ΔT is a spatially uniform temperature offset with − 10°C ≤ 10°C.

2.3. Surrogate model: mapping parameters to model misfit

To facilitate robust parameter estimation via Markov Chain Monte Carlo (MCMC) sampling, we train a Gaussian Process (GP, Rasmussen and Williams, Reference Rasmussen and Williams2006) surrogate model, or emulator, to approximate the output of ISSM for arbitrary parameter inputs. We can think of ISSM as a function that maps a set of parameter inputs to some output of interest. In our case, we are interested in estimating the misfit between modeled and observed terminus positions at Helheim from 2007 to 2019 using ice front positions from Cheng and others (Reference Cheng2021). That is, ISSM can be used to calculate the following residual, which we call the ‘true’ model residual

(9)$$y = {\cal F}( {\bf x}) = {1\over 12} \int_{2007}^{2019} \left(\ell( t) - \ell_o( t) \right)^2 {\rm d}t,\; $$

where $\pmb {x}$ is a vector of parameters and y is a metric measuring the misfit between the modeled and observed terminus positions projected along a central flowline (Fig. 3), ℓ(t) and ℓo(t) respectively, and

(10)$${\bf x} = \left[\,p_{\dot{a}},\; p_{\beta_0},\; \sigma,\; \dot{m},\; \Delta T\right]^T.$$

Since ISSM is expensive for the purposes of parameter estimation using MCMC, we train a computationally efficient surrogate model ${\cal G}( \pmb {x})$ using the GPyTorch module for Python (Gardner and others, Reference Gardner, Pleiss, Bindel, Weinberger and Wilson2021) to approximate the true model residual ${\cal F}( \pmb {x})$. The surrogate model maps a set of input parameters to a Gaussian distribution

(11)$$P( y\vert \pmb{x}) = {\cal G}( \pmb{x}) \sim {\cal N}( m( \pmb{x}) ,\; s^2( \pmb{x}) ) $$

where $m( \pmb {x})$ refers to the mean estimated model residual and $s^2( \pmb {x})$ is the estimated variance for a particular input $\pmb {x}$.

Fig. 3. Location map, bedrock elevation and composite ice velocity map from 1980 to 2018 for Helheim glacier in East Greenland. The cyan line shows the central flowline along which terminus position is projected, while the cyan dot is used for point measurements in certain numerical experiments.

As training data, we perform 3000 evaluations of the full model using Latin-hypercube sampling (Iman, Reference Iman2008) to span the space of admissible parameter combinations. In particular, for the purposes of training the surrogate model, we sample from the prior distribution, which we take to be a Cartesian product of uniform distributions over individual parameters, with bounds given as in Section 2.2. One can think of the surrogate as interpolating the full model between training data points. Additional details on the GP surrogate are provided in Appendix A.

2.4. Likelihood

It is a much easier problem to predict model error than directly predicting terminus position in the surrogate model. However, this means we have to define a sensible likelihood over model error. As a heuristic for defining the likelihood, we seek model parameters that yield a modeled terminus position with an average misfit of ~1 km or less. Therefore we define the ‘observed’ error y o as follows

(12)$$y_o \sim {\cal N}( 0,\; \sigma_o^2) $$

where σo = 1 km. The standard deviation σo should not be confused with the calving stress threshold σ. This yields a Gaussian likelihood of the form

(13)$$P( y_o\vert y) = {1\over \sigma \sqrt{2 \pi}} \exp \left({-y^2\over 2 \sigma_o^2} \right).$$

2.5. Posterior distribution

The surrogate model returns an estimate of the model misfit with an associated uncertainty. For a given point $\pmb {x}$ in parameter space, the misfit $y = {\cal G}( \pmb {x})$ is an imprecise estimate of the true model residual. Intuitively, when sampling from the posterior distribution, uncertainty in the surrogate model contributes to greater uncertainty in parameter estimates. Therefore, uncertainty in the GP surrogate should be incorporated into estimates of the posterior distribution. In order to obtain the posterior distribution over parameters $\pmb {x}$, we first define a joint distribution over inputs $\pmb {x}$ and outputs y

(14)$$P( \pmb{x},\; y \vert y_o) \propto P( \pmb{x}) P( y\vert \pmb{x}) P( y_o \vert y) .$$

Here $P( \pmb {x},\; \, y \vert y_o)$ is the probability of a given input–output pair given the observation y o, $P( \pmb {x})$ is the prior distribution, and $P( y\vert \pmb {x})$ is the probability that the surrogate model predicts y for a given input $\pmb {x}$. To obtain the posterior distribution over $\pmb {x}$, we marginalize over y

(15)$$P( \pmb{x} \vert y_o) \propto \int_y P( \pmb{x}) P( y\vert \pmb{x}) P( y_o \vert y) \; dy.$$

When the likelihood and surrogate prediction are both Gaussian distributed, this integral can be computed analytically (Appendix B).

2.6. Static parameter estimation: results

Estimates of the posterior distribution demonstrate that modeled terminus position at Helheim is most sensitive to parameters affecting frontal ablation, particularly the calving stress threshold, σ, in the von Mises calving law. We achieve close agreement between modeled and observed terminus positions for values of σ in a range between 200 and 600 kPa, with a peak probability at 400 kPa (Fig. 4o).

Fig. 4. Estimates of the posterior distribution obtained by running MCMC sampling using a Gaussian Process surrogate model fitted to an ensemble using static or time-independent parameters. Diagonal subplots a, f, j, m and o show marginal distributions for individual variables obtained using either the most probable surrogate model function (blue line) or accounting for surrogate model uncertainty (red line). Subplots below the diagonal show pairwise marginal distributions accounting for surrogate model uncertainty. Corresponding subplots above the diagonal show the same pairwise marginal distributions estimated using the Gaussian Process mean function (e.g. a.2 shows the same marginal distribution as a.1). Red dots in subplots above the diagonal show 50 ensemble members sub-sampled from the prior based on their posterior probability and plotted in Figure 5.

Fig. 5. Observed terminus position at Helheim from 2003 to 2019 is indicated by the thick red line. Modeled terminus positions for all prior ensemble members are shown as faint black lines. Thin multi-colored lines show terminus positions for an approximation of the posterior distribution, which we obtain from the prior ensemble by weighting the prior ensemble members according to their probability and sub-sampling 50 ensemble members.

The subaqueous melt rate, $\dot {m}$, has a multimodal marginal distribution with one cluster centered at ~0.5 m d−1 and another at ~2 m d−1, but melt rates below ~1 m d−1 are most probable (Fig. 4a). Temperature offsets above 5° C are most probable (Fig. 4f). The sensitivity of the model to ice temperature may be partially due to its influence on ice flow. However, we find that modeled terminus position is fairly insensitive to another ice flow parameter, p β, which rescales the basal traction coefficient obtained by inversion for surface velocities (Fig. 4m). This suggests that the model's sensitivity to the ice temperature may primarily be due to the temperature dependence of the von Mises calving law (see Eqn (3)). In particular, softer, warmer ice corresponds to lower calving rates given the same strain rate and calving stress threshold. We note that inversions for the basal traction coefficient, $\beta _0^2$, are not performed for each temperature offset ΔT. Hence, strain fields differ somewhat between different temperature offsets. However, in the sliding dominant flow regime near the Helheim terminus, we expect that the impact of ice temperature on flow is of secondary importance.

The marginal distribution for the SMB bias parameter $p_{\dot {a}}$ closely resembles the uniform prior distribution, meaning that SMB has a negligible influence on modeled terminus positions on the decadal time scale of interest (Fig. 5j). For the most part, there are not strong correlations between variables. Exceptions include a positive correlation between the undercutting rate $\dot {m}$ and temperature offset ΔT for values of $\dot {m}$ from 0 to 1 m d−1 (Fig. 4b.1), and a negative correlation between $\dot {m}$ and the basal traction scaling parameter p β for $\dot {m} > 1$ m d−1 (Fig. 4d.1).

To understand the impact of surrogate model uncertainty, MCMC sampling is performed using two different strategies. In the first, we account for uncertainty in the GP surrogate model as outlined in Section 2.5. In the second, we neglect uncertainty in the surrogate model and use the most probable or mean GP surrogate model function when sampling the posterior (Fig. 4 all subplots above the diagonal). Accounting for surrogate uncertainty generally leads to broader distributions than using the mean function, but they are otherwise almost identical, suggesting that surrogate model uncertainty does not significantly impact parameter estimates. By validating the surrogate model against test data not included in training, we find that the surrogate reliably estimates the true model misfit within the 2σ uncertainty bounds.

Without performing additional model runs, we sanity check the posterior distribution by weighting prior ensemble members according to their probability and sub-sampling 50 of the prior ensemble members to compare to observations (Fig. 5). While these high probability ensemble members match the observed terminus position within a few kilometers, we find that they do not reproduce observed short-term variations in terminus position. Indeed, beyond imposing a seasonally varying SMB field, we have not introduced any other time-dependent inputs that would be likely to cause such short-term variations. This is addressed in the following sections.

3. Experiments using seasonal forcings

Experiments using time-independent parameter values help us identify the most important parameters controlling terminus position at Helheim. However, they do not allow us to capture its short-term variability. To address this, we test the model's sensitivity to a number of seasonally varying forcings.

3.1. Frontal melt rate

The average subaqueous melt rate depends on both ocean thermal forcing and subglacial freshwater flux. Xu and others (Reference Xu, Rignot, Fenty, Menemenlis and Flexas2013) parameterize the average submarine melt rate as

(16)$$\dot{m} = ( A q_{{\rm sg}}^{\alpha} + B) T_{\rm h}^{\beta}$$

where q sg is the subglacial water flux, T h is ocean thermal forcing and A, B, α, β are uncertain parameters. Equation (16) is based on an empirical fit to a three-dimensional simulation of ice melting and turbulent upwelling from a freshwater plume. Since subglacial discharge is unknown and we do not model subglacial hydrology, we instead parameterize subaqueous melt as a function of surface runoff q s(t). Substituting surface runoff for freshwater flux we have

(17)$$\dot{m} = A T_{\rm h}^{\beta} \left({q_{{\rm s}}( t) \over \max( q_{\rm s}( t) ) } \max( q_{\rm s}( t) ) \right)^{\alpha} + B T_{\rm h}^{\beta}$$

where max(q s(t)) refers to the maximum runoff over the simulated time period. We regroup uncertain terms, letting $\dot {m}_{\max }-\dot {m}_{\min } = A T_{\rm h}^{\beta } \max ( q_{\rm s}) ^{\alpha }$ and $\dot {m}_{\min } = B T_{\rm h}^{\beta }$ yielding

(18)$$\dot{m} = \left({q_{{\rm s}}( t) \over \max( q_{{\rm s}}( t) ) } \right)^{\alpha} ( \dot{m}_{\max} - \dot{m}_{\min}) + \dot{m}_{\min}.$$

Here, $\dot {m}_{\max }$ and $\dot {m}_{\min }$ are uncertain parameters representing the maximum and minimum seasonal subaqueous melt rates respectively, and we let α = 0.5 (Xu and others, Reference Xu, Rignot, Fenty, Menemenlis and Flexas2013). Surface runoff for the Helheim catchment area is estimated using Mankoff and others (Reference Mankoff2020).

To test the sensitivity of terminus position to the subaqueous melt rate, we perform a set of 49 model runs using $\dot {m}_{\max }$ values ranging from 0.5 to 3.5 m d−1. For each value of $\dot {m}_{\max }$, the minimum seasonal melt value $\dot {m}_{\min }$ is obtained by multiplying $\dot {m}_{\max }$ by values between 0 and 1. Other parameter values are determined based on the static parameter estimates already obtained.

3.2. Water pressure

Moon and others (Reference Moon2014) classify Helheim as a type three glacier characterized by acceleration during the early melt season followed by a subsequent slowdown later in the summer. Acceleration during the early melt season is associated with increased input of water into the subglacial drainage system and may result from enhanced sliding caused by high basal water pressure. Velocity declines later in the melt season, ostensibly as the subglacial drainage system adapts to high meltwater input, then gradually rises again during the winter (Davison and others, Reference Davison2020). To simulate this seasonal cycle, we parameterize subglacial water pressure, expressed as a fraction of overburden pressure, by

(19)$$P_{{\rm frac}} = 0.85 + \left(\,p_{_{P_{\rm w}}} \sin \left[-2 \pi t - 2 + {1\over 2} \sin( -2 \pi t - 2) \right]\right)s( x,\; y,\; t) .$$

Here, $0 \leq p_{_{P_{\rm w}}} < 0.15$ is an uncertain parameter determining the magnitude of seasonal pressure oscillations, t is time in years and s(x,  y,  t) is an indicator function that is one when and where surface melt occurs and zero otherwise:

(20)$$s( x,\; y,\; t) = {1\over 1 + \exp( 5 \; \dot{a}( x,\; y,\; t) ) .}$$

Note that in winter, if SMB is uniformly positive, this scaling function is near zero. During the melt season, s(x,  y,  t) ≈ 1 in regions that are melting. Thus, this indicator function confines water pressure oscillations to the ablation zone of Helheim. Other parameter values are determined based on existing static parameter value estimates.

3.3. Seasonal forcings: results

Modeled velocity is sensitive to the mean undercutting rate, yet does not display a seasonal response to variable melt rates. For a high mean seasonal melt rate of 3.5 m d−1, near terminus ice accelerates by ~300% over a matter of ~4 years (Fig. 6b). For mean melt rates below 2 m d−1, velocity remains steady throughout the simulation, whereas melt rates above 2 md−1 cause acceleration of at least 200% for near terminus velocity. The timing of acceleration is controlled by the magnitude of melt (Fig. 6b). In particular, larger mean melt rates correspond to peak acceleration earlier in the simulation.

Fig. 6. (a) Modeled terminus position at Helheim between 2007 and 2017 for experiments using seasonally varying undercutting rates. Line colors correspond to the mean undercutting rate over the duration of the simulation. (b) Percentage change in velocity at the cyan point marked in Figure 3a, relative to initial velocity. (c) Seasonal undercutting rates for a subset of 7 out of the 49 simulations displayed in panels (a) and (b).

Although the model is sensitive to the undercutting rate, seasonal variations in undercutting do not yield discernible seasonal cycles in terminus position (Fig. 6a). A mean melt rate of 3.5 m d−1 ($\dot {m}_{\max } = \dot {m}_{\min } = 3.5$ m d−1) yields a retreat of 8 km over the period from 2007 to 2017, whereas an undercutting rate of 0 m d−1 ($\dot {m}_{\max } = \dot {m}_{\min } = 0$ m d−1) yields a 6 km advance. Acceleration corresponding to higher melt rates is likely a response to terminus retreat. Yet large seasonal variations in undercutting of several meters per day of melt do not reproduce characteristic seasonal oscillations.

Seasonal variations in water pressure result in velocity variations qualitatively similar to near-terminus velocity observations on tidewater glaciers in Greenland (Davison and others, Reference Davison2020). Velocity begins to rise late in the melt season and continues rising through winter, into the following spring (Fig. 7b). This is followed by a sharp decline in velocity during the next melt season.

Fig. 7. (a) Modeled terminus position at Helheim from 2007 to 2017 for runs using different magnitudes of seasonal water pressure oscillations. (b) Seasonal speed anomaly at the cyan point marked in Figure 3. The seasonal speed anomaly is computed by dividing the velocity at a given time by the velocity in a baseline run with no water pressure oscillations. (c) Seasonal water pressure oscillations expressed as a fraction of overburden.

To evaluate the influence of seasonal water pressure variations, we compare near terminus velocities for runs with different magnitudes of seasonal water pressure oscillations to a baseline run wherein water pressure is a fixed fraction of 85% of overburden pressure. The degree of flow acceleration depends on the magnitude of pressure oscillations (Fig. 7). For a seasonal water pressure oscillation of 13% of overburden pressure, in which water pressure peaks at 98% of overburden pressure during the melt season, near terminus velocity at Helheim is up to 225% higher than in the baseline run with no seasonal water pressure variation.

For the largest tested water pressure oscillations of ~13% of overburden pressure, there is an up to a twofold increase in near terminus velocity over the course of the year. Acceleration of this magnitude is inconsistent with observations at Helheim, where seasonal accelerations of up to 50% are more realistic (Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017). However, lower amplitude pressure oscillations of ~8% of overburden pressure yield a seasonal acceleration more consistent with observations. Despite large, near-terminus variations in velocity, there is not a clear seasonal signal in terminus position (Fig. 7a), indicating that pressure variations may not be responsible for observed seasonal variations in terminus position at Helheim.

4. A Markov model for time-dependent calving

There is a strong correlation between observed terminus advance/retreat and surface meltwater runoff at Helheim glacier between 2003 and 2019 (Fig. 1). Generally speaking, high runoff is associated with terminus retreat, while low runoff is associated with terminus advance. As we demonstrated in Section 2.2, modeled terminus position at Helheim is most sensitive to the calving stress threshold σ in the Von Mises Calving law. Hence, we hypothesize that modeled advance or retreat is primarily a function of the imposed stress threshold σ in the model, which directly affects the calving rate.

Motivated by these observations, we define a Markov model that relates σ to surface runoff in a stochastic sense. The Markov model will output one of two states on a biweekly time step. These states are (1) a state in which Helheim is more susceptible to calving for which σ = σmin or (2) a state in which Helheim is less susceptible to calving for which σ = σmax. The calving stress thresholds σmax ≥ σmin are tuned so that if σ = σmin, modeled terminus position at Helheim is likely to retreat due to high calving rates. Similarly, if σ = σmax, modeled terminus position is more likely to advance or remain stagnant. Determination of these calving stress threshold values is discussed in subsequent sections. Given the current state of the terminus (more or less susceptible to calving) and the surface runoff, the aim of the model is to predict whether the glacier will be more or less susceptible to calving the following week.

Conceptually, we hypothesize a link between surface runoff and calving rates due to one of a number of possible physical process, which we address in greater detail in Section 5. For example, one hypothesis is that surface runoff could lead to deeper surface crevasses via hydro fracturing, reducing resistance to tensile stress. This hypothesis was previously proposed and investigated by Cook and others (Reference Cook2014), who found that a surface crevasse-based calving law was sensitive to the modeled water depth in crevasses.

Suppose that S i ∈ {σ = σmin,  σ = σmax} represents the state at time i. Then the Markov model ${\cal M}$ maps the state S i to a new state at time i + 1

(21)$$S_{i + 1} = {\cal M}( S_i,\; r_i) $$

where r i is the surface runoff input (Fig. 8). Transition probabilities from one state to the next, denoted by P(S i+1|S i,  r i), depend on the previous state as well as biweekly surface runoff input. Since surface runoff is a continuous variable, and we have limited observational data, we bin surface runoff data into K = 10 bins B 1,  B 2, · · · ,  B K and seek the following transition probabilities

(22)$$P( S_{i + 1} \vert S_i,\; r_i \in B_k) $$

with S i,  S i+1 ∈ {σ = σmin,  σ = σmax} and 1 ≤ k ≤ N. Transition probabilities are determined using maximum likelihood estimation, wherein transition probabilities are selected to maximize the probability of generating the observed sequences of states for each runoff bin (Appendix C). When determining the transition probabilities, which we also refer to as training the Markov model, we cannot directly observe σ, and instead rely on observations of terminus position assuming that observed terminus retreat corresponds to a lower calving stress threshold state (σ = σmin), and observed terminus advance corresponds to a higher calving stress threshold (σ = σmax). Put otherwise, terminus position observations are converted into a sequence of advancing/retreating states, which are assumed to correspond directly to high/low stress threshold states respectively. Calving stress threshold parameters σmin and σmax need to be estimated to produce realistic results, as the Markov model only determines the transition between states but not what values of σmin and σmax best reproduce observations. We discuss how these parameters are estimated in the next section.

Fig. 8. The Markov model outputs a calving stress threshold on a biweekly time step given calving stress threshold from the previous time step as well as surface runoff. Given the state at time i, denoted S i and predicts the state at the next time step S i+1 given the surface runoff r i.

4.1. Dynamic parameter estimates

We estimate a number of parameters that control time varying model forcings using the same methodology as in Section 2.2. As before, we perform 3000 model runs, which are used to train a surrogate model that approximates the true model residual given arbitrary parameter combinations as input. The surrogate model is subsequently used to sample from the posterior distribution, which allows us to identify parameter combinations that yield a good fit to observations of terminus position.

Our two primary parameters of interest are the calving stress thresholds σmin and σmax, which the Markov model outputs on a biweekly time step. Each ensemble member is forced with different randomly sampled values of these stress threshold parameters as well as a different realization of the Markov model given the same observed surface runoff input data. The parameters σmin and σmax therefore control maximum and minimum calving rates, respectively, while the Markov model determines the seasonality of calving in a stochastic sense by controlling transitions between the high and low calving stress threshold states.

In practice, for each ensemble member, σmax is drawn from a uniform prior distribution with 25 × 103 ≤ σmax ≤ 2000 × 103 Pa. Then, we let $\sigma _{\min } = \sigma _{\max } p_{\sigma _{\min }}$ where the parameter $0 \leq p_{\sigma _{\min }} \leq 1$ is also drawn from a uniform prior distribution. This sampling strategy enforces the constraint σmin ≤ σmax. Once calving parameters have been randomly sampled from the prior distribution at the beginning of each simulation, the Markov model is used to transition between states on a biweekly time step. Thus, once values of σmin and σmax are sampled, the Markov model yields a full time-dependent calving stress threshold input that is used to force ISSM.

Previously, we saw that the model is somewhat sensitive to the subaqueous melt rate. Hence, we let $\dot {m}$ vary seasonally using Eqn (18), and we seek to estimate the minimum and maximum seasonal undercutting rates $\dot {m}_{\min }$ and $\dot {m}_{\max }$ respectively in conjunction with calving parameters. Finally, since static parameter estimates indicated that modeled terminus position was also somewhat sensitive to ice temperature, we estimate the ice temperature offset parameter ΔT. The parameters $\dot {m}_{\min },\; \, \dot {m}_{\max }$ and ΔT are also drawn from uniform prior distributions using identical parameter ranges used for prior distributions for static parameters.

4.2. Markov model: results

Forcing the model with time-dependent values of σ from the Markov model yields more realistic short-term variations in terminus position than using a fixed, time-independent value (Figs 5, 10). By varying the calving stress threshold based on surface runoff, we achieve highly variable seasonal calving rates. This, in turn, yields seasonal advance and retreat cycles similar to those seen in observations (Fig. 10). While the Markov model controls the timing of seasonal fluctuations in terminus position, the maximum and minimum calving stress threshold parameters σmax and σmin must be constrained to produce realistic results. Most posterior samples have σmax values between ~250 and 1000 kPa (Fig. 9m). Probable values of σmin depend on σmax, but mostly fall between 100 and 500 kPa (Fig. 9o).

Fig. 9. Estimates of the posterior distribution obtained by running MCMC sampling using a Gaussian Process surrogate model fitted to an ensemble using dynamic or time-dependent parameters. Diagonal subplots a, f, j, m and o show marginal distributions for individual variables obtained using either the most probable surrogate model function (blue line) or accounting for surrogate model uncertainty (red line). Subplots below the diagonal show pairwise marginal distributions accounting for surrogate model uncertainty. Corresponding subplots above the diagonal show the same pairwise marginal distributions estimated using the Gaussian Process mean function (e.g. a.1 and a.2 represent the same marginal distributions). Red dots in subplots above the diagonal show 50 ensemble members sub-sampled from the prior based on their posterior probability and plotted in Figure 10.

Fig. 10. Observed terminus position at Helheim from 2003 to 2019 is indicated by the thick red line. Modeled terminus positions for all prior ensemble members are shown as faint black lines. Thin multi-colored lines show terminus positions for an approximation of the posterior distribution, which we obtain from the prior ensemble by weighting the prior ensemble members according to their probability and sub-sampling 50 ensemble members.

The most likely value of σmax is ~400 kPa, which is similar to the most probable time-independent value of σ found in Section 2.2. Consider that if we let σmax ≈ σmin, then we essentially identify the same highly probable static values of σ obtained before. However, while this results in a good fit to terminus position observations in the mean sense, this scenario does not yield realistic seasonal fluctuations. Hence, more interesting are posterior samples with σmax > σmin, which show greater temporal variability in terminus position, and we see that σmax values of ~500 kPa combined with σmin values of anywhere from ~100 to 400 kPa are probable as well.

Marginal distributions for the maximum and minimum seasonal undercutting rates, $\dot {m}_{\max }$ and $\dot {m}_{\min }$ respectively, show that undercutting rates likely do not exceed ~1 m d−1 at Helheim (Figs 9a, f). While terminus position is sensitive to ΔT, it is not well constrained by terminus position observations alone. Overall, the marginal distribution for ΔT favors higher temperature offsets above 0°, with a peak probability at ~7.5°C (Fig. 9j). Qualitatively, the marginal distribution for ΔT is similar to that obtained earlier for static parameters, with clusters at ~− 7.5, 0 and 7.5° (Figs 9j, 4f).

Ice temperature shows a number of interesting correlations with other parameters. For example, for a ΔT of ~7.5°C, posterior samples contain a wide range of values for $\dot {m}_{\max }$ and $\dot {m}_{\min }$ (Figs 4c.1, c.2). In contrast, for a ΔT of –7.5°C, posterior samples have a much narrower range of possible melt values. We also see a negative correlation between ice temperature and undercutting rates. This is likely because higher ice temperatures correspond to lower ice hardness and therefore lower calving rates (Eqn (3)). Hence, higher ice temperatures compensate for lower calving stress thresholds.

4.3. Comparison of static and dynamic posterior parameter distributions

Forcing the model with samples from the posterior distribution for time-dependent parameters yields realistic seasonal fluctuations that match observations in a qualitative sense (Fig. 10). In contrast, forcing the model with parameters drawn from the posterior for time-independent parameters (Section 2.2) reproduces the general trend of retreat at Helheim, yet does not appear to capture observed terminus position variabliity (Fig. 5). Here, we are interested in more formally assessing if the Markov model approach better captures the temporal variability of the observations.

To characterize the temporal variability of terminus position at Helheim, we detrend the observed terminus position by subtracting a linear function from the observation so that the initial and final terminus positions are both 0. We then compute the mean and variance of the detrendend observation, which yields a univariate normal distribution ${\cal N}( L_o,\; \, s^2_o)$ with a mean close to 0. The same procedure is repeated for 50 ensemble members from both the time-dependent and time-independent posterior distributions. That is, for each ensemble member, we detrend the modeled terminus position and compute a distribution ${\cal N}( L_i,\; \, s^2_i)$ that characterizes its temporal variability. We then compute the following metric z that compares the temporal variabililty between the model and observations for all ensemble members from both posterior distributions:

(23)$$z = {1\over 50} \sum_{i = 1}^{50} D_{KL} \left({\cal N}( L_i,\; s^2_i) ) \Vert {\cal N}( L_o,\; s^2_o) \right)$$

The Kullback–Leibler (KL) divergence D KL(P||Q) between two probability distributions P and Q can be thought of as a measure of similarity between those distributions. It has a value of 0 when P and Q are identical, and larger values indicate a greater discrepancy between probability distributions. Hence, a lower value of z reflects greater similarity between distributions characterizing the temporal variability for the modeled and observed terminus positions. For two univariate Gaussian distributions, the KL divergence has a simple closed form (Appendix D). Computing this metric for both posterior probability distributionsyields a value of 0.53 for the time-dependent posterior versus 0.87 for the time-independent posterior, indicating that the time-dependent posterior better captures the temporal variability in the observed terminus position.

As another test of the ability of the Markov model to capture seasonal terminus position variability, we examine seasonal rates of modeled advance/retreat for both time-dependent and time-independent parameters. For this purpose, we compare the average rate of advance/retreat for 50 ensemble members from both posterior distributions. For the time-independent parameters, we see a relatively steady average rate of retreat throughout the simulation with little temporal variability (Fig. 11). In contrast, for the time-dependent parameters, we see seasonal fluctuations in average modeled terminus position driven by variable calving rates, and periods of advance/retreat tend to align with the observations.

Fig. 11. The average rate of terminus position advance (positive) or retreat (negative) for 50 ensemble members drawn from the time-independent posterior (Section 2.2) is shown in red. Similarly, the average rate of terminus advance/retreat for 50 ensemble members drawn from the posterior distribution for time-dependent parameters (Section 4) is shown in blue. Shaded regions indicate periods of observed terminus advance.

5. Discussion

Bayesian inference of time-independent model parameters indicates that variations in terminus position at Helheim from 2007 to 2019 are driven predominantly by frontal ablation rather than ice flow or SMB. Terminus position is highly sensitive to the imposed calving stress threshold σ, which directly affects calving rates in a von Mises calving law (Fig. 4o). Modeled terminus position is somewhat sensitive to the undercutting rate $\dot {m}$ and ice temperature offset ΔT (Figs 4a, f). It is likely that the sensitivity of terminus position to temperature is primarily due to its effect on calving rates rather than ice flow, given that terminus position is relatively insensitive to variations in the basal traction field on the timescale of interest.

By constraining static parameter values, we achieve good fit between modeled and observed terminus positions in the mean sense. However, we are not able to reproduce short-term variations in observed terminus position (Fig. 5). Helheim undergoes marked seasonal cycles in terminus position, including advance during winter and retreat during summer, that are not adequately captured in the model (Bevan and others, Reference Bevan, Luckman and Murray2012; Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017). Exploring seasonally varying model forcings, we conclude that seasonal oscillations in subaqueous melt-driven undercutting as well as seasonal water pressure oscillations do not explain seasonal trends in terminus position. As expected, ice velocity is highly sensitive to seasonal water pressure oscillations, yet terminus position does not respond directly to ice acceleration or deceleration. Additionally, terminus position is sensitive to the undercutting rate but seasonal variations in subaqueous melt do not yield seasonal terminus position oscillations.

The sensitivity of modeled terminus position to the calving stress threshold σ and its relative insensitivity to other model parameters motivates the idea of dynamically adjusting calving rates by using a seasonally varying σ. Although this approach is somewhat ad hoc, there is physical justification for the idea of imposing seasonally variable calving rates.

For example, pooling of water in surface crevasses can help drive them downward, causing mechanical weakening of the ice that induces calving (Benn and others, Reference Benn, Warren and Mottram2007). Calving laws based on this physical mechanism have been found to reproduce terminus position observations better than other calving laws in some numerical experiments (Amaral and others, Reference Amaral, Bartholomaus and Enderlin2020). In this work, we adopted a von Mises calving law versus a surface crevasse calving law due to its numerical stability. Nonetheless, a surface crevasse calving model may provide a better physical basis for seasonal terminus position variations at Helheim, explaining why we might expect higher calving rates during the melt season due to pooling of surface runoff in crevasses.

Meltwater in surface crevasses is only one of several possible physical processes that explain highly variable calving rates on weekly to monthly time scales. Undercutting affects glacier calving by changing buoyant forces at the terminus, which can significantly amplify calving rates (Hanson and Hooke, Reference Hanson and Hooke2000; O'Leary and Christoffersen, Reference O'Leary and Christoffersen2013; Slater and others, Reference Slater, Benn, Cowton, Bassis and Todd2021). The pattern and magnitude of undercutting at the glacier front depend on subglacial runoff (Xu and others, Reference Xu, Rignot, Fenty, Menemenlis and Flexas2013; Rignot and others, Reference Rignot2016). Hence, enhanced subglacial runoff during the melt season could lead to greater undercutting and higher calving rates.

Other factors that could affect calving rates on a seasonal timescale include variations in mélange rigidity or basal crevassing. Backstress from the ice mélange may inhibit calving (Amundson and others, Reference Amundson2010; Walter and others, Reference Walter2012). Moreover, the extent and rigidity of the ice mélange vary seasonally, which could play a role in seasonally varying calving rates (Howat and others, Reference Howat, Box, Ahn, Herrington and McFadden2010). Finally, James and others (Reference James, Murray, Selmes, Scharrer and O'Leary2014) suggest that calving at Helheim could be tied to the development of basal crevasses due to buoyancy effects. Buoyancy at the terminus of Helheim in turn depends on subglacial water pressure on a seasonal timescale (Melton and others, Reference Melton2022).

Recently, experiments have been conducted using discrete element models of calving in three dimensions using a particle model (Benn and others, Reference Benn2017; van Dongen and others, Reference van Dongen2020). Notably, these models can resolve individual fractures and relate the stress state of the ice to its complex geometry at high resolution. Such models can reproduce a variety of hypothesized calving mechanisms such as enhanced calving due to undercutting as well as buoyancy-driven calving events (Benn and others, Reference Benn2017). Discrete element models support the idea that stress or strain fields in continuum models may provide a sufficient criteria for estimating calving rates (Benn and others, Reference Benn2017). However, this requires detailed knowledge of frontal geometry as well as the full, three-dimensional stress field (e.g. from a full Stokes stress balance). Hence, commonly used shelfy-stream models, with vertically integrated stress balances, may not resolve the necessary information needed to accurately predict calving rates (Benn and others, Reference Benn, Warren and Mottram2007). Moreover, the ability of models to resolve undercut regions is limited by model resolution.

Discrete models of calving are currently impractical for large-scale modeling studies due to high computational demand, and they are challenging to couple to standard continuum models of ice flow. However, they demonstrate the extreme state dependence of calving – from the fine scale geometry of undercut cavities down to the location of individual fractures. The intractability of representing the complex state of the glacier front means that we cannot expect to model calving in a completely deterministic sense, and indeed, it is the goal of most continuum calving models, including the one formulated in this work, to provide a reasonable approximation of calving rates given limited state information.

The multitude of interrelated processes potentially controlling calving and terminus position at Helheim motivates a stochastic approach to estimating calving rates. Rather than a single time-dependent calving rate, we generate a set of plausible “calving scenarios” using a Markov model. Since we cannot account for all of the physical processes governing calving dynamics, we aim to predict the timing and magnitude of calving events merely in a statistical sense. Many of the hypothesized calving mechanisms depend either directly or indirectly on surface runoff. For example, water depth in crevasses depends directly on surface runoff, while the geometry of undercut regions depends on subglacial discharge. Thus, we opt to relate calving rates to surface runoff, which is not only a proxy for seasonality but a potential driver of calving behavior. This approach is conceptually similar to the calving law in Sikonia (Reference Sikonia1982), which explored a similar relationship between river runoff, taken as a proxy for subglacial discharge, and calving rates at Columbia glacier in Alaska.

For Helheim, forcing the model with time-dependent calving stress threshold values using the Markov model yields realistic short-term variations in terminus position. Applying the model other Greenland marine-terminating glaciers in Greenland would be relatively straightforward given terminus position and runoff observations. Training the Markov model and randomly generating new realizations or ‘calving scenarios’ for use in ice dynamics models is computationally trivial given appropriately tuned values, or preferably probability distributions, for σmin and σmax. At Helheim, we find that σmax values in the vicinity of 500 kPa combined with σmin values of anywhere from ~100 to 400 kPa yield a good fit to observations (Fig. 10). For a fixed strain field, this implies up to a fivefold increase in modeled calving rates at Helheim from one 2-week period to the next using a von Mises calving law. It is likely that other glaciers would require different calving stress threshold values to achieve realistic results, and the main computational expense in this work is in estimating these parameters.

It is unclear if a similar relationship between terminus position and surface runoff exists for other marine-terminating glaciers, or if Helheim is a unique case. Potentially, the Markov model approach to modeling calving could be extended to incorporate additional observable inputs (e.g. ice velocity or the presence/absence of the melange). More sophisticated parametric models relating the calving stress threshold to observed inputs could also be used in place of the Markov model.

The link between calving rates and surface runoff at Helheim implies that it could be highly sensitive to future atmospheric warming. Enhanced surface melt caused by warmer future temperatures could lead to higher calving rates and accelerated retreat. Our results highlight the need to validate and tune calving laws using observations, as without proper tuning, we are not able to replicate seasonal terminus position variability at Helheim. Although we have not tested the stochastic calving model on other glaciers, this may provide an avenue for future work. Once trained, the stochastic calving model could also be used to explore retreat scenarios dependent on future surface runoff.

6. Conclusions

To investigate the physical processes controlling terminus position variability at Helheim between 2007 and 2019, we performed Bayesian inference of a collection of model parameters controlling ice flow, SMB and frontal ablation. Estimates of time-independent parameters show that model is sensitive to parameters controlling frontal ablation, particularly the calving stress threshold in the von Mises calving law, indicating the primary importance of calving dynamics at Helheim. Forcing the model with optimal time-independent parameters reproduces the mean observed terminus position at Helheim, yet fails to capture its short-term variability. An examination of observed seasonal calving behavior at Helheim, as well as a review of the physical processes driving calving led us to hypothesize a link between calving rates and surface runoff. Using a simple statistical model relating calving rates to observed surface runoff, we are able to reproduce Helheim's mean terminus position as well as its characteristic temporal variability.

Acknowledgements

J.D. and D.B. were supported by Heising-Simons Foundation Grant 2019-1157. M.M. was supported by Heising-Simons Foundation Grants 2019-1161 and 2021-3059. We thank Cheng Gong and Jesse Johnson for helpful discussions that aided in the development of this work.

Data availability

Code for the Markov model used in this work is available at https://github.com/JacobDowns/helheim_calving.

Appendix A.

Suppose that we have N potentially noisy evaluations of the full model ${\cal F}$

(A.1)$$y_i = {\cal F}( \pmb{x_i}) + \epsilon _i,\; $$

where i = 1,  2, · · · ,  N. Here, we assume that the noise εi is Gaussian distributed with mean 0 and variance σ2, denoted as

(A.2)$$\epsilon \sim {\cal N}( 0,\; \sigma^2) .$$

Let X be a matrix of training inputs given by

(A.3)$$X = [ \pmb{x_1},\; \pmb{x_2},\; \cdots,\; \pmb{x_N}] $$

and $\pmb {y}$ be a vector of corresponding full model outputs

(A.4)$$\pmb{y} = [ y_1,\; y_2,\; \cdots,\; y_N] ^T.$$

Given these training data, we would like to approximate the full model at new input locations $\pmb {x_i^\ast }$ for i = 1,  2, · · · ,  K. For convenience, we define a matrix containing these points

(A.5)$$X^\ast = [ \pmb{x_1^\ast },\; \pmb{x_2^\ast },\; \cdots,\; \pmb{x_K^\ast }] .$$

We adopt a Gaussian process surrogate ${\cal G}$, which takes in the matrix of inputs $X^\ast$ and returns a multivariate normal distribution of outputs

(A.6)$$\pmb{y^\ast } = \pmb{G}( X^\ast ) \sim {\cal N} \left(\pmb{m},\; P \right).$$

Since the GP surrogate returns a distribution rather than a simple vector of outputs, we are able to estimate uncertainty in surrogate model predictions. Training ${\cal G}$ requires defining a covariance kernel $k( \pmb {x},\; \, \pmb {x'}) \, \colon \, {\opf R}^n \times {\opf R}^n \to {\opf R}$ that encodes our assumptions about how outputs of ${\cal G}$ are correlated given the proximity of the inputs. We use a Matern Kernel of the form

(A.7)$$k( \pmb{x},\; \pmb{x'}) = {1\over \Gamma( \nu) 2^{\nu-1}} \left({\sqrt{2 \nu}\over l} {\rm d}( \pmb{x},\; \pmb{x'}) \right)^{\nu} K_{\nu} \left({\sqrt{2 \nu}\over l} {\rm d}( \pmb{x},\; \pmb{x'}) \right)$$

with ν = 1.5 and l a learned length scale hyper parameter. With this assumption we have

(A.8)$$\left[\matrix{\pmb{y} \cr \pmb{y^\ast } }\right]\sim {\cal N} \left(\pmb{0},\; \left[\matrix{K( X,\; X) + \sigma^2 I & K( X,\; X^\ast ) \cr K( X^\ast ,\; X) & K( X^\ast ,\; X^\ast ) \cr }\right]\right).$$

Here, the entries in the matrices K(X,  X), $K( X,\; \, X^\ast )$, and $K( X^\ast ,\; \, X^\ast )$ are given by

(A.9)$$[ K( X,\; X) ] _{i, j} = k( \pmb{x_i},\; \pmb{x_j}) $$
(A.10)$$[ K( X,\; X^\ast ) ] _{i, j} = k( \pmb{x_i},\; \pmb{x_j^\ast }) $$
(A.11)$$[ K( X^\ast ,\; X^\ast ) ] _{i, j} = k( \pmb{x_i^\ast },\; \pmb{x_j^\ast }) $$

and $K( X,\; \, X^\ast ) = K( X^\ast ,\; \, X) ^T$. Thus, using the properties of multivariate normal distributions, the mean and covariance of the outputs $\pmb {y^\ast }$ are given by

(A.12)$$\pmb{m} = K( X^\ast ,\; X) \left[K( X,\; X) + \sigma^2 I \right]^{-1} \pmb{y}$$

and

(A.13)$$P = K( X^\ast ,\; X^\ast ) - K( X^\ast ,\; X) \left[K( X,\; X) + \sigma^2 I \right]^{-1} K( X,\; X^\ast ) ,\; $$

respectively.

Appendix B.

Suppose that f(y) and g(y) are probability density functions for Gaussian distributions with means m 1 and m 2 and standard deviations σ 1 and σ 2 respectively. Then we have the following identity

(B.1)$$\int g( y) f( y) \; {\rm d}y = {-\frac{\sigma_1^2 \sigma_2^2 ( m_1 - m_2) ^2}{( \sigma_1^2 + \sigma_2^2) ^2} - \left(\frac{\sigma_1^2 \sigma_2^2}{\sigma_1^2 + \sigma_2^2} \right)\log \left(2 \pi ( \sigma_1^2 + \sigma_2^2) \right)\over \frac{2 \sigma_1^2 \sigma_2^2}{\sigma_1 \sigma_2}}.$$

To compute the posterior distribution, the variance of the surrogate prediction $P( y\vert {\cal G}( \pmb {x}) )$ is a function of $\pmb {x}$

(B.2)$$P( \pmb{x}) \propto \int_y P( \pmb{x}) P( y\vert {\cal G}( \pmb{x}) ) P( y) \; {\rm d}y$$
(B.3)$$ = P( \pmb{x}) \int_y P( y\vert {\cal G}( \pmb{x}) ) P( y) \; {\rm d}y.$$
(B.4)$$P( \pmb{x}) \propto F( \pmb{x}) .$$

Appendix C.

Suppose we have an observed sequence of states $s_1^n \equiv s_1,\; \, s_2,\; \, \cdots ,\; \, s_n$ generated by a Markov model with m states {1,  2, · · · ,  m} and an unknown transition matrix T with entries

(C.1)$$t_{i, j} = P( j\vert i) .$$

Using the Markov property, the probability of observing the sequence $s_1^n$ is

(C.2)$$P( s_1^n) = P( s_1) \prod_{2}^{n} P( s_t \vert s_{t-1}) .$$

Rewritten in terms of transition probabilities, we can write the likelihood that the sequence was generated by a given transition matrix

(C.3)$$L( T) = P( s_1) \prod_{2}^{n} t_{s_{t-1}, s_t}.$$

Defining N i,j to be the number of times state i is followed by state j in $s_1^n$, we have:

(C.4)$$L( T) = P( s_1) \prod_{i = 1}^{m} \prod_{\,j = 1}^{m} t_{i, j}^{N_{i, j}}.$$

Given the constraints

(C.5)$$\sum_j t_{i, j} = 1$$

it can be shown that L(T) is maximized when

(C.6)$$t_{i, j} = {N_{i, j}\over \sum_{\,j = 1}^m N_{i, j}}.$$

Appendix D.

Given two univariate normal distributions, ${\cal N}( \mu _1,\; \, \sigma _1^2)$ and ${\cal N}( \mu _2,\; \, \sigma _2^2)$, we have

(D.1)$$D_{KL} \left({\cal N}( \mu_1,\; \sigma_1^2) \Vert {\cal N}( \mu_2,\; \sigma_2^2) \right) = \log{{\sigma_2\over \sigma_1}} + {\sigma_1^2 + ( \mu_1-\mu_2) ^2\over 2 \sigma_2^2} - {1\over 2}.$$

References

Åkesson, H, Morlighem, M, O'Regan, M and Jakobsson, M (2021) Future projections of Petermann Glacier under ocean warming depend strongly on Friction Law. Journal of Geophysical Research: Earth Surface 126(6), e2020JF005921. doi: 10.1029/2020JF005921Google Scholar
Amaral, T, Bartholomaus, TC and Enderlin, EM (2020) Evaluation of iceberg calving models against observations from Greenland outlet glaciers. Journal of Geophysical Research: Earth Surface 125(6), e2019JF005444. doi: 10.1029/2019JF005444Google Scholar
Amundson, JM and 5 others (2010) Ice mélange dynamics and implications for terminus stability, Jakobshavn Isbræ, Greenland. Journal of Geophysical Research: Earth Surface 115(F1), F001405. doi: 10.1029/2009JF001405CrossRefGoogle Scholar
Benn, DI and 7 others (2017) Melt-under-cutting and buoyancy-driven calving from tidewater glaciers: new insights from discrete element and continuum model simulations. Journal of Glaciology 63(240), 691–702. doi: 10.1017/jog.2017.41CrossRefGoogle Scholar
Benn, D, Warren, C and Mottram, R (2007) Calving processes and the dynamics of calving glaciers. Earth-Science Reviews 82, 143–179. doi: 10.1016/j.earscirev.2007.02.002CrossRefGoogle Scholar
Bevan, SL, Luckman, A, Khan, SA and Murray, T (2015) Seasonal dynamic thinning at Helheim Glacier. Earth and Planetary Science Letters 415, 47–53. doi: 10.1016/j.epsl.2015.01.031CrossRefGoogle Scholar
Bevan, SL, Luckman, AJ and Murray, T (2012) Glacier dynamics over the last quarter of a century at Helheim, Kangerdlugssuaq and 14 other major Greenland outlet glaciers. The Cryosphere 6(5), 923–937. doi: 10.5194/tc-6-923-2012CrossRefGoogle Scholar
Bondzio, JH and 6 others (2016) Modelling calving front dynamics using a level-set method: application to Jakobshavn Isbræ, West Greenland. The Cryosphere 10(2), 497–510. doi: 10.5194/tc-10-497-2016CrossRefGoogle Scholar
Box, JE (2013) Greenland ice sheet mass balance reconstruction. Part II: surface mass balance (1840–2010). Journal of Climate 26(18), 6974–6989. doi: 10.1175/JCLI-D-12-00518.1CrossRefGoogle Scholar
Budd, WF, Keage, PL and Blundy, NA (1979) Empirical studies of ice sliding. Journal of Glaciology 23(89), 157–170. doi: 10.3189/S0022143000029804CrossRefGoogle Scholar
Cheng, D and 6 others (2021) Calving Front Machine (CALFIN): glacial termini dataset and automated deep learning extraction method for Greenland, 1972–2019. The Cryosphere 15(3), 1663–1675. doi: 10.5194/tc-15-1663-2021CrossRefGoogle Scholar
Cook, S and 7 others (2014) Modelling environmental influences on calving at Helheim Glacier in eastern Greenland. The Cryosphere 8(3), 827–841. doi: 10.5194/tc-8-827-2014CrossRefGoogle Scholar
Cook, SJ, Christoffersen, P and Todd, J (2021) A fully-coupled 3D model of a large Greenlandic outlet glacier with evolving subglacial hydrology, frontal plume melting and calving. Journal of Glaciology 68, 1–17. doi: 10.1017/jog.2021.109Google Scholar
Davison, BJ and 7 others (2020) Subglacial drainage evolution modulates seasonal ice flow variability of three tidewater glaciers in southwest Greenland. Journal of Geophysical Research: Earth Surface 125(9), 005492. doi: 10.1029/2019JF005492Google Scholar
Downs, J and Johnson, JV (2022) A rapidly retreating, marine-terminating glacier's modeled response to perturbations in basal traction. Journal of Glaciology 1–10. doi: 10.1017/jog.2022.5Google Scholar
Enderlin, EM and Howat, IM (2013) Submarine melt rate estimates for floating termini of Greenland outlet glaciers (2000–2010). Journal of Glaciology 59(213), 67–75. doi: 10.3189/2013JoG12J049CrossRefGoogle Scholar
Fried, MJ and 6 others (2019) Distinct frontal ablation processes drive heterogeneous submarine terminus morphology. Geophysical Research Letters 46(21), 12083–12091. doi: 10.1029/2019GL083980CrossRefGoogle Scholar
Gardner, A, Fahnestock, M and Scambos, T (2021a) ITS_live | Home Page.Google Scholar
Gardner, JR, Pleiss, G, Bindel, D, Weinberger, KQ and Wilson, AG (2021b) GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration. arXiv:1809.11165 [cs, stat], arXiv: 1809.11165.Google Scholar
Hanson, B and Hooke, R (2000) Glacier calving: a numerical model of forces in the calving-speed/water-depth relation. Journal of Glaciology 46, 188–196. doi: 10.3189/172756500781832792CrossRefGoogle Scholar
Holland, DM, Thomas, RH, de Young, B, Ribergaard, MH and Lyberth, B (2008) Acceleration of Jakobshavn Isbræ triggered by warm subsurface ocean waters. Nature Geoscience 1(10), 659–664. doi: 10.1038/ngeo316CrossRefGoogle Scholar
Howat, IM, Box, JE, Ahn, Y, Herrington, A and McFadden, EM (2010) Seasonal variability in the dynamics of marine-terminating outlet glaciers in Greenland. Journal of Glaciology 56(198), 601–613. doi: 10.3189/002214310793146232CrossRefGoogle Scholar
Iman, RL (2008) Latin hypercube sampling. In Encyclopedia of Quantitative Risk Analysis and Assessment, (Eds EL Melnick and BS Everitt), American Cancer Society, doi: 10.1002/9780470061596.risk0299.CrossRefGoogle Scholar
James, TD, Murray, T, Selmes, N, Scharrer, K and O'Leary, M (2014) Buoyant flexure and basal crevassing in dynamic mass loss at Helheim Glacier. Nature Geoscience 7(8), 593–596. doi: 10.1038/ngeo2204CrossRefGoogle Scholar
Kehrl, LM, Joughin, I, Shean, DE, Floricioiu, D and Krieger, L (2017) Seasonal and interannual variabilities in terminus position, glacier velocity, and surface elevation at Helheim and Kangerlussuaq Glaciers from 2008 to 2016. Journal of Geophysical Research: Earth Surface 122(9), 1635–1652. doi: 10.1002/2016JF004133CrossRefGoogle Scholar
Larour, E, Seroussi, H, Morlighem, M and Rignot, E (2012) Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM). Journal of Geophysical Research: Earth Surface 117(F1), F002140. doi: 10.1029/2011JF002140CrossRefGoogle Scholar
MacAyeal, DR (1989) Large-scale ice flow over a viscous basal sediment: theory and application to ice stream B, Antarctica. Journal of Geophysical Research: Solid Earth 94(B4), 4071–4087. doi: 10.1029/JB094iB04p04071CrossRefGoogle Scholar
Mankoff, KD and 10 others (2019) Greenland ice sheet solid ice discharge from 1986 through 2017. Earth System Science Data 11(2), 769–786. doi: 10.5194/essd-11-769-2019CrossRefGoogle Scholar
Mankoff, KD and 9 others (2020) Greenland liquid water discharge from 1958 through 2019. Earth System Science Data 12(4), 2811–2841. doi: 10.5194/essd-12-2811-2020CrossRefGoogle Scholar
Martos, YM (2018) Geothermal heat flux reveals the Iceland hotspot track underneath Greenland. Geophysical Research Letters 45(16), 8214–8222. doi: 10.1029/2018GL078289CrossRefGoogle Scholar
Melton, SM and 7 others (2022) Meltwater drainage and iceberg calving observed in high-spatiotemporal resolution at Helheim Glacier, Greenland. Journal of Glaciology 68(270), 812–828. doi: 10.1017/jog.2021.141CrossRefGoogle Scholar
Moon, T and 6 others (2014) Distinct patterns of seasonal Greenland glacier velocity. Geophysical Research Letters 41(20), 7209–7216. doi: 10.1002/2014GL061836CrossRefGoogle ScholarPubMed
Morlighem, M and 5 others (2010) Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica. Geophysical Research Letters 37(14), 043853. doi: 10.1029/2010GL043853CrossRefGoogle Scholar
Morlighem, M and 6 others (2016) Modeling of Store Gletscher's calving dynamics, West Greenland, in response to ocean thermal forcing. Geophysical Research Letters 43(6), 2659–2666. doi: 10.1002/2016GL067695CrossRefGoogle Scholar
Morlighem, M and 31 others (2017) BedMachine v3: complete bed topography and ocean bathymetry mapping of Greenland from multibeam echo sounding combined With mass conservation. Geophysical Research Letters 44(21), 11051–11061. doi: 10.1002/2017GL074954CrossRefGoogle ScholarPubMed
Motyka, RJ and 5 others (2011) Submarine melting of the 1985 Jakobshavn Isbræ floating tongue and the triggering of the current retreat. Journal of Geophysical Research: Earth Surface 116(F1), 116. doi: 10.1029/2009JF001632CrossRefGoogle Scholar
Mouginot, J and 8 others (2019) Forty-six years of Greenland ice sheet mass balance from 1972 to 2018. Proceedings of the National Academy of Sciences 116(19), 9239–9244. doi: 10.1073/pnas.1904242116CrossRefGoogle ScholarPubMed
O'Leary, M and Christoffersen, P (2013) Calving on tidewater glaciers amplified by submarine frontal melting. The Cryosphere 7(1), 119–128. doi: 10.5194/tc-7-119-2013CrossRefGoogle Scholar
Rasmussen, CE and Williams, CKI (2006) Gaussian processes for machine learning. Adaptive computation and machine learning, MIT Press, Cambridge, Mass, oCLC: ocm61285753.Google Scholar
Rignot, E and 12 others (2016) Modeling of ocean-induced ice melt rates of five west Greenland glaciers over the past two decades. Geophysical Research Letters 43(12), 6374–6382. doi: 10.1002/2016GL068784CrossRefGoogle Scholar
Rignot, E, Fenty, I, Menemenlis, D and Xu, Y (2012) Spreading of warm ocean waters around Greenland as a possible cause for glacier acceleration. Annals of Glaciology 53(60), 257–266. doi: 10.3189/2012AoG60A136CrossRefGoogle Scholar
Rignot, E, Koppes, M and Velicogna, I (2010) Rapid submarine melting of the calving faces of West Greenland glaciers. Nature Geoscience 3(3), 187–191. doi: 10.1038/ngeo765CrossRefGoogle Scholar
Rignot, E and Mouginot, J (2012) Ice flow in Greenland for the International Polar Year 2008–2009. Geophysical Research Letters 39(11), 051634. doi: 10.1029/2012GL051634Google Scholar
Sikonia, WG (1982) Finite element glacier dynamics model applied to Columbia Glacier, Alaska. Professional Paper, (doi: 10.3133/pp1258B.CrossRefGoogle Scholar
Slater, DA, Benn, DI, Cowton, TR, Bassis, JN and Todd, JA (2021) Calving multiplier effect controlled by melt undercut geometry. Journal of Geophysical Research: Earth Surface 126(7), e2021JF006191. doi: 10.1029/2021JF006191Google Scholar
Sundal, AV and 5 others (2011) Melt-induced speed-up of Greenland ice sheet offset by efficient subglacial drainage. Nature 469(7331), 521–524. doi: 10.1038/nature09740CrossRefGoogle ScholarPubMed
Tedesco, M and Fettweis, X (2020) Unprecedented atmospheric conditions (1948–2019) drive the 2019 exceptional melting season over the Greenland ice sheet. The Cryosphere 14(4), 1209–1223. doi: 10.5194/tc-14-1209-2020CrossRefGoogle Scholar
van der Veen, CJ (1998) Fracture mechanics approach to penetration of surface crevasses on glaciers. Cold Regions Science and Technology 27(1), 31–47. doi: 10.1016/S0165-232X(97)00022-0CrossRefGoogle Scholar
van Dongen, ECH and 5 others (2020) Numerical modeling shows increased fracturing due to melt-undercutting prior to major calving at Bowdoin Glacier. Frontiers in Earth Science, 8.CrossRefGoogle Scholar
Voytenko, D and 5 others (2015) Tidally driven ice speed variation at Helheim Glacier, Greenland, observed with terrestrial radar interferometry. Journal of Glaciology 61(226), 301–308. doi: 10.3189/2015JoG14J173CrossRefGoogle Scholar
Walter, JI and 6 others (2012) Oceanic mechanical forcing of a marine-terminating Greenland glacier. Annals of Glaciology 53(60), 181–192. doi: 10.3189/2012AoG60A083CrossRefGoogle Scholar
Williams, JJ, Gourmelen, N, Nienow, P, Bunce, C and Slater, D (2021) Helheim glacier poised for dramatic retreat. Geophysical Research Letters 48(23), e2021GL094546. doi: 10.1029/2021GL094546CrossRefGoogle Scholar
Wright, PJ, Harper, JT, Humphrey, NF and Meierbachtol, TW (2016) Measured basal water pressure variability of the western Greenland ice sheet: implications for hydraulic potential. Journal of Geophysical Research: Earth Surface 121(6), 1134–1147. doi: 10.1002/2016JF003819CrossRefGoogle Scholar
Xu, Y, Rignot, E, Fenty, I, Menemenlis, D and Flexas, MM (2013) Subaqueous melting of Store Glacier, west Greenland from three-dimensional, high-resolution numerical modeling and ocean observations. Geophysical Research Letters 40(17), 4648–4653. doi: 10.1002/grl.50825CrossRefGoogle Scholar
Figure 0

Fig. 1. Terminus position projected along a central flowline (a) and surface runoff (b) at Helheim glacier between 2003 and 2019. Shaded blue regions indicate periods of terminus advance.

Figure 1

Fig. 2. Idealized, conceptual example of summer SMB at a given point in time along a flowline for different values of $p_{\dot {a}}$. Values of − 1, 0 and 1 map to $\dot {a}_{\min }$, $\dot {a}_0$ and $\dot {a}_{\max }$.

Figure 2

Fig. 3. Location map, bedrock elevation and composite ice velocity map from 1980 to 2018 for Helheim glacier in East Greenland. The cyan line shows the central flowline along which terminus position is projected, while the cyan dot is used for point measurements in certain numerical experiments.

Figure 3

Fig. 4. Estimates of the posterior distribution obtained by running MCMC sampling using a Gaussian Process surrogate model fitted to an ensemble using static or time-independent parameters. Diagonal subplots a, f, j, m and o show marginal distributions for individual variables obtained using either the most probable surrogate model function (blue line) or accounting for surrogate model uncertainty (red line). Subplots below the diagonal show pairwise marginal distributions accounting for surrogate model uncertainty. Corresponding subplots above the diagonal show the same pairwise marginal distributions estimated using the Gaussian Process mean function (e.g. a.2 shows the same marginal distribution as a.1). Red dots in subplots above the diagonal show 50 ensemble members sub-sampled from the prior based on their posterior probability and plotted in Figure 5.

Figure 4

Fig. 5. Observed terminus position at Helheim from 2003 to 2019 is indicated by the thick red line. Modeled terminus positions for all prior ensemble members are shown as faint black lines. Thin multi-colored lines show terminus positions for an approximation of the posterior distribution, which we obtain from the prior ensemble by weighting the prior ensemble members according to their probability and sub-sampling 50 ensemble members.

Figure 5

Fig. 6. (a) Modeled terminus position at Helheim between 2007 and 2017 for experiments using seasonally varying undercutting rates. Line colors correspond to the mean undercutting rate over the duration of the simulation. (b) Percentage change in velocity at the cyan point marked in Figure 3a, relative to initial velocity. (c) Seasonal undercutting rates for a subset of 7 out of the 49 simulations displayed in panels (a) and (b).

Figure 6

Fig. 7. (a) Modeled terminus position at Helheim from 2007 to 2017 for runs using different magnitudes of seasonal water pressure oscillations. (b) Seasonal speed anomaly at the cyan point marked in Figure 3. The seasonal speed anomaly is computed by dividing the velocity at a given time by the velocity in a baseline run with no water pressure oscillations. (c) Seasonal water pressure oscillations expressed as a fraction of overburden.

Figure 7

Fig. 8. The Markov model outputs a calving stress threshold on a biweekly time step given calving stress threshold from the previous time step as well as surface runoff. Given the state at time i, denoted Si and predicts the state at the next time step Si+1 given the surface runoff ri.

Figure 8

Fig. 9. Estimates of the posterior distribution obtained by running MCMC sampling using a Gaussian Process surrogate model fitted to an ensemble using dynamic or time-dependent parameters. Diagonal subplots a, f, j, m and o show marginal distributions for individual variables obtained using either the most probable surrogate model function (blue line) or accounting for surrogate model uncertainty (red line). Subplots below the diagonal show pairwise marginal distributions accounting for surrogate model uncertainty. Corresponding subplots above the diagonal show the same pairwise marginal distributions estimated using the Gaussian Process mean function (e.g. a.1 and a.2 represent the same marginal distributions). Red dots in subplots above the diagonal show 50 ensemble members sub-sampled from the prior based on their posterior probability and plotted in Figure 10.

Figure 9

Fig. 10. Observed terminus position at Helheim from 2003 to 2019 is indicated by the thick red line. Modeled terminus positions for all prior ensemble members are shown as faint black lines. Thin multi-colored lines show terminus positions for an approximation of the posterior distribution, which we obtain from the prior ensemble by weighting the prior ensemble members according to their probability and sub-sampling 50 ensemble members.

Figure 10

Fig. 11. The average rate of terminus position advance (positive) or retreat (negative) for 50 ensemble members drawn from the time-independent posterior (Section 2.2) is shown in red. Similarly, the average rate of terminus advance/retreat for 50 ensemble members drawn from the posterior distribution for time-dependent parameters (Section 4) is shown in blue. Shaded regions indicate periods of observed terminus advance.