Calving dynamics at Jakobshavn Isbrae (Sermeq Kujalleq) controlled by local geometry: insights from a 3D Stokes calving model

Abstract

We present the first simulations of Jakobshavn Isbrae (Sermeq Kujalleq), west Greenland, using a 3D Stokes calving model that permits unrestricted advance and retreat. Using the position-based crevasse-depth calving law, the model is applied to simulate the calving dynamics of 2016–2017 season when Jakobshavn Isbrae is assumed to be stable because of the presence of a strong proglacial ice mélange. The calving law needs to be adjusted to avoid an underestimation of calving, but once adjusted the calving model simulates seasonal calving dynamics that reflect observed calving-driven retreat very well. We find that a crevasse penetration threshold of 94.5% best matches observations from satellite imagery. Additional, 2-year transient simulations show that although ice mélange is essential to the glacier's winter readvance, when removed, the glacier only retreats a couple of kilometres before reaching a stable position. While the backstress provided by the ice mélange allows the glacier to advance beyond this point, the retreated terminus position is determined by a combination of bed geometry and glacier dynamics. Ultimately, while the ice mélange allows winter readvance, cessation of the well-documented rapid retreat of Jakobshavn Isbrae will be influenced by the bed geometry.

1. Introduction

Mass loss from the Greenland Ice Sheet (GrIS) has been increasing since the 1980s to 286 ± 20 Gt a⁻¹ between 2010–2018 with tidewater glacier calving accounting for over 40% of the total (Mouginot and others, 2019; Mankoff and others, 2020). Large tidewater glaciers dominate these contributions. For example, Jakobshavn Isbrae (Sermeq Kujalleq), Greenland's largest outlet glacier, accounts for approximately a fifth of the total mass loss from the GrIS (Enderlin and others, 2014). Understanding the complex calving dynamics at Jakobshavn Isbrae, west Greenland, alone will improve estimates of future global sea level rise contributions from the GrIS.

