Inferring time-dependent calving dynamics at Helheim Glacier

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.

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, 2019). 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, 2012). 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, 2008; Motyka and others, 2011; Rignot and others, 2012). In addition to causing direct mass loss, subaqueous melting impacts ice dynamics by undercutting calving fronts, potentially amplifying calving rates (Fried and others, 2019; Slater and others, 2021).

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, 1998). Seasonal variations in meltwater input to the glacier bed cause changes in subglacial water pressure and basal sliding speed (e.g. Sundal and others, 2011). While seasonal acceleration is more pronounced on land-terminating outlet glaciers than marine-terminating glaciers (Davison and others, 2020), seasonal velocity variations may still have an important impact on calving (Cook and others, 2014). More importantly, surface runoff is linked to submarine melt rates via subglacial discharge (Xu and others, 2013). 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, 2021). 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, 2020; Åkesson and others, 2021), 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, 2019). 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, 2015; Kehrl and others, 2017). 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, 2014). Evidence based on geometry and nearby glaciers also suggests that Helheim may soon undergo substantial retreat (Williams and others, 2021).

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, 2012). 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, 2016) 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, 2012). Ice velocity is modeled using the shallow shelf approximation (MacAyeal, 1989). 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, 2021a). 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, 1979). 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, 2015). We therefore assume that water pressure is a fixed fraction (85%) of overburden pressure (Wright and others, 2016), 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, 2017). 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, 2016) to determine the calving rate at the ice-front position, which is tracked using a levelset method (Bondzio and others, 2016). 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, 2020). 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, 2010). To maximize spatial coverage of velocity observations, we again use the 1985–2018 composite ice velocity field from ITS_LIVE (Gardner and others, 2021a). The basal traction coefficient β ² is then given by

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

Uncertainty in β ² 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 (2022), 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 × 10³ ≤ σ ≤ 2000 × 10³ Pa. A lower threshold of 25 × 10³ Pa represents a physical extreme, where ice calves nearly immediately when entering the sea. The upper threshold of 2000 × 10³ 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, 2022), 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 (2013) estimated a mean submarine melt rate, $\dot {m}$, of ~0.56 m d⁻¹ at Helheim between 2000 and 2010. Rignot and others (2010) estimated melt rates of between 1 and 4 m d⁻¹ 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⁻¹. 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 ₀, 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 (2013), and a uniform geothermal heat flux of 70 mW m⁻² at the basal boundary based on Martos and 5 others (2018). 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, 2006) 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 (2021). 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, 2021) 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, 2008) 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⁻¹ and another at ~2 m d⁻¹, but melt rates below ~1 m d⁻¹ 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⁻¹ (Fig. 4b.1), and a negative correlation between $\dot {m}$ and the basal traction scaling parameter p _β for $\dot {m} > 1$ m d⁻¹ (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 (2013) 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, 2013). Surface runoff for the Helheim catchment area is estimated using Mankoff and others (2020).

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⁻¹. 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 (2014) 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, 2020). 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⁻¹, near terminus ice accelerates by ~300% over a matter of ~4 years (Fig. 6b). For mean melt rates below 2 m d⁻¹, velocity remains steady throughout the simulation, whereas melt rates above 2 md⁻¹ 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⁻¹ ($\dot {m}_{\max } = \dot {m}_{\min } = 3.5$ m d⁻¹) yields a retreat of 8 km over the period from 2007 to 2017, whereas an undercutting rate of 0 m d⁻¹ ($\dot {m}_{\max } = \dot {m}_{\min } = 0$ m d⁻¹) 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, 2020). 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, 2017). 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 (2014), 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 ₁,  B ₂, · · · ,  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 × 10³ ≤ σ_max ≤ 2000 × 10³ 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⁻¹ 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, 2012; Kehrl and others, 2017). 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, 2007). 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, 2020). 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, 2000; O'Leary and Christoffersen, 2013; Slater and others, 2021). The pattern and magnitude of undercutting at the glacier front depend on subglacial runoff (Xu and others, 2013; Rignot and others, 2016). 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, 2010; Walter and others, 2012). 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, 2010). Finally, James and others (2014) 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, 2022).

Recently, experiments have been conducted using discrete element models of calving in three dimensions using a particle model (Benn and others, 2017; van Dongen and others, 2020). 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, 2017). 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, 2017). 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, 2007). 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 (1982), 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.