Modelling calving front dynamics using a level-set method: application to Jakobshavn Isbræ, West Greenland

Abstract

Calving is a major mechanism of ice discharge of the Antarctic and Greenland ice sheets, and a change in calving front position affects the entire stress regime of marine terminating glaciers. The representation of calving front dynamics in a 2-D or 3-D ice sheet model remains non-trivial. Here, we present the theoretical and technical framework for a level-set method, an implicit boundary tracking scheme, which we implement into the Ice Sheet System Model (ISSM). This scheme allows us to study the dynamic response of a drainage basin to user-defined calving rates. We apply the method to Jakobshavn Isbræ, a major marine terminating outlet glacier of the West Greenland Ice Sheet. The model robustly reproduces the high sensitivity of the glacier to calving, and we find that enhanced calving triggers significant acceleration of the ice stream. Upstream acceleration is sustained through a combination of mechanisms. However, both lateral stress and ice influx stabilize the ice stream. This study provides new insights into the ongoing changes occurring at Jakobshavn Isbræ and emphasizes that the incorporation of moving boundaries and dynamic lateral effects, not captured in flow-line models, is key for realistic model projections of sea level rise on centennial timescales.

Introduction

Calving of icebergs is a major mean of ice discharge for marine terminating glaciers around the world. It accounts for about half of the ice discharge of the Greenland and Antarctic ice sheets . This process causes calving front retreat, which leads to reduced basal and lateral resistive stress and results in upstream flow acceleration.

In order to assess the impact of calving on the dynamics of outlet glaciers using an ice sheet model, we need to include a dynamically evolving calving front. This requires tracking the calving front position and adjusting the boundary conditions accordingly. Addressing these issues is rather straightforward for 1-D flow-line or 2-D flow-band models , where the calving front is tracked along the flow line. However, this type of model lacks the consistent representation of lateral momentum transfer and lateral ice influx from tributaries for example, which have to be parameterized instead. This parameterization may neglect feedback effects important for simulations on decadal to centennial timescales, e.g. catchment area entrainment .

It is therefore critical to include a front tracking scheme in 2-D horizontal and 3-D models, which has been addressed by only a few ice sheet models e.g.. Various approaches to model the evolution of the shape of ice exist. Explicit methods track the position of a set of points, which represent the calving front. They require a complex technical framework to allow for geometric operations like folding and intersection of the continuum boundary, tracking singularities in curvature, and determining the position of a point in space relative to the modelled continuum.

Alternatively, the level-set method LSM; represents the continuum boundary implicitly by a contour, or “level set”, of a so-called “level-set function” (LSF). It easily handles topological changes of the modelled continuum, like splitting and merging. The LSM is based on a partial differential equation similar to the mass transport equation solved by ice sheet models. This makes the method straightforward to implement and allows for the application to continental-scale ice sheet simulations. Although the method does not necessarily conserve volume accurately, it is well established in continuum fluid mechanics . A LSM has been applied to represent the ice surface in 2-D flow-band models but not to model real ice sheets yet.

Observed ice surface velocities 2008/2009 of the Jakobshavn Isbræ drainage basin (logarithmic scale). Background image from Google Earth^©.

[Figure omitted. See PDF]

Understanding calving dynamics remains challenging because of the diversity of factors involved in calving events. Bathymetry, tides, and storm swell, as well as sea ice cover, ice mélange, and temperatures of both sea water and air are possible factors influencing calving rates. However, their effect, their respective share, and their interplay seem to vary from glacier to glacier and are not well understood . Therefore, no universal calving rate parameterization exists to date , and we rely here on user-defined calving rates. However, incorporating calving rate parameterization in the LSM should be straightforward.

Jakobshavn Isbræ is a major marine terminating glacier in West Greenland, which drains about 6.5 % of the Greenland Ice Sheet . It is characterized by two branches, which today terminate into a 30 km long ice-choked fjord (Figs.  and ). The southern branch exhibits high flow velocities, which are confined to a narrow, deep trough of about 5 km width. The trough retrogradely slopes inland to a maximum depth of about 1700 m below sea level and discharges most of the ice of the drainage basin . We refer to the fast-flowing area as “ice stream” and to the surrounding slow-moving ice as “ice sheet”. Those areas are separated by pronounced shear margins on either side of the ice stream.

Observations have shown that the fast flowing areas of Jakobshavn Isbræ exhibit a weak bed with a basal layer of temperate, soft ice . Basal sliding and shear in this layer cause most of these areas' horizontal motion. A large fraction of the ice stream's momentum is transferred to the adjacent ice sheet by lateral stress. It is thus well justified to use the 2-D shelfy-stream approximation SSA, to simulate this glacier.

Winter (Feb–Mar) calving front positions from 2009 to 2014 overlaid on a TerraSAR-X scene from 7 February 2015 (DLR^©). Dashed lines are used in case of ambiguous calving front positions.

[Figure omitted. See PDF]

Until the late 1990s, Jakobshavn Isbræ had a substantial floating ice tongue, which extended well into the fjord and was fed by both branches. The calving front position remained fairly constant from 1962 to the 1990s , and the glacier exhibited negligible seasonal variations in flow speed . In the 1990s, the glacier started a phase of acceleration, thinning, and retreat that followed the breakup of its ice tongue. Seasonal variations in calving front position and flow velocity increased sharply . Today, the glacier is one of the fastest ice streams in the world. It is still far from equilibrium and is a major contributor to global sea level rise . Observations suggest that the calving front position is a major control on the ice stream dynamics .