The complexity of the tidewater glacier system makes developing detailed models of the evolution of Greenlandic glaciers challenging. Firstly, the wide range of interacting controls, including submarine melting (O'Leary and Christoffersen, 2013), buttressing provided by a proglacial ice mélange (Cassotto and others, 2021), and bedrock and fjord geometry (Enderlin and others, 2013; Todd and others, 2018; Benn and others, 2023) makes the separation of the influence of individual processes difficult. Secondly, including these processes within models is an additional challenge.

Large-scale studies have shown a sensitivity between ice sheet and ocean temperatures (Holland and others, 2008; Christoffersen and others, 2011; Cowton and others, 2018). The ocean heat is often assumed to be transferred to glaciers through background fjord conditions and circulation (How and others, 2019) or by turbulent meltwater plumes (Slater and others, 2017). Plumes have been hypothesised to directly influence calving through local melting (Cowton and others, 2019; Wagner and others, 2019). However, at fast flowing glaciers, such as Jakobshavn Isbrae, modelled frontal melt rates are far less than glacier speeds, so the impact of frontal melting is much less than for the case of slowly moving glaciers (Joughin and others, 2020; Kajanto and others, 2023). Instead the sensitivity to background ocean temperature is hypothesised to be conveyed through the backstress from mixtures of sea ice and calved icebergs known as ice mélange (Joughin and others, 2020). Background water temperatures control the presence of a proglacial ice mélange which has been associated with the stabilisation of tidewater glaciers. The buttressing provided by the ice mélange can suppress calving by opposing longitudinal stretching and iceberg overturn, but the granular nature of the ice mélange makes this hard to model or quantify (e.g. Amundson and others, 2010; Todd and others, 2018; Cassotto and others, 2021). The most referenced ice-mélange backstress quantification is 30–60 kPa through the full height of Store Glacier (also known as Sermeq Kujalleq in Greenlandic), from the field observations of Walter and others (2012). This equates to 235–470  kPa resistive pressure over the contact area between the glacier and the ice mélange (Todd and others, 2018, 2019). When used in modelling studies, the complexity of ice-mélange dynamics (Cassotto and others, 2021) is typically reduced to a binary boundary pressure that is turned off or on between simulations or over the course of a simulation (e.g. Todd and others, 2018; Cook and others, 2022).

Shallow or narrow reaches of a fjord tend to suppress calving, and thus act as pinning points within the system such that glacier termini tend to stabilise at such locations. If subjected to a sufficiently large perturbation, such as an increase in submarine melting or a reduction in buttressing from ice mélange, the glacier may rapidly move from one stable position to another (Bevan and others, 2019; Benn and others, 2023). To be effective, therefore, calving models need to exhibit appropriate sensitivity to submarine melting and the presence or absence of ice mélange, as well as the tendency for ice fronts to stabilise at pinning points when conditions permit.

There is a large spectrum of methods to model calving, from the most computationally expensive first principle particle models (e.g. Aström and others, 2014) or damage mechanics (e.g. Mercenier and others, 2019) to simpler ice-sheet scale modelling. Here, we focus on the simplified calving functions that can easily be transferred to larger continuum ice-sheet scale modelling. Modelling calving in a continuum model requires a so-called ‘calving law’ to predict either the calving rate or the position of the calving front at each timestep. Calving laws make predictions based on the geometry, stress state or velocity field of the glacier, but most require tuning on a case-by-case basis (e.g. Choi and others, 2018; Amaral and others, 2020). In order to replicate the behaviour of glacier terminus movement between stable states, a position-based calving law must be used. In this study, we use the crevasse-depth calving law (CDL), which predicts the position of the calving margin from the estimated penetration of basal and surface crevasses based on the stress field (Benn and others, 2007a; Nick and others, 2010). Crevasses are not physically represented in these simulations; rather, the calving prediction is based on the extent of regions within the ice which are assumed to be crevassed to greater or lesser extent on the basis of the stress state. In its original form, CDL predicts calving when either surface crevasses penetrate to the waterline or surface and basal crevasses together penetrate the full thickness of the glacier. Several studies have found that surface crevasses must be assumed to contain water to allow the CDL to match observed terminus positions (Cook and others, 2014; Lea and others, 2014; Enderlin and Bartholomaus, 2020), and it has been rightly noted that the required water depths are not physically realistic (Amaral and others, 2020). This problem can be overcome if the CDL is modified to predict calving at some lesser amount of crevasse penetration. Previous Stokes modelling studies at Store Glacier argued that this step is not necessary and that the model performs well if the original criteria are used (e.g. Todd and others, 2018, 2019; Benn and others, 2023). However, the CDL underestimated calving at Store Glacier (Cook and others, 2022, 2023) suggesting that a slightly adjusted calving criterion may have improved results. Calving models must also overcome significant technical hurdles, such as remeshing, in order to implement the calving law through a calving algorithm. The difficulty of calving implementation is less often discussed than calving laws themselves, but consideration of the calving algorithm is essential to addressing model uncertainties. Until now, the intricacy of implementing calving in 3D has limited modelling studies to stable glaciers such as Store Glacier (Todd and others, 2018, 2019; Cook and others, 2022). The study presented here takes advantage of a new calving algorithm that makes it possible to simulate the 3D calving dynamics for a much greater range of glacier circumstances (Wheel and others, 2024). The algorithm's capabilities that allow this study are:

1. 1. unlimited advance or retreat in 3D,

2. 2. unrestricted 3D calving geometries can be utilised by the model, and

3. 3. any vertically or horizontally varying melt field can be applied to the glacier front.

Together, the CDL and new calving algorithm allow the influence of submarine melt, ice mélange and fjord geometry to be modelled explicitly within a single framework. In this paper, we apply the calving algorithm to seasonal fluctuations of Jakobshavn Isbrae in the first study to use a 3D Stokes model to investigate the controls on the calving dynamics of a historically unstable glacier. A first set of simulations were run to optimise the crevasse penetration threshold in the calving law by finding the best match to observed seasonal patterns of calving front positions. Further simulations were conducted to investigate the possible role of ice mélange in controlling the position of the ice front. The results show that although ice mélange allows the ice front to advance in winter, its absence would not allow unlimited retreat of the ice in subsequent summers. Rather, the present summer minima of Jakobshavn Isbrae appear to be controlled by a topographic pinning point. The close match between model results and observations indicates that the CDL implemented in a Stokes model provides a powerful tool for predicting the future evolution of tidewater glaciers in Greenland.

2. Data and methods

The simulations presented here were conducted with the open-source, 3D Stokes finite-element model Elmer/Ice (Gagliardini and others, 2013). Elmer/Ice can solve a range of problems including the grounding dynamics (Favier and others, 2012) and ice temperature evolution (Gagliardini and others, 2013) and has a range of methods to tackle inverse problems (Gillet-Chaulet and others, 2012). We implement the calving algorithm developed by Wheel and others (2024), which builds on the work of Todd and others (2018). In this section we detail the modelling procedure, boundary conditions and input data particular to the simulations of Jakobshavn Isbrae. More information on the calving algorithm and examples of the model input files can be found in Wheel and others (2024) and from the official source code (https://github.com/ElmerCSC/elmerfem).

2.1 Study site

Jakobshavn Isbrae is a large, well-studied tidewater glacier in West Greenland (e.g. Iken and others, 1993; An and others, 2017; Khazendar and others, 2019; Xie and others, 2019; Cassotto and others, 2021). It underwent a rapid and widely documented retreat since the turn of the millennium (Joughin and others, 2012) although since 2013 the retreat came to an unexpected halt and annual minima have maintained similar positions since that time (Joughin and others, 2020). It has been argued that the retreat phase was triggered by warmer ocean conditions (Holland and others, 2008; Rignot and others, 2015), and was arrested by buttressing from proglacial ice mélange, emphasising the role of ocean drivers in the recent behaviour of the glacier (Joughin and others, 2020). The technical advances made to the Elmer/Ice calving model as described in Wheel and others (2024) give an unparalleled opportunity to investigate such influences on calving dynamics at Jakobshavn Isbrae. Model results are compared to Sentinel-2 or Landsat 8 and 9 images or from the synthetic aperture radar (SAR) records provided by Joughin and others (2020).

2.2 Model setup

We use a 3D model domain of Jakobshavn Isbrae that extends 20 km inland from the terminus position on 2 October 2016 (Fig. 1). This start date was chosen because the summer minimum position of the glacier has remained relatively stable since that time despite large seasonal fluctuations (Joughin and others, 2020). During late 2016 and 2017, the glacier underwent a substantial annual cycle of advance and retreat of several kilometres, providing very challenging calving dynamics to replicate in a model. Interested readers are pointed to Wheel and others (2024) for a full discussion of these challenges and how the new Elmer/Ice calving algorithm addresses them. Fjord boundaries beyond the terminus were defined to allow the model to advance beyond its starting domain (Fig. 1). The starting mesh was composed of linear tetrahedral elements, producing a mesh that was fully unstructured both vertically and horizontally. The vertical resolution was roughly 50 m across the domain while horizontal resolution varies from 100 m at the terminus to 1500 m further inland. The upper surface was defined using a digital elevation model from Joughin and others (2020) with the ArcticDEM used for the remaining upstream area. The basal surface followed the bed topography provided by BedMachine (Morlighem and others, 2017). The domain had six boundary conditions that consist of the calving front of the glacier (Γ_front), both lateral margins (Γ_left, Γ_right), the interior inflow (Γ_inflow), the base (Γ_base) and the top surface (Γ_surf). The conditions applied to each boundary are discussed below.

Figure 1. Model mesh of Jakobshavn Isbrae. (a) Top view of the surface elements of the 3D mesh in black with a Sentinel-2 image in the background. The mesh density increases towards the terminus. The red lines are the fjord walls that are used to project lateral advance. (b) Grounding zone of the model after the spin up shown as a grounded mask viewed from the top. Red areas (values of 1) indicate grounded ice, blue areas (values of −1) are ungrounded ice and values of 0 are the grounding lines.

2.2.1 Ice dynamics

Glen's flow law is applied to specify the glacier rheology, which is then used to solve the Stokes equations (Cuffey and Paterson, 2010) – often referred to in literature as ‘full-Stokes’ modelling. The temperature dependency of Glen's flow law is calculated using the Arrhenius equation following Gagliardini and others (2013) and Eqns (3) and (4) therein. The rate factors are detailed in Cuffey and Paterson (2010). The model accounts for temperature increases due to internal dissipation and friction. The ice temperature was calculated during the spin-up of the model and coupled into the ice dynamics. The short time scale of the transient experiments undertaken means it is not necessary to continually solve for the temperature. Instead, the temperature field is assumed to remain unchanged from the spin-up in forward simulations. This is similar to the method detailed in Todd and others (2018). Borehole data from Iken and others (1993) were used to provide a temperature profile for the inflow boundary. The surface values from Iken and others (1993) are used as the surface boundary condition. Geothermal heating of 60 mW m⁻² is specified on the base boundary (Martos and others, 2018).

The lateral boundaries adhere to the fjord walls such that a ‘no penetration’ condition is applied. A simple linear friction law was applied at the lateral boundaries using a constant slip coefficient of 1e-3 MPa m⁻¹ a⁻¹, while a linear friction law is applied on the glacier base with the slip coefficient here calculated using the Adjoint method (Gillet-Chaulet and others, 2012) based on MEaSUREs monthly velocity maps (Joughin and others, 2021). The initial estimate of basal slip was derived using the shallow ice approximation (SIA) and the topography of the starting domain. Since Elmer/Ice solves the grounding dynamics, this condition is only applied where the glacier is grounded. The lateral and grounded ice boundary conditions are then

(1)$$u_{\perp} = 0,\; \rm{ on } \;\Gamma_{{\rm left}},\; \Gamma_{{\rm right}},\; \Gamma_{{\rm base}},\; $$

(2)$$\sigma_{\Vert } = -u_{\Vert }\beta,\; \quad \rm{on } \;\Gamma_{{\rm left}},\; \Gamma_{{\rm right}},\; \Gamma_{{\rm base}},\; $$

where u is the velocity component, σ the stress component, β is the slip coefficient and the perpendicular and tangential components are shown by ⊥ and || respectively. For basal areas where the ice becomes ungrounded, no friction is applied and the glacier is free to move vertically. Akin to ungrounded ice, the terminus face experiences no friction, but an external water pressure is applied as a normal stress:

(3)$$\sigma_{\perp} = \hbox{min}( -\rho_{\rm w} gh,\; 0) ,\; \quad \rm{on } \;\Gamma_{{\rm base}},\; \Gamma_{{\rm term}},\; $$

(4)$$\sigma_{\Vert } = 0,\; \quad \rm{on } \;\Gamma_{{\rm base}},\; \Gamma_{{\rm term}},\; $$

where ρ _w is the density of the water, g is the gravitational acceleration and h is the depth below the water level. A stress-free condition is applied to the surface boundary, while the velocity of the inflow boundary is set to match velocity observations from the MEaSUREs dataset (Joughin and others, 2021).

2.2.2 Free surfaces

The free surfaces at both the base and surface of the glacier are solved at each timestep. The surface evolves in response to a surface mass-balance boundary condition and ice dynamics. The surface mass balance is provided by daily RACMO 2.3p2 data which have a horizontal resolution of 1 km and is linearly interpolated onto the surface of the glacier mesh (Noël and others, 2018). The majority of the basal boundary is grounded and as such has a no penetration condition. Ungrounded ice is free to evolve in response to ice dynamics and no additional melt is provided to the base boundary.

2.2.3 Submarine melting and ice mélange

A background melt field and a meltwater plume are applied as a 2D melt field to the calving face. The melt rates are based on the detailed MITgcm simulations undertaken by Kajanto and others (2023), and the meltwater plume is 2 km wide and applied to the central portion of the calving face. The plume is not coupled into the subglacial hydrology like Cook and others (2022), rather a predefined melt field varies temporally based on the predefined MITgcm output.

A backstress acting on the terminus from proglacial ice mélange has often been included in calving models (e.g. Todd and others, 2018, 2019; Cook and others, 2022). However, estimations of the exact normal stress exerted by the ice mélange vary greatly. The values of backstress that are most commonly used are from Walter and others (2012), who used velocity observations at Store Glacier to estimate a mélange backstress of 30–60 kPa averaged for the full height of the glacier calving face. When used in 3D modelling studies at Store Glacier, this was converted to a 240–480 kPa pressure over the contact area (Todd and others, 2018) which stretched across the full terminus width and had a vertical thickness of 75 m. Here, we assume that the ice mélange at Jakobshavn Isbrae has a greater contact height of 100 m. Using the Walter and others (2012) estimation, this equates to a total contact backstress of 144–288 kPa. The ice-mélange pressure is applied between 1 November and 1 May each year based on the observations of Cassotto and others (2015).

2.2.4 Calving criteria and modelled calving

The calving model is split into the calving law that predicts calving and the calving algorithm that implements the calving prediction. In its original form, the CDL predicts calving under two conditions: (1) the surface crevasse field extends to the waterline (Benn and others, 2007a, 2007b) and (2) basal and surface crevasse fields extend through the full thickness of the glacier (Nick and others, 2010). Using a modified Nye (1957) criterion using the principal stress (Otero and others, 2010), the vertical extents of crevasse fields are predicted by

(5)$$\sigma_{\,p, {\rm surf}} = \sigma_1,\; $$

(6)$$\sigma_{\,p, {\rm basal}} = \sigma_1 + P_{\rm w},\; $$

where σ ₁ is the largest Cauchy stress and P _w is the water pressure at basal crevasses. Since the full stress tensor is used the first principal stress varies vertically and locally. Water pressure in basal crevasses is calculated following Todd and others (2018). Vertical ray casting is used to map surface and basal crevasse penetration onto a 2D planar mesh (Todd and others, 2018). This produces three key outputs:

(7)$$H_{{\rm crev}, {\rm surf}} = {d_{{\rm surf}}\over z},\; $$

(8)$$H_{{\rm crev}, {\rm basal}} = {d_{{\rm surf}} + d_{{\rm basal}}\over H},\; $$

(9)$$H_{{\rm crev}} = \hbox{max}( H_{{\rm crev}, {\rm surf}},\; H_{{\rm crev}, {\rm basal}}) ,\; $$

where the depths of basal and surface crevasses are d _surf and d _basal respectively, z is the ice freeboard and H is ice thickness. Additionally, the ‘Calving index’ is 1 − H _crev and represents the relative proportion of intact ice. Although the vertical ray casting assumes the terminus immediately exposed by calving is vertical, the terminus is free to deform in any direction based on the velocity profile and submarine melting.

A new term, the crevasse penetration threshold, is introduced to control the sensitivity of the calving law and has a real value between 0 and 100% (or a value of 0 to 1). The crevasse penetration threshold is intended to account for processes not included in the model (e.g. stress concentrations or inherited ‘damage’), without invoking physically unrealistic water depths in surface crevasses. For example, assigning a value of 75% to the crevasse penetration threshold will mean that calving will be induced when either: (1) a surface crevasse reaches 75% of the distance to the water line (H _crev,surf) or (2) surface and basal crevasses extend through 75% of the entire ice column (H _crev,basal). The CDL is unchanged from its original form when assigned a value of 100%, while values below this increase the sensitivity of the calving law.

Contours delineating where the crevasse penetration threshold is met are isolated and undergo postprocessing to check whether an isolated volume of ice (‘iceberg’) is capable of being evacuated from the terminus (Wheel and others, 2024). Any potential bottlenecks are removed. A signed distance variable is then calculated on the 3D glacier mesh where negative values indicate areas isolated by crevassed ice and positive values indicate the main glacier body. The signed distance variable is then used to create a new domain through the remeshing presented in Wheel and others (2024). The key methodological advances in the postprocessing of crevasses that allow the simulation of calving at Jakobshavn Isbrae are:

1. 1. the calving front can potentially take any form, as there is no requirement for projectability in the calving law, and

2. 2. the presence of the fjord walls is taken into account when the final calving prediction is made in the algorithm.

These advances are detailed and visualised in Wheel and others (2024). The input parameters for the calving algorithm were set to maintain the mesh density of 100 m at the terminus before increasing to 500 m upstream. The area within 1500 m of the predicted new terminus position was remeshed at each timestep. Consequently, the area being remeshed moved as the terminus advanced or retreated. Calving was implemented at each timestep, with up to three extra calving cycles following large calving events (adaptive time stepping: Wheel and others, 2024), to allow stabilisation of the calving front before calculating a new stress and velocity solution.

2.3 Model procedure

2.3.1 Model spin up

A model spin-up was undertaken to allow model components to relax before the transient simulations were begun. The workflow was as follows:

1. 1. an initial steady-state inversion based on an initial estimate of basal slip,

2. 2. a steady-state simulation to produce converged temperature and velocity fields,

3. 3. a secondary steady-state inversion using the temperature field from (2) to produce a basal slip map, and

4. 4. a 5-year transient simulation with a fixed terminus geometry to allow the surface and bed to relax.

Steps 1–3 were replicated for each month between 2016–10 and 2017–09 to produce monthly maps of basal slip conditions. For the inversions, monthly velocity maps were used from the MEaSUREs dataset as the observed velocity for the surface boundary (Joughin and others, 2021). The relatively short model relaxation in (4) is justified with the assumption that the geometry of Jakobshavn Isbrae over October 2016 was not stable. The timestep of the relaxation varied from 0.001a to 0.005a. Calving is a process resulting from an unstable geometry and acts to move the terminus to a more stable position. Therefore, there is potential that a prolonged relaxation may produce an unrealistically stable geometry. The model spin-up produced a complex grounding zone that shows similar patterns to other modelling studies at Jakobshavn Isbrae where the majority of ungrounded ice is on the north side of the glacier (Bondzio and others, 2016, Fig. 1).

2.3.2 Refinement of calving law

The final relaxation from the spin-up procedure described above was used as the starting point for all transient experiments. Initial simulations with the 100% crevasse penetration threshold caused the CDL to underestimate observed calving, which guided further simulations in which the crevasse penetration threshold was adjusted to provide a better match with observed terminus changes. To determine the optimal crevasse penetration threshold, 1-year transient simulations with 1 d timesteps were run using the starting model domain described above. The calving model was active during these simulations, allowing for terminus advance and retreat. The crevasse penetration threshold was adjusted between each simulation based on over- or under-prediction of the observed change in the previous simulation. In addition to the adjustment of the CDL, the magnitude of ice-mélange back pressure was varied to assess its impact on the winter readvance of the glacier. Using Walter and others (2012), values of 144 and 288 kPa of ice-mélange pressure were applied for the upper-most 100 m of the submerged terminus when calculating the crevasse penetration threshold.

2.3.3 Two-year experiments

To investigate calving dynamics at Jakobshavn Isbrae, several 2-year transient experiments were conducted once optimal crevasse penetration thresholds and ice-mélange back pressure had been identified. The input parameters such as basal conditions, surface mass balance, melt and inflow velocity remained unchanged annually. In effect, this means that the conditions between October 2016 and October 2017 were repeated for a subsequent 365 d period. The only change between the simulations was the presence or absence of an ice-mélange pressure to the calving face over the winter periods as defined by Cassotto and others (2015). The experiment variations are outlined in Table 1. For experiments 3 and 4 the geometry was taken from the first year of experiments 1 and 2, respectively.

Table 1. Summary of ice-mélange backstress application over 2-year transient experiments

3. Results

3.1 Refinement of the calving law

Using a crevasse penetration threshold of 100% with an ice-mélange backstress of 30kPa (144kPa contact pressure) during the winter period, the glacier terminus continually advanced throughout winter and summer (Fig. 2a). Although the whole calving face advanced, a large central peninsula formed and advanced throughout the simulation, undergoing little calving activity. The final terminus position was far more advanced than the observed terminus position in October 2017, and the terminus shape was a poor match to observations.

Figure 2. Terminus positions for 1-year model runs from 2016 to 2017 with varying crevasse penetration thresholds. (a) Crevasse penetration threshold = 100% and ice mélange pressure of 30 kPa. (b) Crevasse penetration threshold = 95% with a 30 kPa ice-mélange pressure. (c) Crevasse penetration threshold = 95% with a 60 kPa ice-mélange pressure. (d) Crevasse penetration threshold =  93.5% with a 60 kPa ice-mélange pressure. The full simulations can be seen in S2.1 to S2.4.

In order to obtain a better match with observed terminus positions, the crevasse penetration threshold was progressively reduced. Calving substantially increased when the crevasse penetration threshold was assigned a value of 95% (Fig. 2b). The glacier terminus advanced during winter conditions when the ice mélange was present before retreating to roughly its starting position over the summer months. The terminus followed a more realistic concave profile, indicative of a compressive arch.

If the ice-mélange backstress was doubled to 60 kPa (contact pressure of 288 kPa), the high-end of the Walter and others (2012) estimates, the central portion of the glacier underwent a large advance over the winter months (Fig. 2c). This pattern of terminus advance was not seen in the satellite observations. Once the ice-mélange pressure was removed, the glacier calved back to a similar position as the simulation using a lower estimate of ice-mélange backstress. To determine whether the unrealistic winter readvance shown in Fig. 2c was a consequence of poor ice-mélange back pressure estimations or the crevasse penetration threshold assigned in the calving law, a further simulation was conducted with a crevasse penetration threshold of 93.5% (Fig. 2d). This simulation had an ice-mélange back pressure of 60 kPa and again formed a large peninsula during the winter months that calved off at the start of summer. The retreat in summer was much larger than that observed at Jakobshavn Isbrae. The most significant calving occurred on the northern margin of the glacier where the majority of ungrounded ice was present (Figs 1b and 2d).

The simulation that matched best with observations was found using a crevasse penetration threshold of 94.5% and an ice-mélange backstress of 30 kPa (Fig. 3). This simulation tracked the annual terminus pattern at Jakobshavn Isbrae (Fig. 3a). The majority of the ice lost from calving was from full thickness calving events with limited ice lost from smaller serac calving events. Seasonal velocity patterns were similar between the model and the observations but the observed velocity was up to 2 km a year quicker. Changes in the model velocity conincided with changes in either calving, the presence of the winter ice-mélange or basal conditions (Fig. 3a). The terminus advanced during the winter when the ice mélange was present, reaching a seasonal maxima in May (Fig. 3). At the end of the summer, the terminus retreated to within 100 m of the width averaged observed glacier position (Fig. 3b). Whereas the terminus positions agree well both at the end of winter period and end of summer, the major discrepancy between the modelled terminus and observations were at the start of summer. Observations show the terminus continued to advance until mid-June when the glacier began its summer retreat (Joughin and others, 2020). The simulation instead showed retreat as soon as the ice-mélange buttressing was removed on 1 May.

Figure 3. Simulation for 2016–2017 using a crevasse penetration threshold of 94.5% and an ice-mélange pressure of 30 kPa. (a) Modelled (red lines) and observed (blacked dots) terminus position, surface velocity and surface elevation. Observations are from Joughin and others (2020). Both observations and modelled variables for the location denoted by the black dot in panel (b) which sits on the flowline from Joughin and others (2020). The inversion parameter beta is also shown as proxy for seasonal basal conditions. (b) Terminus positions from the simulation above a Sentinel-2 satellite image of the calving front on 2017-09-27 at the end point of the modelled timespan. The full simulation can be seen in S3.

3.2 Two-year simulations

The first year of the transient 2-year simulation with ice mélange present (experiment 1, Table 1) used the same conditions described when the calving law had been optimised (Fig. 4a). The second year of this simulation showed the same pattern with a winter readvance when an ice-mélange pressure was applied, before retreating to roughly the start position (Fig. 4c). Conversely, the glacier underwent a sustained retreat throughout the year when the ice-mélange was not applied in the first year (Fig. 4b). There was no seasonality to the terminus evolution with a uniform retreat through the winter and start of summer. Retreat was concentrated on the north side of the glacier where retreat was greater than 2 km. The area of retreat again coincided with the location of the majority of the ungrounded ice. The glacier's south margin remained stable – retreating and advancing a few hundred metres around one point. Retreat slowed towards the end of the year and the glacier reached a stable position. The stable position was maintained (Fig. 4d) over the second year when still no ice mélange was present (experiment 2, Table 1). Very little movement of the terminus occurred over the year. When the terminus did advance beyond this position, calving eroded it back to the stable point but never caused it to retreat farther.

Figure 4. Terminus outputs from 2-year experimental model runs using a 94.5% crevasse penetration threshold and 30 kPa ice-mélange pressure. The background Sentinel-2 satellite image is from 2017-09-27. For each simulation the experiment number and year as noted in Table 1 is in brackets. (a) Terminus positions during the first year when an ice mélange was applied (experiment 1, year 1; experiment 3, year 1). (b) Terminus positions during the first year when an ice mélange was not applied (experiment 2, year 1; experiment 4, year 1). (c) Terminus positions during the second year when an ice mélange was applied for both years (experiment 1, year 2). (d) Terminus positions during the second year when an ice mélange was not applied for both years (experiment 2, year 2). (e) Terminus positions during the second year when an ice mélange was applied for first year but not the second (experiment 3, year 2). (f) Terminus positions during the second year when an ice mélange was applied for second year but not the first (experiment 4, year 2). The star indicates the year presented. (g) Width averaged terminus position over the course of each experiment. The full simulations can be seen in S4.1 to S4.6.

When the geometry from the end of experiment 1, year 1 was run for another year, but this time without ice mélange (experiment 3, Table 1), the glacier rapidly retreated (Fig. 4e). After the initial rapid retreat over the first 4 months, it reached the same stable position observed in the no ice-mélange simulation outlined above. The terminus remained stable for the remainder of the simulation. The addition of an ice-mélange backstress to an already retreated simulation (experiment 4, Table 1) allowed the glacier to undergo a winter readvance (Fig. 4f). However, the winter advance was not substantial enough to produce a different final terminus position compared to experiment 3 when no ice-mélange pressure was applied. Following the removal of the ice-mélange pressure, it only took the terminus 3 months to return to the stable retreated position.

3.3 Stable retreated position

All simulations where no ice mélange was present retreated to the same stable position (Fig. 4). Taking a closer look at the final terminus position after 2 years with no ice mélange, it can be seen that the glacier was grounded on its south side on an area of relatively shallow bed topography (Fig. 5a). Behind this portion of grounded ice was a large area of ungrounded ice. The CDL predicted no substantial crevassing across the terminus, resulting in the stable terminus position (Fig. 5b). On the south lateral margin, the ice column was modelled to be almost intact (Calving index = 1, H _crev = 0). However, a large area of extension was present slightly upstream where full thickness crevassing was modelled. This can be attributed to the basal crevassing shown in Figure 5c. This region does not calve as it was neither present on the terminus calving face nor can be evacuated as a calved iceberg. Across the full terminus, the majority of the crevasses modelled can be attributed to basal crevassing. Areas of the terminus which are grounded show little basal crevassing and no calving is predicted.

Figure 5. Analysis of the predicted crevassing at the end of a 2-year simulation with no ice mélange. (a) The grounding zone of the glacier terminus where red areas (values of 1) are grounded ice, blue areas (values of −1) are ungrounded ice and values of 0 are the grounding lines. (b) The predicted ‘Calving index’ which is 1 − H _crev and shows the percentage of intact ice. (c) The predicted penetration of basal crevasses. The basal crevasse index is H _crev, basal. (d) The predicted penetration of surface crevasses. The surface crevasse index is H _crev, surf. Note the small area of grounding on the left of the terminus that is suppressing basal crevassing and so calving. The N and S mark the north and south margins of the glacier. This is the same orientation as all other figures.

4. Discussion

The simulations presented here show that, when implemented in the new calving algorithm in Elmer/Ice, the CDL produces a close match to observed seasonal terminus variations of Jakobshavn Isbrae, once the calving law was adjusted through the crevasse penetration threshold and an appropriate ice mélange backstress was applied (Fig. 4). Similarly, the model replicated the observed calving behaviour seen at Jakobshavn Isbrae with full-thickness calving dominating the mass loss from calving. The results thus provide an opportunity to assess both the performance of the calving law and calving algorithm, and the relative importance of oceanic and topographic controls on the recent behaviour of Jakobshavn Isbrae. The high velocities compared to melt rates present at Jakobshavn Isbrae meant submarine melting was not considered important in modelling calving.

4.1 Evaluation of the crevasse-depth calving law

In this study, the crevasse-depth calving law was refined to optimise fit between model output and observed ice-front positions, in contrast with applications of the CDL at Store Glacier which used an unaltered model (Todd and others, 2018). This reflects differences between Store Glacier and Jakobshavn Isbrae, as well as differences between the calving algorithms that were used. First, the proportion of frontal ablation at Jakobshavn Isbrae resulting from calving (as opposed to melting) is much greater compared to Store Glacier (Mouginot and others, 2019). Therefore, the terminus position is more sensitive to the simulated calving magnitude, and any potential underestimation in calving will result in a more pronounced impact on the simulated terminus position. Second, the algorithm used by Todd and others (2018) had fixed lateral margins, in contrast with the new calving algorithm which allows the lateral margins to advance (Wheel and others, 2024). Importantly, the lateral margins provide the majority of support for the terminus of Greenlandic glaciers (Benn and others, 2023), so the new algorithm includes an important stabilising mechanism unavailable in the old model, reducing predicted calving. Combining these considerations, it is to be expected that implementation of the new calving model at Jakobshavn Isbrae will predict less calving compared to the work of Todd and others (2018) at Store Glacier. Although Todd and others (2018) represented the overall patterns of calving at Store Glacier well, Cook and others (2022) found the model underestimated calving, particularly small events, suggesting that some optimisation of the CDL could have resulted in further improvement of the results.

It is unsurprising that the CDL requires some tuning, for several reasons. First, the ‘zero-stress’ method for calculating crevasse depth is very simple and ignores factors such as fracture toughness of the ice. Second, crevasse depths are calculated based on the stresses in intact ice, and any stress concentrations that could lead to further fracture propagation are not included. Third, predicted crevasse depths are based on the instantaneous stress state, with no form of ice or crevasse history. Berg and Bassis (2022) showed a large increase in calving by adding a basic history in the CDL on a 2D Stokes synthetic geometry. In the study presented here, calving observations could be replicated for the 2016–2017 season at Jakobshavn Isbrae using a crevasse penetration threshold of 94.5% (Fig. 3). The differences in the modelled magnitude of winter maxima and the observations is due to the incapability of the model to fully capture the differing ice properties at the shear margins. As a consequence, the modelled glacier velocity is lower than that observed, resulting in a smaller winter readvance. Importantly, the seasonal advance and velocity patterns were captured by the model. Giving the missing processes not accounted for in the current implementation of the CDL, it is entirely reasonable that fracturing through 94.5% of the ice column would lead to full face calving. The adjustment of the CDL should not be considered to be a ‘tuned’ calving law but rather the adjustment increases sensitivity to account for missing processes. Ideally, this adjustment should be the same between locations and time periods but may need to be altered for numerical model parameters such as mesh resolution (Wheel and others, 2024). Further work is needed to investigate these processes and the potential for location dependency of the crevasse penetration threshold. Importantly, the use of the new 3D calving algorithm has made such work technically possible and could be combined with a full-stress analysis such as that described by Benn and others (2023).

Once adjusted, the calving law performed very well when compared to observations. The seasonal cycle was represented well with a large winter readvance, which is often difficult for calving laws to replicate (Benn and others, 2017). These results show the potential of the CDL as a candidate for further glacier modelling that will lead to improved understanding of glacier dynamics (Benn and others, 2023; Cook and others, 2023).

4.2 Importance of ice mélange

An ice-mélange backstress was essential to produce a simulated winter readvance of Jakobshavn Isbrae (Fig. 4). Without the presence of an ice mélange, the glacier continued to retreat throughout the winter months. Field estimations of ice-mélange backstress vary greatly (Walter and others, 2012), so values between 30 and 60 kPa were tested. Even for simulations where the calving law was very sensitive (93.5% CDL), an ice-mélange back pressure value of 60 kPa produced implausible winter terminus advance with a large peninsula forming (Fig. 2d), which are not seen in the observations. In contrast, values at the lower end of the Walter and others (2012) estimations produced more realistic results.

Although seasonal patterns are replicated well by adding a simple binary ice-mélange backstress (Fig. 3), there is a discrepancy of around 2 km between the maximum seasonal advance of the modelled terminus and observations (Fig. 3a). Furthermore, the modelled winter maxima occurs in May but observations suggest this occurs in April (Joughin and others, 2020). This can be directly attributed to the binary buttressing condition because calving occurs immediately once the ice mélange is removed. This implementation contrasts with the natural environment in which the ice mélange can go through several breakup and reforming events before becoming dislodged completely in summer (Cassotto and others, 2021). The granular and dynamic nature of ice mélange make more realistic dynamics difficult to incorporate in glacier modelling (Cassotto and others, 2021). It would be possible to include a sinusoidal pattern to the ice-mélange back pressure (Todd and others, 2018) but this does not capture the frequent breakup and reforming where any back pressure is removed entirely for short periods of time of less than a day (Cassotto and others, 2021). Additionally, where strong winter ice-mélange is present, it may be difficult to define the interface between calved ice and heavily crevassed glacier ice, a situation that will be impossible to resolve in a continuum model. Given these considerations, the incorporation of a simple binary buttressing force produces a satisfactory winter readvance and captures the suppression of calving extremely well (Fig. 3). Nevertheless, further work is needed to quantify ice-mélange backstress, especially considering most studies rely on the estimations based on one field campaign at Store Glacier (Walter and others, 2012).

The winter readvance acted as buffer allowing the glacier to advance beyond a stable minimum terminus position. As soon as the buttressing was removed, retreat ensued immediately. If the winter readvance, enabled by the presence of ice-mélange, was greater than the summer retreat then the glacier could maintain a more advanced position (Fig. 4). This produced a seasonal cycle which, from first viewing, suggests a stable glacier where annual terminus positions remain consistent (Figs 4a, c). If the ice mélange is removed, the terminus retreats continually until it reaches a stable position. Yet, once retreated, the addition of ice-mélange is not enough to allow a more advanced terminus position (Fig. 4f). Instead, the winter advance quickly calves off after a few months into the summer season. The results indicate that although the ice mélange can act as a buffer enabling the terminus to advance over unstable terrain, the ultimate retreated position of the glacier is defined by ice dynamics and geometry (Fig. 4).

4.3 Analysis of the stable position

For simulations where no ice mélange is present, the glacier does not go into a self-sustaining retreat. Instead, a stable position is reached towards the end of the first year that is maintained throughout the second year (Figs 4b, d). The retreat of around 2 km on the north side of the glacier is much more pronounced than the couple of hundreds of metres retreated on the south side. The larger retreat on the north side corresponded to an area of ungrounded ice (Fig. 2). Taking a closer look at the stable position reached after 2 years of simulation, an area of grounding on the glacier's south side suppresses simulated crevasse penetration and calving (Fig. 5). An area of extensional tensile stress sits behind this high point on the bed where full thickness basal crevassing occurs, but the ice here is unable to be evacuated as icebergs. A high point on the glacier bed also exists on the north side of the glacier (Fig. 6), indicating that the stable retreated position occupies a pinning point. Although uncertainties surrounding the bed elevation can lead to modelled features that do not exist in reality, this seems unlikely in this case given the size of the elevated bed area compared to the error of the bed product (Fig. 6; Morlighem and others, 2017). The model results are consistent with observational evidence that Jakobshavn Isbrae has occupied a consistent minimum position since ~2013 (Fig. 7), and we conclude that the stable minimum position in the model is the result of topographic control, and not time-dependent oceanographic factors such as the duration of ice mélange.

Figure 6. Bedrock elevation compared to modelled terminus positions at Jakobshavn Isbrae. The background is the bedrock height from BedMachine (Morlighem and others, 2017). The red line is the terminus position on 2017-04-29 just before the model ice mélange breaks up. The green/yellow line is the terminus position of the model setup. The grey line is the stable retreated position of the model taken at the end of the 2-year simulation with no ice mélange. Note the stable retreated terminus is positioned on bedrock highs at both lateral margins.

Figure 7. End of summer terminus positions from 2017 to 2023 with the modelled stable retreated position shown by the dashed line. The dates for each position are: 20-09-2017, 07-09-2018, 28-09-2019, 30-09-2020, 25-09-2021, 29-08-2022 and 01-10-2023. The background image is the Landsat-9 image from 01-10-2023.

The difference in terminus shape between the modelled stable position and observed summer positions is likely down to the coarse nature of the bed geometry (Figs 6 and 7). This study's results indicate that this stable position is independent of the presence of a strong ice-mélange. Instead, environmental conditions different to those in 2016–2017 (as used here) would be required to push the glacier beyond this point.

5. Conclusions

This study represents the first 3D Stokes modelling of Jakobshavn Isbrae. From a technical perspective, this study is a great step forward showcasing that unrestricted calving can be implemented with mesh adaptation at a large and dynamic glacier. The new Elmer/Ice calving model represents seasonal changes well, if a crevasse penetration threshold of 94.5% is incorporated into the crevasse-depth calving law. This adjustment compensates for missing processes such as stress concentrations or ice history, which are not currently incorporated in the calving law.

Winter readvance is only possible through the presence of ice-mélange backstress, which acts as a buffer that allows the glacier to advance beyond a stable position. Best-fit solutions were obtained using backstress of 30 kPa, at the low end of reported values. Once the ice-mélange is removed, the glacier calves towards a stable point but may not reach this position during one summer if there is a substantial winter advance. This allows the glacier to appear annually stable at a more advanced position while undergoing seasonal fluctuations. However, the current stable position of Jakobshavn Isbrae is determined by a combination of bed geometry and glacier dynamics rather than the presence of proglacial ice mélange. Further work could confirm this interpretation by a detailed stress analysis, similar to that of Benn and others (2023).