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

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

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

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

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

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

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

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.

#### 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

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 × 10^{3} ≤ σ ≤ 2000 × 10^{3} Pa. A lower threshold of 25 × 10^{3} Pa represents a physical extreme, where ice calves nearly immediately when entering the sea. The upper threshold of 2000 × 10^{3} 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:

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

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

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

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}$.

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

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

### 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*

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*

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

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

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

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

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

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:

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.

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.

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

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

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.

### 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^{3} ≤ σ_{max} ≤ 2000 × 10^{3} 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).

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:

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.

## 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}$

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

Let *X* be a matrix of training inputs given by

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

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

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

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

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

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

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

and

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

To compute the posterior distribution, the variance of the surrogate prediction $P( y\vert {\cal G}( \pmb {x}) )$ is a function of $\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

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

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

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

Given the constraints

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

## Appendix D.

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