Various hypotheses have been proposed to explain the mechanisms behind this change. All identify the breakup of the floating ice tongue as the initial trigger of this dramatic chain of events, but different mechanisms have been proposed to explain the sustained acceleration, thinning, and retreat of the glacier. On the one hand, studies by and propose loss of buttressing and changes in basal conditions as the main cause behind the ongoing acceleration. On the other hand, argue that weakening of the lateral shear margins has significantly amplified the upstream acceleration. Several modelling studies of the glacier, which use 1-D flow-line and 2-D flow-band models, project unstable retreat of the glacier along its southern trough for up to 60 km inland within the next century . Other modelling studies argue that this type of ice stream is stable as long as it is fed by the surrounding ice sheet . However, numerical 2-D plan-view modelling efforts of Jakobshavn Isbræ so far lacked the representation of a dynamically evolving calving front. Hence, the hypotheses could not be tested in a satisfactory manner.

We present here a LSM-based framework to model the dynamic evolution of a calving front. This method is a step towards better physical representation of calving front dynamics in 2-D and 3-D ice sheet models. We describe the implementation of the method into the Ice Sheet System Model ISSM,, a parallel, state-of-the-art ice sheet model, and apply it here to Jakobshavn Isbræ in order to model its dynamic response to perturbations in calving rate.

Theory

Ice flow model

We employ the SSA on both floating and grounded ice. It neglects all vertical shearing but includes membrane stresses. The ice viscosity, $\mathit{\mu }$, follows Glen's flow law : $$\mathrm{2}\mathit{\mu }=B{\dot{\mathit{\epsilon }}}_{\text{e}}^{\frac{\mathrm{1}-n}{n}}.$$ Here, $n=\mathrm{3}$ is Glen's flow law coefficient, $B$ the ice viscosity parameter, and ${\dot{\mathit{\epsilon }}}_{\text{e}}$ the effective strain rate. We apply a Neumann stress boundary condition at the ice–air and ice–water interface, corresponding to zero air pressure and hydrostatic water pressure respectively. A linear friction law links basal shear stress, ${\mathit{\sigma }}_{\text{b}}$, to basal sliding velocity, ${\mathit{v}}_{\text{b}}$, on grounded ice: $${\mathit{\sigma }}_{\text{b}}=-{\mathit{\alpha }}^{\mathrm{2}}N{\mathit{v}}_{\text{b}},$$ where $\mathit{\alpha }$ denotes the basal friction parameter. We calculate the effective basal pressure, $N$, assuming that sea water pressure applies everywhere at the glacier base, which is a crude approximation far from the grounding line. The ice thickness, $H$, evolves over time according to the mass transport equation: $${\displaystyle \frac{\partial H}{\partial t}}=-\mathrm{\nabla }\cdot \left(H,\mathit{v}\right)+a_{\text{s}}+a_{\text{b}}.$$ Here, $\mathit{v}$ is the depth-averaged horizontal ice velocity, and $a_{\text{s}}$ and $a_{\text{b}}$ are the surface and basal mass balance respectively. We determine the grounding line position using hydrostatic equilibrium and treat it with a sub-element parameterization . We refer the reader to for details on the treatment of these equations in ISSM.

Symbols and model parameters.

Symbol	Quantity
$\mathit{\mu }$	Ice viscosity
$B$	Ice viscosity parameter
${\dot{\mathit{\epsilon }}}_{\text{e}}$	Effective strain rate
$n$	Glen's flow law parameter
$\mathit{\alpha }$	Basal friction parameter
$N$	Effective basal pressure
$\mathit{v}$	Depth-averaged horizontal ice velocity
$H$	Ice thickness
$a_{\text{s}}$	Surface mass balance
$a_{\text{b}}$	Basal mass balance
$\mathrm{\Omega }$	Computational domain
${\mathrm{\Omega }}_{\mathrm{i}}$	Ice domain
${\mathrm{\Omega }}_{\mathrm{c}}$	Ice-free domain
$\mathrm{\Gamma }$	Ice boundary
${\mathrm{\Gamma }}^{\mathrm{h}}$	Numerical ice boundary
$\mathit{\phi }$	Level-set function
$\mathit{n}$	Unit surface normal
$\mathit{w}$	Level-set velocity
$a^{⟂}$	Ablation rate
$c^{⟂}$	Calving rate
$m^{⟂}$	Melting rate
$s$	Scaling function
$p$	Perturbation function
$\mathrm{\Delta }t$	Perturbation duration
$p_{\mathrm{0}}$	Perturbation strength
$L$	Seasonal calving period length
${\mathit{\varphi }}_{\mathrm{0}}$	Phase shift
$Q_{\mathrm{cf}}$	Calving flux

Level-set method

Let $\mathrm{\Omega }$ be a computational domain in 2-D or 3-D space and $\mathit{\phi }$ a real, differentiable function on $\mathrm{\Omega }\times {\mathbb{R}}_{+}$, called LSF. For any $c\in \mathbb{R}$, we define the contour, or “$c$ level set”, of $\mathit{\phi }$ as $\mathit{\phi }(\mathbf{x},t)=c$. Taking its material derivative yields the “level-set equation” (LSE): $${\displaystyle \frac{\partial \mathit{\phi }}{\partial t}}+\mathit{w}\cdot \mathrm{\nabla }\mathit{\phi }=0.$$ This Hamilton–Jacobi type partial differential equation describes how level sets move with the local value of the velocity $\mathit{w}$, which is called level-set velocity. We need to provide an initial condition ${\mathit{\phi }}_{\mathrm{0}}\left(\mathbf{x}\right)=\mathit{\phi }(\mathbf{x},t=\mathrm{0})$ to solve the LSE.

