Controls on calving at a large Greenland tidewater glacier: stress regime, self-organised criticality and the crevasse-depth calving law

Abstract

We investigate the physical basis of the crevasse-depth (CD) calving law by analysing relationships between glaciological stresses and calving behaviour at Sermeq Kujalleq (Store Glacier), Greenland. Our observations and model simulations show that the glacier has a stable position defined by a compressive arch between lateral pinning points. Ice advance beyond the arch results in calving back to the stable position; conversely, if melt-undercutting forces the ice front behind the stable position, it readvances because ice velocities exceed subaqueous melt rates. This behaviour is typical of self-organising criticality, in which the stable ice-front position acts as an attractor between unstable super-critical and sub-critical regimes. This perspective provides strong support for a ‘position-law’ approach to modelling calving at Sermeq Kujalleq, because any calving ‘rate’ is simply a by-product of how quickly ice is delivered to the critical point. The CD calving law predicts ice-front position from the penetration of surface and basal crevasse fields, and accurately simulates super-critical calving back to the compressive arch and melt-driven calving into the sub-critical zone. The CD calving law reflects the glaciological controls on calving at Sermeq Kujalleq and exhibits considerable skill in simulating its mean position and seasonal fluctuations.

1. Introduction

Recognition of the important role of marine-terminating glaciers in the mass balance of the cryosphere has prompted rapid progress in understanding calving processes and their controls, and considerable advances have been made in both observational and modelling capability (e.g. Benn and Åström, 2018; Fried and others, 2018; Cook and others, 2022). To accurately predict future mass changes, however, this understanding needs to be translated into practical, robust methods for representing calving processes in ice-sheet models. A long-standing problem is the complexity of calving processes: fracture propagation and iceberg detachment can occur by multiple processes with a wide range of possible controls (Benn and others, 2007; Bassis and Jacobs, 2013). This means that simple calving parameterisations tuned for particular glaciers and epochs will likely yield inaccurate results if applied to other locations and time periods (Amaral and others, 2020). To avoid this problem, calving laws ideally should be based on first principles, avoiding the need for tuning on a case-by-case basis.

Calving rate (or frontal ablation rate) can be defined as the difference between vertically averaged ice velocity at the glacier terminus and glacier length change over time. The problem of framing calving laws can therefore be approached in two ways. First, the rate at which calving erodes the ice front can be treated as a function of variables such as water depth or near-terminus strain rates, and then combined with ice velocity to determine changes in terminus position (rate laws); and second, the position of the glacier terminus at any given time can be predicted from glaciological variables such as the height of the ice cliff above buoyancy or crevasse depths (position laws). Good reviews of both types of calving law can be found in Choi and others (2018) and Amaral and others (2020). In a comparison of six calving laws, Amaral and others (2020) found that position laws performed better than rate laws in reproducing observations of 50 marine-terminating glaciers in Greenland. Overall, the crevasse depth (CD) calving law offered the best balance of high accuracy and low sensitivity to parameter calibration.

The CD calving law is based on the idea that ice-front position can be predicted from the depth of penetration of surface and basal crevasse fields, which in turn can be estimated from the state of stress in the glacier (Benn and others, 2007; Nick and others, 2010; Todd and others, 2018). It has been shown to perform well in a range of situations (e.g. Otero and others, 2010; Todd and others, 2019; Cook and others, 2022; Holmes and others, 2023) but a number of issues remain to be resolved, including the most appropriate functions used to estimate crevasse depths, the absence of ‘memory’ in the system and the interpretation of model calibration parameters (Amaral and others, 2020; Enderlin and Bartholomaus, 2020). Furthermore, to date there has been no rigorous demonstration that the physical basis of the CD calving law accurately reflects observed conditions on a tidewater glacier.

In this paper, we investigate the physical basis of the CD calving law by assessing controls on terminus position of Sermeq Kujalleq (Store Gletsjer/Glacier), west Greenland. We characterise calving regimes using a dense time series of Sentinel-1 imagery and uncrewed aerial vehicle (UAV) surveys, then use the full-Stokes continuum model Elmer/Ice and the Helsinki Discrete Element Model (HiDEM) to examine the relationship between glaciological stresses and calving processes. We investigate the degree to which near-terminus stress fields determine the position of the calving front, and consider how the stress fields are modified by ice advance and melt-undercutting. We conclude by evaluating the ability of the CD calving law to determine summer ice-front positions of the glacier.

2. Sermeq Kujalleq (Store Gletsjer/Glacier)

Sermeq Kujalleq is the second largest outlet glacier in West Greenland in terms of ice flux (Weidick and Bennike, 2007). The glacier catchment extends 280 km inland to the ice divide, has a typical width of 50 km narrowing to 5 km at the terminus. Sermeq Kujalleq advances and retreats several hundred metres on annual and sub-annual timescales, although the mean terminus position has remained largely unchanged since at latest 1948 (Weidick, 1995), a period of over 70 years. The mean ice-front position coincides with a narrow and shallow reach of the fjord, which provides a very effective pinning point (Todd and Christoffersen, 2014). Upstream of this pinning point, the glacier occupies a 30 km long overdeepening, reaching a depth of 900 m below sea level. Mean annual velocity is ~5800 m a⁻¹ (16 m d⁻¹) across most of the glacier terminus (Fig. 1), with seasonal fluctuations due to changing fjord conditions and subglacial hydrology (Walter and others, 2012; Cook and others, 2022).

Figure 1. (a) Location of Sermeq Kujalleq (Store Gletsjer/Glacier) in West Greenland; (b) surface velocities on the lower glacier tongue in typical summer conditions; (c) surface velocities under typical spring conditions with frozen ice mélange in the fjord. Background images: (a) Sentinel 1A from 04/09/2018; (b, c) TerraSAR-X images from 29/06/2013 and 29/03/2015.