We use $\mathit{\phi }$ to partition $\mathrm{\Omega }$ into three disjoint subdomains: the ice domain, ${\mathrm{\Omega }}_{\mathrm{i}}\left(t\right)$, its complement, ${\mathrm{\Omega }}_{\mathrm{c}}\left(t\right)$, and their common boundary, $\mathrm{\Gamma }\left(t\right)$, $$\left\{\begin{array}{rl}\mathit{\phi }(\mathbf{x},t)<\mathrm{0}& \iff \mathbf{x}\in {\mathrm{\Omega }}_{\mathrm{i}}\left(t\right)\\ \mathit{\phi }(\mathbf{x},t)=\mathrm{0}& \iff \mathbf{x}\in \mathrm{\Gamma }\left(t\right)\\ \mathit{\phi }(\mathbf{x},t)>\mathrm{0}& \iff \mathbf{x}\in {\mathrm{\Omega }}_{\mathrm{c}}\left(t\right).\end{array}\right.$$ We omit the time dependence of these sets in the remainder of this article. By construction, $\mathrm{\Gamma }$, the 0 level set of $\mathit{\phi }$, separates ${\mathrm{\Omega }}_{\mathrm{i}}$ and ${\mathrm{\Omega }}_{\mathrm{c}}$.

We extend $\mathit{n}$, the outward-pointing unit normal of $\mathrm{\Gamma }$, onto $\mathrm{\Omega }$ using the LSF by $$\mathit{n}={\displaystyle \frac{\mathrm{\nabla }\mathit{\phi }}{\left|\mathrm{\nabla }\mathit{\phi }\right|}}.$$ For details on the level-set method and its applications, we refer to and .

Schematic of the numerical ice margin. The red dashed lines denote different contour lines (level sets) of the LSF $\mathit{\phi }$. The thick red line marks the 0 level set, $\mathrm{\Gamma }$, and the yellow line the numerical calving front ${\mathrm{\Gamma }}^{\text{h}}$. Dark blue triangles are ice-free elements, white ones are ice-filled, and the light blue ones are the front elements. The three vectors show an example of the level-set velocity $\mathit{w}=\mathit{v}+\mathit{c},\mathit{c}=-c^{⟂}\mathit{n}$ at a finite element node.

[Figure omitted. See PDF]

The boundary position of an ice sheet evolves with the sum of the ice velocity $\mathit{v}$ along $\mathit{n}$ and an ablation rate, $a^{⟂}$: $$\mathit{w}\cdot \mathit{n}=\mathit{v}\cdot \mathit{n}-a^{⟂}.$$ It follows that the ice boundary is stationary if and only if $a^{⟂}=\mathit{v}\cdot \mathit{n}$, i.e. the ablation rate matches the ice velocity perpendicular to the ice boundary. The ablation rate needs to be prescribed, either based on observations or through a parameterization.

Note that no limitations have been made so far with respect to the dimension of the problem in this section. Accordingly, the method could be applied to model the evolution of the glacier thickness and lateral extent simultaneously . However, in this 2-D plan-view model study, we use the LSM to model only the horizontal extent of the ice sheet. Its vertical extent is described by the mass transport Eq. ().

We model lateral ablation as the sum of a melting rate, $m^{⟂}$, and a calving rate, $c^{⟂}$: $a^{⟂}=c^{⟂}+m^{⟂}$. For simplicity, we assume in the remainder of the article that lateral ablation occurs in the form of calving exclusively, i.e. $m^{⟂}=\mathrm{0}$. Calving itself is assumed to be a quasi-continuous process, consisting of frequent, but small, calving events. With Eqs. () and (), the LSE (Eq. 4) becomes $${\displaystyle \frac{\partial \mathit{\phi }}{\partial t}}+\mathit{v}\cdot \mathrm{\nabla }\mathit{\phi }=c^{⟂}\left|\mathrm{\nabla }\mathit{\phi }\right|,$$ which is also known as “kinematic calving front condition” KCFC; . Both the calving rate and ice velocity need to be provided on the entire 2-D computational domain $\mathrm{\Omega }$ in order to solve the KCFC. An example of a calving rate field will be given in Sect.  and is shown in Fig. . The KCFC implies that all level sets of $\mathit{\phi }$, including the calving front $\mathrm{\Gamma }$, move at a given time with the local sum of the horizontal ice velocity and calving rate along the normal $\mathit{n}$ (Fig. ). We define the calving flux, $Q_{\mathrm{cf}}$, as the ice flux crossing the calving front: $$Q_{\mathrm{cf}}=\underset{\mathrm{\Gamma }}{\int }{c}^{⟂}\left(r\right)H\left(r\right)\text{d}r.$$

The ice velocity components and the ice thickness are only defined on ${\mathrm{\Omega }}_{\mathrm{i}}$ and need to be extended onto ${\mathrm{\Omega }}_{\mathrm{c}}$ (see also Sect. ). Any scalar field, $S$, is extrapolated onto ${\mathrm{\Omega }}_{\mathrm{c}}$ by solving $$\mathit{n}\cdot \mathrm{\nabla }S=\mathrm{0},$$ while keeping $S$ fixed on ${\mathrm{\Omega }}_{\mathrm{i}}$. This type of extrapolation has the tendency to preserve $\left|\mathrm{\nabla }\mathit{\phi }\right|=\mathcal{O}\left(\mathrm{1}\right)$ when we use the extrapolated ice velocity field to solve the KCFC .

Implementation

ISSM relies on the finite element method (FEM) to solve partial differential equations. It applies a continuous Galerkin FEM using triangular (2-D) and prismatic (3-D) Lagrange finite elements and uses anisotropic mesh refinement to limit the number of degrees of freedom while maximizing spatial resolution in regions of interest.

We discretize the KCFC (Eq. ) and extrapolation equation (Eq. ) using linear finite elements on the same mesh as the one used to model the ice dynamics. We stabilize both equations with artificial diffusion , which after thorough testing proved to be the most robust stabilization scheme. We integrate over time using a semi-implicit time-stepping scheme. We solve the KCFC and the field equations for ice flow modelling in a decoupled fashion. The KCFC is solved first with input data from the previous time step. We then update the numerical ice domain using the new LSF as described below and update boundary conditions accordingly. Finally, we solve the momentum balance and the mass transport equation on the updated ice domain.

The 0 level set of $\mathit{\phi }$, $\mathrm{\Gamma }$, does not in general coincide with the finite element mesh edges due to its implicit representation. It intersects a number of elements (“front elements”) with a hyperplane, which divides them into an ice-filled and an ice-free part (Fig. ). This has various implications on the numerical level. When assembling the system stiffness matrices for ice flow modelling, exclusive integration over the ice-filled part of the element would be required. The stress boundary condition at the calving front would have to be applied at the intersecting hyperplane. Currently, ISSM is not capable of resolving those submesh-scale processes.

Therefore, we either fully activate or deactivate a mesh element at every time step. Only active elements are considered for the numerical discretization of the respective field equations. We activate an element if at least one of its vertices is in ${\mathrm{\Omega }}_{\mathrm{i}}$ or $\mathrm{\Gamma }$, and the element is then considered to be entirely filled with ice. We flag the element as ice free if it lies entirely inside ${\mathrm{\Omega }}_{\mathrm{c}}$, and it is deactivated. As a consequence, the numerical calving front, ${\mathrm{\Gamma }}^{\mathrm{h}}$, runs along mesh edges and is updated in a discontinuous manner (Fig. ). We apply the stress boundary condition along ${\mathrm{\Gamma }}^{\mathrm{h}}$ for numerical consistency. Calving front normals on $\mathrm{\Gamma }$ and ${\mathrm{\Gamma }}^{\mathrm{h}}$ may differ significantly in direction. However, stress components tangential to $\mathit{n}$ cancel out along ${\mathrm{\Gamma }}^{\mathrm{h}}$, so that the integrated stress exerted at the calving front is close to the one applied along $\mathrm{\Gamma }$. For all further calculations where a normal is involved, like extrapolation, the normal to the LSF (Eq. ) is used.

The numerical calving front is by definition further downstream than $\mathrm{\Gamma }$. This may lead to slightly higher resistive lateral stress at the calving front, the magnitude of which depends on the excess ice area of the intersected front element and the front geometry. We use a fine mesh resolution in the vicinity of the calving front to limit this effect. We extrapolate the calving front thickness onto the ice-free domain using equation (Eq. ). This yields realistic ice thickness and ice thickness gradients across the front elements that would otherwise lead to overestimated driving stress and underestimated water pressure at the ice–ocean interface. If not corrected, those two effects unrealistically increase ice velocities at the calving front, which then feed back into the mass transport and LSM schemes.

Calving rate field $c_{\mathrm{0}}^{⟂}$ in the region of fast flow, derived from modelled ice velocities at the end of the relaxation run. The red line indicates the 0 level set of the initial LSF used for geometry relaxation and as start position for the calving front during the experiments. The turquoise line marks the grounding line. Purple contours indicate zero bedrock elevation. Black lines are the “along-trough” (a) and “across-trough” (b) profiles used in Fig. . The finite element mesh is displayed in grey.

[Figure omitted. See PDF]

We present two experiments for validation of the LSM implemented here in the Appendix. The first experiment shows that the ice margin is advected with the prescribed level-set velocity $\mathit{w}$. The linear representation of the LSF on an unstructured mesh causes a small error in the exact level-set position, which depends on element size and cancels out over time. The second test shows that errors in volume conservation introduced by the LSM decrease with finer mesh resolution and are below 0.2 % after 100 years for a mesh resolution of 1 km. In the application to Jakobshavn Isbræ, we use a front element size of 0.5 km. The potential volume loss inherent to this implementation of the LSM is thus far below current uncertainties of model input data.

Inclusion of the LSM requires additional computational effort for the extrapolation of field variables, to solve the KCFC and for extra iterations of the momentum balance solver, since the stress boundary conditions at the calving front change frequently. Its amount depends on the flow approximation and especially on whether the model setup is close to a stable configuration or not. Using SSA, the additional computational cost reaches up to 25 %, of which 11 % is caused by the solution of the KCFC.

Data and model setup

Jakobshavn Isbræ model setup

We use Jakobshavn Isbræ's drainage basin from to generate a 2-D horizontal finite element mesh with element size varying from 500 m in the fjord and areas of fast flow to 10 km inland (Fig. ). We choose this high mesh resolution to minimize calving front discretization errors, and to resolve the fjord and the deep trough accurately in the model. The resulting mesh has about 10 000 vertices and 19 000 elements. Due to high flow velocities, the Courant–Friedrich–Lewy condition CFL; dictates a time step on the order of days for the solution of the momentum balance equation, the mass transport equation, and the KCFC.

We use bed topography from , derived using a mass conservation approach . The ice surface elevation is taken from GIMP , and ice thickness is the difference between ice surface and ice base elevation. Bathymetry of the ice-choked fjord of Jakobshavn Isbræ is difficult to measure and currently poorly known. As a first-order estimate, we apply a parabolic profile of 800 m depth along the ice fjord, fitted via spline interpolation to known topography data. We rely on for the surface mass balance and, as a first approximation, use their surface temperatures to initialize our ice viscosity parameter, $B$, based on the table from . Surface temperatures range from $\sim -\mathrm{8}$ ${}^{\circ }$C near the coast to $\sim -28$ ${}^{\circ }$C near the divide. Basal mass balance is set to 0 and no thermal model is run. All these forcings are kept constant over time.

We infer a basal friction coefficient, $\mathit{\alpha }$, in Eq. () using an adjoint-based inversion of surface velocities from 2009 . In regions like the fjord, where there is no ice today, we apply an area-averaged value of $\mathit{\alpha }=30{\text{ a}}^{\mathrm{1}/\mathrm{2}}{\text{m}}^{-\mathrm{1}/\mathrm{2}}$. At the margins of the computational domain we prescribe zero horizontal ice velocities in order to prevent mass flux across this boundary. The friction parameter is kept fixed over time for all model simulations.

Calving front positions for experiments A, B1, B2, and C4 at the start of each year, plotted over basal topography (grey).

[Figure omitted. See PDF]

Inconsistencies in model input data cause sharp readjustments of the glacier state at the beginning of each simulation, which would make it difficult to distinguish between such effects and those of the applied forcing . Therefore, we relax the model prior to the experiments using a fixed, piecewise linear LSF ${\mathit{\phi }}_{\mathrm{0}}$, whose 0 level set corresponds to the mean annual calving front position of 2009 (Fig. ). Since the glacier in this configuration is far from steady state, model relaxation causes considerable thinning across the glacier's catchment area. In order not to deviate too much from present day's geometric setting we choose a 100-year relaxation time period. Note that the grounding line retreats during the relaxation due to dynamic thinning, so that the glacier forms a new floating ice tongue. This ice tongue extends about $15\text{km}$ to a local topographic maximum in the southern trough and $\mathrm{3}\text{km}$ into the northern one (Fig. ). The relaxed geometry constitutes the initial state for our experiments. Due to this deviation in geometry, providing quantitative insights into Jakobshavn Isbræ is beyond the scope of this study. However, the main characteristics of the ice stream (e.g. its large drainage basin and the narrow outlet channel) are preserved, so that the results presented in this paper qualitatively represent the behaviour of Jakobshavn Isbræ.

Description of experiments

We set $$c_{\mathrm{0}}^{⟂}=q\left|{\mathit{v}}_{\mathrm{0}}\right|$$ as a basic calving rate estimate, where ${\mathit{v}}_{\mathrm{0}}$ is the velocity field at the end of the geometry relaxation run, extended onto ${\mathrm{\Omega }}_{\mathrm{c}}$ (Fig. ). The continuous function $q$ is equal to $\mathrm{1}$ in areas where the bed lies below $-$300 m, and it linearly drops to 0 in areas of positive bed elevation. It prevents calving occurring in areas with a glacier bed above sea level, as suggested by observations of tidewater glaciers . The choice of the calving rate estimate is motivated by the small observed angle between ${\mathit{v}}_{\mathrm{0}}$ and $\mathit{n}$ at the calving front (${\mathit{v}}_{\mathrm{0}}\approx \left|{\mathit{v}}_{\mathrm{0}}\right|\mathit{n}$). Then ${\mathit{w}}_{\mathrm{0}}\cdot \mathit{n}={\mathit{v}}_{\mathrm{0}}\cdot \mathit{n}-c_{\mathrm{0}}^{⟂}\approx \mathrm{0}$, so that we can expect this calving rate estimate to yield a stationary calving front, if applied to a geometry that is in steady state. We scale $c_{\mathrm{0}}^{⟂}$ over time with a scaling function, $s$, which allows for the representation of seasonal cycles, and a perturbation function, $p$, to modify the calving rate for some period of time. The applied calving rate is then $c^{⟂}(\mathbf{x},t)=c_{\mathrm{0}}^{⟂}\left(\mathbf{x}\right)s\left(t\right)p\left(t\right)$.

We perform three suites of experiments in order to analyse the impact of the calving rate on the glacier's dynamics. The calving front is now allowed to freely evolve in response to $c^{⟂}$. All experiments run for $120$ years.

In experiment A, we keep the calving rate constant over time, i.e. we set both $s\left(t\right)=p\left(t\right)=\mathrm{1}$. Hence, $c^{⟂}(\mathbf{x},t)=c_{\mathrm{0}}^{⟂}\left(\mathbf{x}\right)$. This experiment, although not physically motivated, is used to evaluate whether a stable calving front position can be reached using the LSM and for comparison to the experiments described below.

In experiment suites B and C, we represent the seasonal cycle by scaling $c_{\mathrm{0}}^{⟂}$ by $s\left(t\right)=max(\mathrm{0},\mathit{\pi }$ $\sin \left(\mathrm{2},\mathit{\pi },,(,t,/,L,-,{\mathit{\varphi }}_{\mathrm{0}},)\right))$, with a phase shift ${\mathit{\varphi }}_{\mathrm{0}}=\mathrm{4}/12$ and a period $L=\mathrm{1}\text{ a}$. We perturb the calving rate during a limited duration, $\mathrm{\Delta }t$, with a perturbation strength $$p_{\mathrm{0}}\ge \mathrm{0}:p\left(t\right)=\mathit{}\left\{\begin{array}{ll}{p}_{\mathrm{0}},& \text{if}{t}_{\mathrm{0}}<t<t_{\mathrm{0}}+\mathrm{\Delta }t,\text{and}\\ \mathrm{1},& \text{else.}\end{array}\right.$$ We start the perturbation at $t_{\mathrm{0}}=20\text{a}$ for all experiments. In experiment suite B, we perform five experiments with $\mathrm{\Delta }t=\mathrm{1}\text{a}$, while varying $p_{\mathrm{0}}$ from 0 to 4 by increments of 1. In experiment suite C, we keep $p_{\mathrm{0}}=\mathrm{2}$ fixed and set $\mathrm{\Delta }t$ as 2, 4, and 8 years. We use the notation B $<p_{\mathrm{0}}>$ and C $<\mathrm{\Delta }t>$ to identify single experiments, e.g. B1 for experiment B with perturbation strength $p_{\mathrm{0}}=\mathrm{1}$, which represents the case of unperturbed periodic calving. B1 is used as the control run to which the other experiments are compared. Table  lists all the experiments performed here.

Table of experiments.

Name	$p_{\mathrm{0}}$	$\mathrm{\Delta }t$	Name	$p_{\mathrm{0}}$	$\mathrm{\Delta }t$	Name	$p_{\mathrm{0}}$	$\mathrm{\Delta }t$
A	1	0	B2	2	1	C2	2	2
B0	0	1	B3	3	1	C4	2	4
B1	1	1	B4	4	1	C8	2	8

Results

Figure  shows calving front positions for experiments A, B1, B2, and C4. Under constant calving rate forcing, the calving front remains at a stable position after minor readjustments in the first decade of the simulation. In experiment A, the calving front undergoes gradual retreat over time due to the slowdown of the glacier caused by its ongoing thinning. When we perturb the calving rate, the calving front migrates, and higher calving rates lead to larger retreats. The retreat is highest in areas of fast flow and strongly decreases towards the ice stream margins. This yields the characteristic concave shape of a retreating calving front. The modelled calving front positions and their shape are in good agreement with observations (Fig. ). The retreat rate during continued phases of calving decreases to 0, so that the calving front reaches a new stable position 9 km upstream of its initial position (Fig. d). In experiments B and C, the calving front returns to a similar position as in the unperturbed experiment B1 within 10 to 20 years after the perturbation stops.

Section profiles of ice geometry (a) and ice velocity (b) along-trough, as well as ice geometry (c) and effective strain rates (d) across-trough for experiment C4 at the end of the calving season in October each year. Positions of the lines are given in Fig. .

[Figure omitted. See PDF]

Figure  shows ice velocity, geometry, and strain rates for experiment C4 along two lines, which go along and across the southern trough respectively (Fig. ). During the first 20 years prior to the perturbation, the ice thickness in the floating part decreases by about 100 m (Fig. a). As the calving front retreats during the perturbation, the ice velocity increases (Fig. b), and the ice thickness adjusts accordingly (Fig. a). The ice thinning leads to a fast retreat of the grounding line in the regions of locally retrograde bed and temporarily stabilizes over local along-trough topographic maxima, referred to as “local highs”. The southern trough has many local highs, which act as pinning points and are critical for flow dynamics, in agreement with earlier results from . The acceleration of the ice stream extends tens of kilometres upstream to areas of grounded ice (Fig. b). Thinning and acceleration are strongest over the ice stream and spread out to the surrounding ice sheet in a dampened fashion. These thinning and acceleration patterns increase surface and velocity gradients, especially in the shear margins (Fig. c), where the effective strain rates gradually increase up to a factor of 4 in experiment C4 (Fig. d), which corresponds to a drop in viscosity of about 60 % (Eq. ). This substantially weakens the mechanical coupling between the ice stream and the surrounding ice sheet.

Calving front and grounding line positions along-trough (left column), calving front thickness, ice velocity, and effective strain rate relative to their initial value along-trough (right column) over time for experiments A, B1, B2, and C4. Perturbation intervals are marked in grey. Relative values for ice velocity have been shifted up by 0.5 for better visibility (red $y$ axis).

[Figure omitted. See PDF]

Figure  shows the intra-annual variability of ice properties at the calving front and grounding line for experiments A, B1, B2, and C4. All shown variables reflect the characteristics of the applied calving rate forcing. The constant calving rate applied in experiment A leads to a steady configuration (Fig. a and b). For an unperturbed periodic calving rate forcing (Fig. c and d), the calving front position oscillates around a constant annual mean value by $\pm$3 km, while the grounding line position remains unchanged at kilometre 29. Ice velocities and thickness at the calving front act in phase with the calving front position, while the response of strain rates at the grounding line is delayed by about a month. The ice velocity varies over a year by $\pm$20 %, which corresponds to about $\pm$2 km a${}^{-\mathrm{1}}$, the ice thickness by $\pm$13 %, or $\pm$100 m, and effective strain rates at the grounding line by $\pm$7 %, or $\pm$0.1 a${}^{-\mathrm{1}}$.

The response to a calving rate perturbation scales with $p_{\mathrm{0}}$ and $\mathrm{\Delta }t$. When the calving rate doubles (B2, C4), the calving front retreats initially at an average rate of 4.5 km a${}^{-\mathrm{1}}$. The calving front stabilizes 9 km upstream for longer perturbations (Fig. g). The intra-annual variability of the calving front position doubles to $\pm$6.5 km. The grounding line position is hardly affected by small calving rate perturbations, but large perturbations trigger fast retreats of several kilometres, which in turn cause drastic, but short-lived, flow accelerations (Fig. g and h). The annual average ice velocity increases by 10 %, and its intra-annual variability doubles to $\pm$38 % (Fig. h). The mean calving front thickness decreases by 30 % towards the end of the perturbation of experiment C4 and experiences large variations up to $\pm$75 %. This high thickness variability is due to the front retreating into areas of thick ice in summer followed by stretching and thinning during calving front advance in winter. For small perturbations, variations of effective strain rates at the grounding line quadruple to $\pm$25 % (Fig. f). Once the calving rate perturbation stops, all variables display remarkable reversibility.

When calving is temporarily turned off (experiment B0, not shown here), the response of the glacier is reversed: the calving front advances, creating a convex ice tongue. Meanwhile, the ice stream decelerates, thickens, and the grounding line advances. After the perturbation, the glacier retreats into a state slightly thicker and faster than the one of experiment B1.

(a) Absolute difference in ice volume for the different simulations with respect to experiment B1. The non-oscillating ice volume profile of experiment A causes its difference to experiment B1 to oscillate. (b) The volume differences from experiments B and C divided by $\mathrm{\Delta }t(\mathrm{1}-p_{\mathrm{0}})$, the measure of the time-integrated calving rate perturbation.

[Figure omitted. See PDF]

Figure a shows the evolution of the ice volume with respect to experiment B1, the control run. The glacier in experiment B1 continues to lose volume at an average rate of $-22.8{\text{km}}^{\mathrm{3}}{\text{a}}^{-\mathrm{1}}$ due to the ongoing geometry relaxation. In experiment A, Jakobshavn Isbræ loses an additional $0.4{\text{km}}^{\mathrm{3}}{\text{a}}^{-\mathrm{1}}$, which corresponds to the gradual retreat of the calving front. Enhanced calving causes additional volume loss proportional to $\mathrm{\Delta }t(\mathrm{1}-p_{\mathrm{0}})$, the measure of the time-integrated calving rate perturbation (Fig. b). If the calving rate is doubled, the additional volume loss reaches $-35.7{\text{ km}}^{\mathrm{3}}{\text{a}}^{-\mathrm{1}}$ in the first year but decreases with time, as the calving front thins and retreats into areas of lower calving rates. Those numbers agree well with recent ice discharge observations . Over the first decade after the perturbation, all modelled glaciers recover 40 to 60 % of the volume deviation to the control run.

Discussion

The applied calving rate determines the behaviour of the calving front and the ice stream. In our simulations, larger perturbation strengths $p_{\mathrm{0}}$ lead to faster calving front retreats. In the case of long perturbations (experiments C4 and C8), the calving front reaches a new stable position. A stable calving front position requires the calving rate to be larger than the ice velocity if the calving front advances; similarly, if the calving front retreats, the calving rate needs to be lower than the ice velocity.

Several mechanisms determine how the model responds to the calving rate forcing. First, a calving rate increase leads to a retreat of the calving front position, ice stream acceleration and dynamic thinning in the vicinity of the terminus. Second, this dynamic thinning increases surface slopes and therefore the driving stress. As the glacier locally speeds up, the ice thinning propagates upstream. The ice stream thins much faster than the surrounding ice sheet, which steepens the surface across the shear margins. Lateral inflow of ice into the ice stream hence increases until it balances the calving flux. This limits the thinning of the ice stream. Thinning of the ice stream in turn leads to grounding line retreat and reduction in basal effective pressure, both of which reduce basal drag significantly in the vicinity of the grounding line. We showed that grounding line retreat leads to short-lived but drastic increases in ice flux. This mechanism is qualitatively the same as the one described in and . Several pinning points along the retrograde trough of the southern branch, as well as the lateral stress transfer and mass influx, prevent the modelled ice stream from being prone to the marine ice sheet instability , a hypothesis which states that grounding line positions are unstable on retrograde slopes. This corroborates earlier results by , who presented examples of stable grounding line positions on retrograde beds. However, due to large uncertainties in the input data, and since some physical processes are not represented in our experiments, evaluation of this question for Jakobshavn Isbræ is beyond the scope of this study.

A third mechanism is related to the calving front lengthening during its retreat (e.g. Fig. ). The lengthening causes tributaries of the main ice stream to calve directly into the fjord, thereby increasing the calving flux $Q_{\mathrm{cf}}$ (Eq. ) and thinning of the terminus vicinity.

Finally, the ice stream accelerates faster than the surrounding ice sheet, which increases strain rates at the shear margins. This reduces the ice viscosity in these areas, which mechanically decouples the ice stream from the ice sheet, allowing the ice stream to accelerate further. This positive feedback confines the initial thinning to the ice stream and is controlled by the rate at which ice enters the ice stream. This mechanism is essential for enabling ice stream acceleration tens of kilometres upstream of the grounding line, since large fractions of the ice stream's driving stress are balanced by lateral stress. This corroborates force balance arguments produced earlier by .

In experiments A and B1, we apply the same annual mean calving rate. However, due to the lack of seasonal cycle in calving rate the mechanical coupling between the ice stream and ice sheet is higher in experiment A. The ice stream velocity is therefore lower, causing gradual calving front retreat and additional ice volume loss. This illustrates that volume change estimates from models with and without seasonal cycles of calving may differ. Our results suggest that including both a dynamically evolving calving front as well as seasonal cycles are critical for accurate projections of future contributions of ice sheets to global sea level rise on decadal to centennial timescales.

Response mechanisms not covered here will likely include feedbacks with damage mechanics and thermodynamics due to the increased strain rates. During longer perturbations, ice surface lowering will probably affect the surface mass balance and the drainage basin outline.

The reversibility of the calving front configuration after the calving rate perturbation is a robust feature across all experiments. The short duration of the perturbation, the prescribed calving rates, and the geometry of the glacier are responsible for this behaviour. The volume change in all experiments never exceeds 0.1 % of the initial glacier volume in the experiments shown here. Once the perturbation stops, the surrounding ice sheet continues to replenish the ice stream, which allows for its quick recovery.

The modelled glacier response to enhanced calving is in good qualitative agreement with observations, which corroborates that calving is a major control on this glacier. The similar shape of the modelled and observed calving front suggests that calving rates are indeed proportional to its flow speed during the glacier's current retreat. However, the reversibility of the modelled calving front position is in contrast to Jakobshavn Isbræ's actual behaviour. Sustained high calving rates are therefore necessary to explain the continued retreat of the glacier, as our results suggest that the glacier would have readvanced otherwise. Accurate model input data, representation of all relevant physical processes, and incorporation of a suitable calving rate parameterization will be necessary for quantitative analysis of this dynamic ice stream.

Conclusions

In this study, we present the theoretical framework for coupling a LSM to ice dynamics and implement it into ISSM. The LSM proves to be a robust method for modelling the dynamic evolution of a calving front. We apply this technique to Jakobshavn Isbræ using prescribed calving rates, and we find that the glacier is highly sensitive to this forcing, which agrees well with observations.

Calving rate perturbations strongly affect the ice stream through several linked mechanisms. First, changes in calving rate cause calving front migration and alter the ice discharge. Second, the resulting thickness change at the calving front spreads out to the surrounding ice sheet. Third, thinning-induced grounding line retreat causes further ice stream acceleration and creates a positive feedback. Finally, shear margin weakening caused by the ice stream acceleration decreases lateral drag resisting ice flow. This positive feedback mechanism sustains significant acceleration of the ice stream tens of kilometres upstream of the grounding line.

The surrounding ice sheet is barely affected by short periods of enhanced calving. It stabilizes the ice stream and allows for quick reversibility of the calving front position through lateral ice influx and stress transfer once the calving rates are set back to their initial values.

Since the calving front position and dynamic lateral effects are critical to simulate and understand the behaviour of marine terminating glaciers, the inclusion of moving boundaries in 2-D plan-view and 3-D models is key for realistic sea level rise projections on centennial timescales. This method is a step towards better physical representation of calving front dynamics in ice sheet models.

Validation of the level-set method

We present two simple test setups to validate the LSM. The first is designed to show the accurate advection and shape preservation properties of the method. The second setup aims to give an estimate for the volume change introduced by the LSM for different mesh resolutions.

Advection

Let $\mathrm{\Omega }$ be a 50 km square with the initial LSF as $${\mathit{\phi }}_{\mathrm{0}}\left(\mathbf{x}\right)=\Vert \mathbf{x}-{\mathbf{x}}_{\mathrm{0}}{\Vert }_{\mathrm{2}}-R,$$ where ${\mathbf{x}}_{\mathrm{0}}=(25,25)$ km and $R=12.5$ km, so that the initial 0 level set describes a circle in the middle of the domain. We prescribe a constant velocity $\mathit{v}=\left(\cos ,(,\mathit{\pi },/,\mathrm{4},),,,\sin ,(,\mathit{\pi },/,\mathrm{4},)\right)$ km a${}^{-\mathrm{1}}$ everywhere. We advect ${\mathit{\phi }}_{\mathrm{0}}$ over 10 years with time steps of $0.1$ a and keep track of the 0 level set.

Figure  shows the 0-level-set position at the beginning of every year. The LSM preserves the initial circular shape and can be used to model both advance and retreat of a calving front. We measure the advection speed of the 0 level set along the diagonal marked in white in Fig. . Figure  shows the standard deviation of the numerical error relative to the prescribed advection speed taken over time for different element sizes. The numerical error is due to the linear interpolation of the curved shape, which causes variations of the level-set velocity around the prescribed value. The standard deviation of the error linearly decreases with mesh resolution and drops below 1 % for elements sizes below 0.5 km. We therefore choose a mesh resolution of 0.5 km in the vicinity of the calving front in our simulations.

Volume conservation

Let $\mathrm{\Omega }$ be a $200\times 20{\text{ km}}^{\mathrm{2}}$ rectangle with an initial LSF given by $${\mathit{\phi }}_{\mathrm{0}}\left(\mathbf{x}\right)=\left(\mathrm{1},,,\mathrm{0}\right)\cdot \mathbf{x}-100\text{km}.$$ The initial lateral extent is thus a $100\times 20{\text{ km}}^{\mathrm{2}}$ rectangle. The geometry corresponds to the ice shelf ramp presented in . The ice thickness linearly decreases from 400 m at the grounding line ($x=\mathrm{0}$ km) to 200 m at the calving front ($x=100$ km). We apply zero surface accumulation and basal melt, as well as zero grounding line velocity and free slip boundary conditions at $y=\mathrm{0}$ km and $y=20$ km. The ice sheet spreads under its own weight for 100 years.

Figure  shows the evolution of the ice volume for different element sizes. All simulations show volume loss due to the free flux boundary condition at the numerical ice front, which is not entirely balanced by the volume added through the ice thickness extrapolation. The volume loss decreases with element size and is below 0.2 % of the initial ice volume after 100 years for an element size of 1 km. This volume loss is far below current uncertainties of other model input data.

Zero-level-set positions at the start of every year, plotted over ${\mathit{\phi }}_{\mathrm{0}}$, which is in grey scale. An example of the mesh with element size 1 km is marked in black. The white diagonal marks the line along which the velocity of the 0 level set is tracked.

[Figure omitted. See PDF]

Standard deviation of the relative numerical error in advection velocity of the 0 level set depending on element size.

[Figure omitted. See PDF]

Evolution of the relative ice volume change for different element sizes. The red line shows volume conservation.

[Figure omitted. See PDF]

Acknowledgements

A. Humbert acknowledges support of the DLR proposal HYD2059 which provides TerraSAR-X data for the project HGF-Alliance Remote Sensing and Earth System Dynamics. H. Seroussi, M. Morlighem, and E. Y. Larour are supported by grants from the National Aeronautics and Space Administration, Cryospheric Sciences, and Modeling, Analysis and Prediction programs. The authors thank the referees G. Jouvet and J. Bassis as well as the editor O. Gagliardini for their helpful and insightful comments. Edited by: O. Gagliardini