The absence of long-term climate-forced trends at Sermeq Kujalleq makes it an ideal location to investigate how glaciological factors and fjord conditions control calving processes and terminus position. A growing body of work provides important data on the calving behaviour of the glacier and the characteristics of the adjacent fjord. Dynamic response of the glacier to variations in ice mélange and tidal backpressure has been investigated by Walter and others (2012). The morphology of the submerged glacier front has been surveyed by Rignot and others (2015), showing that the glacier tongue is floating over the deepest part of the trough and is grounded elsewhere. The grounded part of the tongue is undercut by up to at least 50 m as a result of subaqueous melting, which is locally enhanced by upwelling meltwater plumes (Xu and others, 2013; Chauché and others, 2014). Short-term changes in glacier terminus morphology, including records of calving events, have been documented using repeat high-resolution UAV surveys, terrestrial radar interferometry and time-lapse photography (Ryan and others, 2015; Chudley and others, 2019; Cook and others, 2021). Calving at Sermeq Kujalleq has been modelled by Morlighem and others (2016), using a rate law in the ‘shallow shelf’ implementation of the Ice Sheet System Model (ISSM). Most recently, detailed modelling studies using the full-Stokes continuum model Elmer/Ice and the CD calving law have investigated controls on the annual cycle of ice-front position (Todd and others, 2018, 2019) and the influence of subglacial hydrology on the location of meltwater plumes and melt of the submerged ice front (Cook and others, 2020). A fully coupled model of the glacier, including ice flow, subglacial hydrology, frontal melting and calving, has been presented by Cook and others (2022). In this contribution, we use observational data and simulations with Elmer/Ice and HiDEM to investigate the relationships between glacier stress regime, calving and ice-front positions, and elucidate the physical basis for the CD calving law.

3. Methods

3.1 Observations

Glacier-front positions were manually digitised from Sentinel-1A imagery from 11 October 2014 to 17 June 2020 and Sentinel-1B imagery from 10 June 2016 to 23 June 2020. Image spatial resolution is ~20 m. All images were re-projected to the Polar Stereographic coordinate system before the fronts were digitised. Ice-front positions were then recorded using the intersection points of the digitised fronts with a series of linear flowlines with ~150 m spacing at the terminus (cf. Luckman and others, 2015; Fig. 2).

Figure 2. Sampling frame for measuring changes in terminus position, showing the range of observed ice-front positions for 2016 and the extent of the North, Central and South sectors (coloured lines). Also shown are the positions used for the ‘advanced’ (white line) and ‘retreated’ (black line) model domains. The ‘advanced’ and ‘retreated’ ice-front positions differ slightly from those shown in Figures 8 and 9 due to changes that occurred during model relaxation.

To highlight across-glacier variations in ice-front behaviour, mean positions were calculated for three subsets of the ice front (Fig. 2): north (~850 m width), south (~1350 m width) and central (~2500 m width). These sectors were defined following manual identification of recurrent patterns of ice-front evolution on the satellite imagery.

To determine patterns of strain at the glacier surface, ice surface velocities were derived from feature tracking of TerraSAR-X image pairs in slant range using correlation windows of 200 × 200 pixels at every 20 pixels, and subsequently ortho-rectified to a pixel size of 40 m using a DEM (Luckman and others, 2015). Uncertainties in surface velocity are estimated to be 0.4 m d⁻¹ and comprise a co-registration error (~0.2 pixels) and errors arising from unavoidable smoothing of the velocity field over the feature-tracking window.

For 2D strain in the horizontal x y plane, the strain rate components are:

(1a)$$\dot{\varepsilon }_x = \displaystyle{{{\rm d}u} \over {{\rm d}x}}$$

(1b)$$\dot{\varepsilon }_y = \displaystyle{{{\rm d}v} \over {{\rm d}y}}$$

(1c)$$\dot{\varepsilon }_{xy} = \displaystyle{1 \over 2}\left({\displaystyle{{{\rm d}v} \over {{\rm d}x}} + \displaystyle{{{\rm d}u} \over {{\rm d}y}}} \right)$$

where u is the velocity in the x direction, and v is the velocity in the y direction.

The principal strain rates are:

(2)$$\dot{\varepsilon }_{1, 2} = \displaystyle{{{\dot{\varepsilon }}_x + {\dot{\varepsilon }}_y} \over 2} \pm \sqrt {{\left[{\displaystyle{{{\dot{\varepsilon }}_x-{\dot{\varepsilon }}_y} \over 2}} \right]}^2 + {\dot{\varepsilon }}_{xy}^2 } $$

The principal strain rates lie at an angle θ from the coordinate axes:

(3)$${\rm tan}2\theta = \displaystyle{{2{\dot{\varepsilon }}_{xy}} \over {{\dot{\varepsilon }}_x-{\dot{\varepsilon }}_y}}$$

A minimum value of 0.0005 d⁻¹ (one-tenth of the observed maximum value) was adopted for plotting and interpreting strain orientation data.

Data on the magnitude and location of individual calving events were derived from UAV surveys conducted in 2017 (Chudley and others, 2019). Overlapping imagery was captured using a Sony α6000 camera mounted on a Skywalker X8 2 m fixed-wing UAV flown at an altitude of ~450 m above ground level. Forward overlap was 80% and sidelap was 60%, targeting a ground sampling distance of ~11 cm. Multi-View Stereo (SfM-MVS) photogrammetry was used to produce 3D models using Agisoft Metashape software. Models were geolocated via aerial triangulation using an L1 carrier-phase GPS receiver mounted on the UAV, post-processed kinematically against a bedrock-mounted GPS base station. For a full description of the methods, see Chudley and others (2019).

3.2 Modelling

Model experiments were conducted using the continuum model Elmer/Ice and the discrete element model HiDEM. Elmer/Ice was used to (1) determine values of basal shear stress from observed velocity data; (2) analyse the stress regimes associated with a range of glacier configurations; and (3) predict changes in ice-front position using the CD calving law. HiDEM was used to simulate patterns of fracture and calving for selected glacier configurations imported from Elmer/Ice.

Elmer/Ice is an open-source, finite-element model that solves the unaltered Stokes equations (aka ‘full-Stokes’), treats grounding-line dynamics as a contact problem, and provides inverse methods for deriving basal and englacial conditions (Gagliardini and others, 2013). Model code is freely available (https://github.com/ElmerCSC/elmerfem), and detailed descriptions of parts of the ice flow, calving and remeshing algorithms are provided in Todd and others (2018). HiDEM is described in detail by Åström and others (2013) and van Dongen and others (2020), and is freely available at https://github.com/joeatodd/HiDEM. HiDEM represents ice as particles connected by breakable elastic beams, stacked together to form specified 3D glacier geometries. Although individual particles are rigid, in bulk the material is compressible and elastic. The model is initialised with specified cumulative strain thresholds for beam breakage and density of randomly scattered small pre-existing cracks. Values of breakage threshold = 0.0003 and damage density = 0.1 were chosen to yield bulk ice properties consistent with observations (Young's modulus ~1 GPa, Poisson's ratio ~0.2). At the beginning of a simulation, the domain was allowed to relax without fracture, creating an initial static elastic stress field. Thereafter, if cumulative strain on beams reaches the failure threshold, the beam breaks and particles (or aggregates) become disconnected but continue to interact as long as they are in contact. In this way, cracks in the modelled ice body are formed and propagated and may result in the detachment and release of calved blocks.

For the simulations reported here, the model domains extend 12 km (Elmer/Ice) and 5 km (HiDEM) upstream from the terminus. To represent the range of observed terminus positions, we used two configurations of the glacier terminus as initial geometries for model simulations: (1) a ‘retreated’ configuration based on data from 8 July 2016; and (2) an ‘advanced’ configuration using data from 26 April 2016. These configurations do not bracket the entire range of observed terminus positions, but provide alternative starting points from which the ice may advance or retreat during simulations. The temperature structure for the model domain was obtained from the simulations of Todd and others (2018), which solved the thermo-mechanical problem over the whole catchment.

Glacier surface DEMs were derived from ArcticDEM products (nominal date July 2015). The DEM was geoid corrected, low-pass filtered at 90 m to remove surface crevasses and extrapolated where necessary. The final product was up-sampled to 30 m. Surface mass balance was from RACMO 2.3 (mean of 1958–2013 data). We generated a glacier bed DEM using a mass-conservation approach, constrained by thickness data from Operation IceBridge flight lines (https://espo.nasa.gov/missions/oib/) (Todd and Christoffersen, 2014), downsampled to 10 m to account for the finer mesh resolution in this study. We generated our own DEM in preference to the available BedMachine product (Morlighem and others, 2017), which exhibits a number of probable artefacts in the terminal zone of Sermeq Kujalleq. Values of basal shear stress were obtained by using Elmer/Ice to assimilate MEaSUREs velocity data (https://nsidc.org/data/nsidc-0478) for August 2016, which were closest to the mean of the full range of available data (September 2000 to August 2016).

We use an unstructured, vertically extruded mesh in Elmer/Ice, with a horizontal edge length varying from 45 m at the terminus to 300 m at the inflow boundary. The mesh is vertically extruded to 15 terrain-following layers. We begin by initialising the model in Elmer/Ice, use observed velocities to invert for basal shear stress and allow the domain to relax slightly (0.1 years) to avoid issues arising from imprecision in the input data. Inversions were run only on the ‘advanced’ domain, and the resulting slip coefficient map was used for both domains in forward simulations. The optimisation routines could only reproduce the observed centreline velocities by softening the ice in lateral shear margins. We therefore imposed a Glen Enhancement Factor of up to 6 in imposed marginal shear sectors 100 m wide, the positions of which were guided by observations (cf. Cuffey and Paterson, 2010). The enhancement factor drops off as a Gaussian curve from the centre of the shear zone to a background value of 1.0. Relaxed Elmer/Ice geometries were used to produce HiDEM domains with hexagonal close-packed (hcp) lattices and particle size of 30 m.

Calving in Elmer/Ice was implemented using the CD calving law (Benn and others, 2007; Nick and others, 2010), modified to include all stress components from the solution to the Stokes equations (Todd and others, 2018). Fracturing is assumed to be possible wherever the largest principal stress σ ₁ (the largest eigenvector of the Cauchy stress) is extensional (positive). Water pressure can promote fracture propagation by opposing ice overburden pressure. The effect of water on the depth of surface crevasses is likely highly localised (Chudley and others, 2021) and is not included here as a control on crevasse depths. In contrast, all basal crevasses are assumed to be accessed by pressurised water from the subglacial hydrological system. Near the terminus, subglacial water pressure P _W is assumed to be equal to the pressure exerted by the sea on the glacier front, which acts as hydrological base level. The fracture criteria can therefore be expressed as:

(4a: surface crevasses)$$\sigma _1 > 0$$

(4b: basal crevasses)$$\sigma _1 + P_{\rm W} > 0$$

It is important to note that these expressions are not intended to predict the depths of individual crevasses, which are strongly influenced by local stress intensity factors (cf. Jimenez and Duddu, 2018). Rather, Eqns (4a) and (4b) are used to approximate the extent of regions within the glacier where ice is assumed to be damaged to a greater or lesser degree, and are employed as a consistent and efficient method for predicting the extent of damage for every time step across the entire model domain, without the need for a priori knowledge of crevasse spacing and other factors. In this study, we follow previous convention and predict calving events where the region of surface damage meets the waterline, or where the regions of surface and basal damage meet and isolate portions of the ice front (cf. Todd and others, 2018, 2019; Cook and others, 2022).

The form of the CD calving law used here is based on the principal Cauchy stresses. For illustrating the relationships between stress and fracture patterns, however, principal deviatoric stresses (i.e. the difference between the Cauchy eigenvalues and the mean or hydrostatic stress) offer greater clarity. The deviatoric stresses are that part of the total stress that can cause irrecoverable strain, and hence the possibility of fracture. In Figures 8–12, the first and second principal deviatoric stresses $\sigma _1^{\rm ^{\prime}} $ and $\sigma _2^{\rm ^{\prime}} $ are visualised at the glacier surface and at 50 and 90% of the glacier thickness.

4. Ice-front oscillations and calving processes

Throughout the 2014–2020 period covered by the Sentinel-1 data, the terminus of Sermeq Kujalleq underwent seasonal and shorter-term oscillations about a stable mean position (Fig. 3). The terminus position at the lateral margins remained almost constant, while the greatest variation occurred in the central part of the glacier. Seasonal cycles are most clearly evident in the North and South sectors, with net advance of ~400 m during the autumn, winter and spring, followed by rapid retreat over a similar distance in June and July. In the Central sector, seasonal cycles are weak and inconsistent, and large amplitude (~600 m) oscillations record repeated large calving events and readvances throughout the year. Some, but not all of the oscillations in the Central sector correspond to events in the South and (to a lesser extent) North sectors.

Figure 3. Time series of ice-front positions for the north, central and south sectors of Sermeq Kujalleq, determined from Sentinel-1 imagery. The absolute positions of each series are offset for clarity.

The advance-retreat cycles that occur in the summer months exhibit characteristic recurrent spatial patterns. A persistent embayment develops in the North sector, while the South sector undergoes repeated minor advances and calves back to a consistent position with a linear ice front. In the Central sector, irregular promontories form by ice flow followed by episodic detachment of icebergs, including large, tabular bergs. Figure 4 shows TerraSAR-X images from 28th June, 9th July and 20th July 2014. During the period spanned by these images, a large section of the glacier front was lost, extending from the south side of the glacier to a promontory near the boundary between the Central and North sectors. The promontory and part of the southern ice front calved first, the latter following the line of a rift visible on the first image. A large portion of the Central sector calved next, creating a new promontory which was advected forward by ice flow. It is not known whether the calving occurred in just two events or via a number of smaller events between the dates of the imagery.

Figure 4. TerraSAR-X images of Sermeq Kujalleq on (a) 28^th June, (b) 9^th July and (c) 20^th July 2014. The red line in panels (b) and (c) indicates the ice-front position in (a).

Insight into smaller-scale calving processes is provided by imagery derived from UAV surveys (cf. Ryan and others, 2015; Chudley and others, 2019; Cook and others, 2021). An example from the South sector is shown in Figures 5a and b. In the image for 12 July 2017, open surface crevasses can be observed corresponding with the post-calving ice front the following day, indicating that failure was preceded by longitudinal extension at the glacier surface, consistent with forward rotation or extension of the glacier front. A second example from the North sector of the glacier is shown in Figures 5c and d. The two images were acquired 6.5 h apart on 11 July 2017, and show calving along a line of surface crevasses, consistent with fracture propagation associated with forward rotation of the ice front. A similar calving event at the same location was described by Ryan and others (2015) based on UAV surveys on 1 and 2 July 2013 (their calving event A).

Figure 5. (a, b) Map-view orthoimages produced from UAV surveys on 12 and 13 July 2017, showing a calving event in the southern part of the glacier front, along the line of a surface crevasse. (c, d) Map-view orthoimages from UAV surveys on 11 July 2017 at 10.20 (a) and 16.50 UTC (b), showing a calving event in the north sector of the glacier. Note surface crevasses along and close to the line of calving failure (arrowed).

These examples, and other examples described by Cook and others (2021), are consistent with calving triggered by melt-undercutting, in which forward rotation of the upper part of the glacier front occurs in response to increased torque following loss of ice below the waterline (cf. van Dongen and others, 2020; Slater and others, 2021). Analysis of radar interferometry data and time-lapse imagery of the North sector by Cook and others (2021) shows a bimodal distribution of calving event size, with smaller events (~10⁴ m³) representing the detachment of portions of the subaerial ice front and larger events (~10⁵ m³) representing forward rotation and deep propagation of surface crevasses. No very large tabular-style calving events were observed in their data series. Cook and others (2021) also found a sudden increase in the number and volume of calving events coincident with ice mélange break-up, and a greater number (but not total volume) of calving events in the vicinity of active meltwater plumes.

5. Basal stress and patterns of strain

To investigate the physical basis of the CD calving law, we begin by investigating relationships between large-scale patterns of stress and strain and the position of the ice front. Figure 6 shows the subglacial topography, extent of floating ice and basal shear stress τ _B below the glacier tongue. Extent of flotation and basal shear stress were determined in Elmer/Ice, using surface and bed DEMs and observed velocity data as inputs. Over large parts of the lower tongue of the glacier, the basal shear stress is <10% of the local driving stress τ _D (Fig. 6b). Conversely, there are regions along both margins (outlined in blue) where τ _B > τ _D, indicating that the lateral margins play a large role in supporting the driving stress. There is zero basal drag below the floating southern portion of the tongue, whereas around 5–10% of resistance to the driving stress is provided by a bedrock ridge below the northern side of the glacier terminus.

Figure 6. (a) Basal topography below the lower tongue of Sermeq Kujalleq, with the extent of floating ice shown by the white line. (b) Basal stress τ _B, expressed as a proportion of the driving stress τ _D. The regions outlined in blue indicate where τ _B > τ _D. Note logarithmic scale.

Patterns of strain at the glacier surface, determined from velocity data using Eqns 1–3, are shown in Figure 7. In ‘summer’ conditions, when frozen ice mélange is absent, the first principal strain rate $\dot{\varepsilon }_1$ is positive (extensional) almost everywhere due to acceleration and transverse spreading of the ice as it approaches the terminus, and strong lateral shear. The second principal strain rate $\dot{\varepsilon }_2$ is negative (compressional) over much of the glacier tongue (Fig. 7b). The compressional strain vectors are typically aligned obliquely up-glacier from the margins to the centreline, forming a chevron-like pattern on which local variations are superimposed. Near the glacier terminus, the pattern of strain changes. Orientations of both $\dot{\varepsilon }_1$ and $\dot{\varepsilon }_2$ are variably longitudinal, transverse and oblique, and $\dot{\varepsilon }_2$ is extensional or weakly compressive. The approximate up-glacier limit of this region is highlighted in yellow in Figures 7a and b. These patterns can be seen more clearly in the modelled first and second principal deviatoric stresses ($\sigma _1^{\rm ^{\prime}} \;$ and $\sigma _2^{\rm ^{\prime}} $) for the ‘retreated’ glacier geometry (Fig. 8), which is similar to the configuration shown in Figure 7. In particular, the $\sigma _2^{\rm ^{\prime}} $ vectors display a clear concave-downglacier arcuate pattern, which breaks down near the terminus. This arcuate pattern of mainly compressive stress (and strain) is interpreted as a compressive arch. Within the arch, drag at the lateral margins is transmitted by membrane stresses in the ice, allowing much of the driving stress to be supported by non-local boundary stresses (cf. Doake and others, 1998). Downglacier of the arch, loss of lateral support results in complex patterns of stress and strain including large areas where $\sigma _1^{\rm ^{\prime}} \;$and $\sigma _2^{\rm ^{\prime}} $ (and $\dot{\varepsilon }_1$ and $\dot{\varepsilon }_2$) are positive.

Figure 7. First and second principal strain rates at the glacier surface under (a, b) ‘summer’ conditions (no mélange); and (c, d) ‘winter’ conditions with frozen mélange. Colours represent strain rate magnitudes and the white lines indicate vector orientations. The bold yellow line in panels (a) and (b) shows the approximate downglacier limit of ice supported by lateral drag.

Figure 8. First $( \sigma _1^{\rm ^{\prime}} ) $ and second $( \sigma _2^{\rm ^{\prime}} ) $ principal deviatoric stress magnitudes (colours) and orientations (green lines) for the ‘retreated’ configuration.

Patterns of strain during March 2015, when frozen ice mélange was present in the fjord, are shown in Figures 7c and d. The first principal strain rate remains positive (extensional) over most of the lower tongue, but the second principal strain rate displays a very different pattern to the summer state. The areas of extensional strain downglacier of the compressive arch have largely disappeared, although the arch itself can still be discerned. The main effect of the ice mélange, therefore, appears to be to support the ice downstream of the compressive arch, suppressing extensional strain and potentially inhibiting calving.

6. Model results: stress, fracture and calving

6.1 The ‘retreated’ and ‘advanced’ geometries

The first set of model experiments uses Elmer/Ice and HiDEM to explore the patterns of stress, fracture and calving associated with the ‘retreated’ and ‘advanced’ model geometries. For the ‘retreated’ geometry, the first principal deviatoric stress $\sigma _1^{\rm ^{\prime}} \;$at the surface is positive (extensional) almost everywhere except in the southern, floating part of the terminus where there is an oval region of compressive stress at the surface immediately downglacier of the grounding line (Fig. 9a). At 50 and 90% depth, the first principal deviatoric stress in this region becomes extensional (Figs 9c, e). This pattern of compressive stress at the surface and extensional stress at depth is consistent with a ‘bottom-out’ torque, interpreted as the result of buoyant forces acting on the ice as it adjusts towards hydrostatic equilibrium downglacier of the grounding line (cf. James and others, 2014; Murray and others, 2015; Benn and others, 2017). In HiDEM, several transverse basal crevasses are initiated in this area (blue lines in Fig. 9), but do not propagate the full thickness of the ice.

Figure 9. First $( \sigma _1^{\rm ^{\prime}} ) $ and second $( \sigma _2^{\rm ^{\prime}} ) $ principal deviatoric stress magnitudes (colours) and fracture patterns (lines) for the ‘retreated’ configuration. Stress magnitudes were computed in Elmer/Ice and fractures simulated using the same geometry in HiDEM. Stresses are shown for the glacier surface (a, b), 50% depth (c, d) and 90% depth (e, f). Major fractures are shown in red (surface), blue (basal), black (full-depth) and green (internal fractures that do not intersect the surface or the bed). Locations of fracture initiation are indicated by stars. Areas of ungrounded ice are delineated with grey dashed lines.

The second principal deviatoric stress is negative (compressional) over much of the glacier tongue for the ‘retreated’ geometry, as is the case for the strain data (Fig. 9b). The areas of positive (tensile) $\sigma _2^{\rm ^{\prime}} $ display similar patterns to those in the strain data (Fig. 7b), but are resolved in greater detail. At all depths, areas of tensile $\sigma _2^{\rm ^{\prime}} $ occur up-glacier of the grounding line (where the ice passes over a bedrock bump) and close to the glacier terminus. The near-terminus zone of positive $\sigma _2^{\rm ^{\prime}} $ is consistent with loss of lateral support from the glacier margins as ice flows past the compressive arch at the pinning point. At the glacier surface, this pattern is interrupted by a region of compressive $\sigma _2^{\rm ^{\prime}} $ in the northern half of the glacier terminus, possibly associated with ice flow against the subglacial bedrock ridge.

In the HiDEM simulation with the ‘retreated geometry’, a series of arcuate crevasses formed across the South and Central sectors of the glacier terminus. Most are surface crevasses, and only in a few places penetrated the full thickness of the glacier, and no calving had occurred anywhere in the domain by the end of the simulation. Thus, although the patterns of $\sigma _2^{\rm ^{\prime}} $ indicate an increase in tensile membrane stress as ice flows past the pinning point, this is insufficient to initiate calving, suggesting that the ‘retreated’ configuration of the glacier is stable.

For the ‘advanced’ geometry, the patterns of first and second principal deviatoric stress are similar to the ‘retreated case’, except for the occurrence of a much larger region near the terminus where both $\sigma _1^{\rm ^{\prime}} $ and $\sigma _2^{\rm ^{\prime}} \;$ are mainly positive (tensile) (Fig. 10). This shows that, in the absence of backstress from ice mélange in the fjord, ice in the ‘advanced’ position is subject to tensile membrane stresses in both horizontal dimensions due to lack of support from the glacier margins and, in the northern part of the glacier, the subglacial bedrock ridge. The HiDEM simulation for the ‘advanced’ geometry produced full-depth fractures across almost the entire glacier, triggering widespread calving and ice retreat back to a position similar to the ‘retreated’ configuration.

Figure 10. As for Figure 9, but for the ‘advanced’ glacier geometry. The cross-hatched regions represent ice that calved in the HiDEM simulations.

In summary, the experiments with the ‘retreated’ and ‘advanced’ configurations highlight the importance of the near-terminus force balance as a control on ice fracture and calving. For the ‘retreated’ configuration, the effect of ice flow past the compressive arch is evident in the near-terminus region where both $\sigma _1^{\rm ^{\prime}} $ and $\sigma _2^{\rm ^{\prime}} \;$ are mostly tensile at all depths. In this region, lateral support is lost and basal drag is small or zero (Fig. 6b), but tensile membrane stresses are not sufficiently high to trigger widespread full-depth fracture and calving. For the ‘advanced’ configuration, however, both $\sigma _1^{\rm ^{\prime}} $ and $\sigma _2^{\rm ^{\prime}} \;$ are tensile over a much larger area, and full-depth fracture and calving are widespread. In the absence of backstress from frozen ice mélange in the fjord, therefore, the ‘advanced’ configuration appears to be unstable.

6.2 Ice advance from the ‘retreated’ position

A second set of experiments was conducted to determine the effects of ice advance from the stable ‘retreated’ position. Elmer/Ice was used to simulate ice advance, with the calving function switched off, then ice geometries at timesteps ranging from 2 to 18 d of ice advance were imported into HiDEM to determine patterns of ice fracture and calving. No calving occurred after 2 d of advance, but after 4 d (i.e. ~60 m of ice advance) a calving event in the South and part of the Central sector caused ice retreat back to near the initial position. As the period of ice advance was increased, the lateral extent of calving extended progressively northward, but even after 18 d it did not extend into the North sector.

Results for 8 and 16 d of ice advance are shown in Figures 11 and 12. Both cases exhibit similar stress patterns. In the southern parts of the glacier front, $\sigma _1^{\rm ^{\prime}} $ remains tensile and increases in magnitude, while regions of tensile $\sigma _2^{\rm ^{\prime}} \;$ increase in area and magnitude compared with the initial ‘retreated’ position. These changes are greater, and extend farther north, after 16 d of advance than for 8 d of advance, and it is notable that the major fractures and failure surfaces simulated in HiDEM closely follow areas where extensional principal deviatoric stresses have increased (Figs 11c, d, 12c, d).

Figure 11. (a, b) Magnitude and orientation of first and second principal deviatoric stresses for 8 d of ice advance from the ‘retreated’ configuration. (c, d) Difference in principal deviatoric stresses between the 8 d advance case and the initial ‘retreated’ configuration, together with fractures modelled in HiDEM. Hatched areas indicate ice that calved during the simulation.

Figure 12. (a, b) Magnitude and orientation of principal deviatoric stresses for 16 d of ice advance from the ‘retreated’ configuration. (c, d) Difference in principal deviatoric stresses between the 16 d advance case and the ‘retreated’ case, together with fractures modelled in HiDEM. Hatched areas indicate ice that calved during the simulation.

In the northern part of the terminus, $\sigma _1^{\rm ^{\prime}} $ is tensile along-flow, and the tensile stress becomes larger (more positive). However, $\sigma _2^{\rm ^{\prime}} \;$ remains compressional (though less negative) reflecting the influence of the bedrock ridge, and calving does not occur.

6.3 Melt undercutting at the ‘retreated’ position

A third set of experiments was conducted to investigate the effect of undercutting on the stable ‘retreated’ configuration. The experiments shown in Figure 13 illustrate the impact of distributed undercutting, with undercuts increasing linearly from zero at the waterline to maxima at the bed. The figure shows the difference in $\sigma _1^{\rm ^{\prime}} $ compared with the non-undercut case, at the glacier surface and 90% depth for undercuts of 60, 100 and 200 m. These undercuts are equivalent to 12, 20 and 40 d of melting, respectively, for typical melt rates of 5 m d⁻¹ (Cook and others, 2020).

Figure 13. Difference in first principal deviatoric stress for undercut ‘retreated’ configurations compared with the non-undercut case, and fractures modelled in HiDEM. Stresses shown for the glacier surface (a, b, c) and at 90% depth (d, e, f), with linear undercuts of 60 m (a, d), 100 m (b, e) and 200 m (c, f). A visual impression of undercut size is given by the area of ‘missing’ ice in the panels for 90% depth (uncoloured stippled areas). Key for fractures and calving same as in Figure 9.

In all cases, undercutting increases tensile first principal deviatoric stresses at the glacier surface, with the effect becoming stronger for deeper undercuts. At 90% depth the absolute values of $\sigma _1^{\rm ^{\prime}} $ (not shown) remain tensile but with decreasing magnitude (i.e. less tensile) as undercutting increases. These patterns reflect ‘top-out’ torque in response to asymmetric forces acting on the front (cf. Benn and others, 2017; Slater and others, 2021). For relatively small undercuts, calving occurs only in the South sector of the glacier, and with increasing undercutting the calved area expands northwards and affects the entire front for the deepest undercut.

Experiments were also run for undercuts with conical geometries to represent the effect of localised undercutting associated with upwelling meltwater plumes at known locations. However, these had very little effect on calving, and the results are not shown here.

6.4 Simulated and observed calving processes

The styles of calving in the HiDEM simulations are similar in many respects to those described in Section 4 above. In all of the ice-advance simulations, the ice calves back to positions similar to the ‘retreated’ front geometry, consistent with the repeated advance-retreat cycles shown in Figure 3. The simulations, together with the Elmer/Ice visualisations of the associated stress fields, support the idea that ice advance beyond the compressive arch and loss of basal and lateral drag results in increased membrane stresses in the unsupported ice, which in turn promotes fracture propagation and calving. The repeated advance-retreat cycles seen in the Sentinel-1 data are thus interpreted as a manifestation of quasi-periodic cycling of stress build-up and failure. It is notable that in the HiDEM simulations ice advance beyond the ‘retreated’ position for 4 d or more results in calving; this is similar to observed calving front oscillations in the central portion of the glacier in summer and autumn, in which calving occurs after only minor advances (~100 m). Greater advances of the front tend to occur only in the winter and spring (January to June), when additional support is provided by ice mélange. The presence of ice mélange, however, is insufficient to completely suppress calving in this sector.

Calving events in the undercut simulations are similar to those imaged in the UAV data. Figure 14 shows a snapshot of a calving event in the HiDEM run with a linear undercut of 60 m (this is the same run shown in the left panels in Fig. 13), which is strikingly similar in both location and magnitude to the observed calving event shown in Figures 5a and b. The presence of open surface crevasses in the UAV image prior to the observed calving event indicates that failure was preceded by longitudinal extension of the glacier surface, consistent with forward rotation of the glacier front in response to melt-undercutting (cf. Benn and others, 2017; van Dongen and others, 2020).

Figure 14. (a, b) Calving event in the south sector of the glacier (from Fig. 5); (c) snapshot of HiDEM simulation for the 60 m undercut case, showing close similarity in both location and magnitude to the observed event.

7. Calving predicted by the CD calving law

The model experiments described above have demonstrated that the position of the ice front is located where stresses are sufficient to propagate fractures through the ice, as a result of either ice advance beyond the pinning point or melt-undercutting. This is consistent with the idea that underpins the CD calving law: that glaciological stresses determine patterns of fracture, which in turn determine where calving will occur. We now compare the HiDEM model results with ice-front positions predicted by the CD calving law implemented in Elmer/Ice, to determine whether the two models converge on similar solutions and to consider the reasons for any similarities and differences. Figure 15 shows the results for two of the ice-advance experiments and two of the undercut experiments. For the ice-advance cases, the calving margin predicted by Elmer/Ice is less regular than the HiDEM prediction, with more promontories and embayments. This is especially the case for the 16 d advance experiment, in which Elmer/Ice predicts a large embayment in the floating part of the glacier tongue. The two models show better agreement in the melt-undercutting experiments, although Elmer/Ice predicts that promontories remain in the floating region, whereas HiDEM predicts a more linear ice-margin configuration.

Figure 15. Simulated calving front positions using HiDEM (purple line) and the crevasse-depth calving law in Elmer/Ice (yellow line). (a) After 8 d of advance from ‘retreated’ position; (b) after 16 d of advance; (c) ‘retreated’ position with 100 m undercut; (d) ‘retreated’ position with 200 m undercut. The background colours indicate the ratio of crevassed to uncrevassed ice in the ice column calculated using the CD calving law in Elmer/Ice.

The differences in the results reflect the contrasting structure of the two models. The calving law in Elmer/Ice is based on theoretical extent of crevassing calculated from the stress field in undamaged ice, and thus neglects prior damage and feedbacks between fracture propagation and stress evolution. In contrast, the HiDEM simulation computes evolving inter-particle stresses as bonds stretch and break, and thus includes feedbacks between fracture propagation and the state of the surrounding ice. This allows HiDEM to represent processes to which the Elmer/Ice model is insensitive, such as the ability of fractures to propagate across initially intact but vulnerable regions (e.g. the major promontory in Fig. 15b), or the ability of surrounding areas to support and prevent calving of fractured ice (e.g. the large embayment in Fig. 15b). It is notable, however, that in the HiDEM simulation full-depth crevasses extend around much of the area occupied by the embayment (Fig. 13), although not enough to trigger calving.

Despite differences in detail, however, there is generally close agreement between the two models, and the predicted ice fronts typically differ by only a few tens of metres. These differences were quantified by measuring the difference between the Elmer/Ice and HiDEM calving positions along the 33 longitudinal lines shown in Figure 2. The overall mean difference is +37 m (i.e. the HiDEM position is on average 37 m more advanced than the Elmer/Ice position), with a std dev. of 134 m. Of all the possible positions that the two models might have predicted (from no calving at all to complete calving of the entire terminal zone) and the radical differences in model formulation, this result indicates a high degree of agreement regarding the predicted location of the calving front. Furthermore, this location is closely similar to the ‘retreated’ configuration of the glacier, to which it repeatedly returns over the course of each summer (Fig. 3).

8. Discussion

The model results indicate that the unperturbed ‘retreated’ model geometry (which approximates the mean frontal position of Sermeq Kujalleq in the summer months; Fig. 3) represents a stable configuration of the glacier. This position reflects the location of a compressive arch at the narrowing of the fjord, and a bedrock ridge beneath the North sector of the glacier. Calving occurs when either one of two perturbations is imposed: (1) advance of the ice beyond the stable position, or (2) undercutting of the ice front by subaqueous melting. In the ice-advance case, calving reflects loss of the stabilising influence of lateral drag where the fjord widens (South and Central sectors) and the additional loss of basal drag from the bedrock ridge (North sector). Loss of resistance at the lateral and basal boundaries causes membrane stresses to increase, with both $\sigma _1^{\rm ^{\prime}} $ and $\sigma _2^{\rm ^{\prime}} $ becoming tensile over increasingly large areas. When frozen ice mélange is present in the fjord (typically January until June), the glacier is able to advance beyond the ‘retreated’ position, although it remains susceptible to large calving events (Fig. 3). The stabilising influence of ice mélange was not modelled in this study, but was investigated in the model experiments of Todd and others (2018, 2019) and Cook and others (2022).

Melt-undercutting perturbs the near-terminus stress field by amplifying the top-out torque at the glacier front. This is associated with an increase in tensile stresses at the surface and a decrease in tensile stresses at the base, encouraging the propagation of surface crevasses and suppressing basal crevassing (e.g. O'Leary and Christoffersen, 2013; van Dongen and others, 2020; Slater and others, 2021). Observations in the field (Fig. 5; Chudley and others, 2019; Cook and others, 2021) indicate that ‘top-forward’ calving frequently follows growth of surface crevasses at Sermeq Kujalleq, although small-scale failures of the subaerial part of the ice cliff are more common (but lesser in total volume). Seasonal mean melt rates are below 5 m d⁻¹, compared with ice velocities up to 16 m d⁻¹ at the glacier centreline (Cook and others, 2022). Thus, any ice-front retreat driven by melt-undercutting can be rapidly reversed by ice advection. Previous modelling studies have shown that very high subaqueous melt rates would be required to push the ice front permanently into the deeper water behind the bedrock ridge (Morlighem and others, 2016; Todd and others, 2019).

The fluctuations of Sermeq Kujalleq around a stable position suggest that the system exhibits self-organised criticality (cf. Åström and others, 2014; Chapuis and Tetzlaff, 2014; Cook and others, 2021). The classic example of a self-organised critical system is a sandpile to which grains are added slowly from above (Bak and others, 1987). As grains accumulate, the pile steepens until a critical slope angle is reached, whereupon avalanches transfer grains downslope. In response to the opposing processes of grain accumulation and avalanching, therefore, the system spontaneously evolves towards a stable ‘critical point’ between unstable sub-critical and super-critical regimes (Bak and Paczuski, 1995). In the case of Sermeq Kujalleq, ice advance moves the glacier into a super-critical state, whereupon large-scale calving shifts the ice front back to or beyond the critical point. Conversely, melt-driven calving events can cause glacier retreat into a sub-critical state, whereupon ice advance returns it to or past the critical point. The pinning point thus acts as an attractor, which in this case has maintained stability of the ice-front position for over 70 years (Weidick, 1995). Systems exhibiting self-organised criticality typically have a characteristic power-law distribution of event sizes. Reanalysis of the calving event-size data of Cook and others (2021) by Christoffersen and others (in prep.) shows that, for calving events smaller than 10⁵ m³, the data have a power-law distribution with an exponent of −1.2, the same as that found by Åström and others (2014) for glaciers in Svalbard, Alaska, Greenland and Antarctica. Above that size, calving events exhibit log-normal or exponential distributions, consistent with super-critical states (cf. Åström and others, 2021). Thus, both large-scale system behaviour and the statistical properties of calving events are indicative of self-organised criticality.

Self-organisation of Sermeq Kujalleq's calving front at the pinning point provides strong support for a position law approach to modelling calving. The calving rate at Sermeq Kujalleq is a secondary property governed by rate of ice delivery to the critical position; in such cases, it makes little sense to attempt to use a calving ‘rate’ to predict ice-front position. Calving position laws, such as the CD calving law used in this study, are thus conceptually in much closer alignment with how calving actually works at Sermeq Kujalleq. Moreover, the CD calving law implemented in the full-Stokes model Elmer/Ice exhibits considerable skill in predicting both the ‘retreated’ position of Sermeq Kujalleq and the ‘advanced’ position permitted by the presence of frozen ice mélange (Todd and others, 2018, 2019).

The role of the compressive arch in determining the critical calving front position is related to the concept of first-order calving introduced by Benn and others (2007). First-order calving was originally defined in terms of longitudinal extension in one dimension, but it is more usefully viewed as a critical stability threshold in three dimensions. Compressive arches have long been recognised as key to understanding ice-shelf stability (Doake and others, 1998; Levermann and others, 2012) and, with suitable modification, the concept offers a powerful organising principle for understanding the seasonal and longer-term behaviour of Greenlandic tidewater glaciers (e.g. Cowton and others, 2019). Melt-undercutting (designated a second-order process by Benn and others, 2007) influences the transient ice-front position at Sermeq Kujalleq, but due to the large ice flux, melt rates are insufficient to force the ice front up-glacier of the pinning point on a permanent basis (Cook and others, 2022). At glaciers where ice flux is small, melt-undercutting can be the rate-limiting calving process (e.g. How and others, 2019). In such cases, calving rate laws based on correlations between water temperature and frontal ablation may provide simple and accurate means of predicting ice-front evolution (Luckman and others, 2015).

Finally, this study provides a useful perspective on the question of what a calving law is expected to do. At the event scale, calving is a stochastic process, and the magnitude and timing of individual events may be poorly correlated with environmental variables (Cook and others, 2021). Although the CD law treats calving as a succession of discrete events, its power lies in its ability to determine the limiting position of the ice front where sub-critical behaviour gives way to super-criticality. To provide reliable predictions of tidewater glacier behaviour, a calving law should be capable of determining the location of critical points from glaciological variables (thickness, ice flux, stress, etc.), and predicting the timing and rate of transition to new states (if such exists). In many situations, it will also be desirable to simulate seasonal fluctuations, to determine the magnitude of perturbations to which the ice-front may be subject. In combination with the previous works of Todd and others (2018, 2019) and Cook and others (2022), this paper has shown that, when implemented in a 3D, full-Stokes flow model, the CD law accurately predicts the stable, ‘summer’ ice-front position of Sermeq Kujalleq and the observed seasonal fluctuations of the glacier.

9. Summary

1. (1) The long-term stability of Sermeq Kujalleq's mean ice-front position reflects a major threshold in the stress regime of the glacier at a pinning point. Up-glacier of the pinning point, much of the driving stress is supported by lateral drag transmitted via a compressive arch, whereas downglacier of the pinning point, tensile membrane stresses increase. Stress analysis in Elmer/Ice and calving simulations in HiDEM show that fracture initiation and propagation occur in regions of high tensile stress when ice advances beyond the compressive arch.

2. (2) The location of the mean summer position of the ice front is accurately predicted by both the HiDEM and the CD calving law implemented in Elmer/Ice. Although the two models have completely different formulations, both account for the state of stress in the ice down-flow of the compressive arch and converge on similar solutions.

3. (3) The mean summer ice-front position of Sermeq Kujalleq represents an attractor or critical point between unstable regimes. If ice advances forward from the critical point it moves into the super-critical regime and rapidly calves back to, or beyond the critical point. Conversely, if the ice-front retreats up-glacier of the critical point by melt-driven calving, high ice velocities cause the glacier to readvance. This oscillation around a critical point is typical of a self-organising critical system. The size-frequency distribution of calving events follows a power-law form across a range of scales, indicating that the ice front of Sermeq Kujalleq exhibits self-organising criticality.

4. (4) This study provides strong support for the idea underpinning the CD calving law, that ice-front position can be predicted from the penetration of surface and basal crevasse fields, which in turn are controlled by glaciological stresses. It also demonstrates the ability of the CD law to predict observed locations of the terminus of Sermeq Kujalleq, when implemented in the full-Stokes 3-D model Elmer/Ice. Future work will explore the applicability of the model in other settings.