Fragmentation theory reveals processes controlling iceberg size distributions

Abstract

Iceberg calving strongly controls glacier mass loss, but the fracture processes leading to iceberg formation are poorly understood due to the stochastic nature of calving. The size distributions of icebergs produced during the calving process can yield information on the processes driving calving and also affect the timing, magnitude, and spatial distribution of ocean fresh water fluxes near glaciers and ice sheets. In this study, we apply fragmentation theory to describe key calving behaviours, based on observational and modelling data from Greenland and Antarctica. In both regions, iceberg calving is dominated by elastic-brittle fracture processes, where distributions contain both exponential and power law components describing large-scale uncorrelated fracture and correlated branching fracture, respectively. Other size distributions can also be observed. For Antarctic icebergs, distributions change from elastic-brittle type during ‘stable’ calving to one dominated by grinding or crushing during ice shelf disintegration events. In Greenland, we find that iceberg fragment size distributions evolve from an initial elastic-brittle type distribution near the calving front, into a steeper grinding/crushing-type power law along-fjord. These results provide an entirely new framework for understanding controls on iceberg calving and how calving may react to climate forcing.

1. Introduction

Iceberg production is a key mass loss process for glaciers and ice sheets around the world (Benn and others, 2007). Changes in iceberg calving influence the glacier stress balance and must be accurately incorporated into numerical ice flow models in order to properly constrain mass loss in a changing climate (Pollard and others, 2015; Choi and others, 2018). One complicating factor is that iceberg calving can occur as the result of several mechanisms (Benn and others, 2007), from the release of large tabular icebergs from floating ice shelves (Joughin and MacAyeal, 2005), to the toppling of smaller icebergs from fast-moving outlet glaciers (Benn and others, 2017a). Each of these processes produces icebergs of a different size and multiple studies have examined the size distribution of icebergs produced at different sites (e.g. Tournadre and others, 2016; Sulak and others, 2017; Barbat and others, 2019; Neuhaus and others, 2019). Iceberg size distributions can yield a range of useful information: they provide insight on the fundamental processes driving iceberg calving (Benn and Åström, 2018; Cook and others, 2018); they can be useful for assessing iceberg risks to shipping and hydrocarbon infrastructure (Fuglem and others, 1996; Habib and others, 2016); iceberg size distributions have been shown to affect fresh water transport in global ocean models (Stern and others, 2016); and shifts in iceberg size distributions have the potential to alter the timing, magnitude and spatial distribution of ocean freshwater fluxes near glaciers and ice sheets (Enderlin and others, 2016; Rezvanbehbahani and others, 2020).

The most common fragment-size distribution (FSD) observed for icebergs is a power law. Power law distributions have been observed for pan-Antarctic tabular icebergs (Tournadre and others, 2016; Barbat and others, 2019), datasets from Greenland outlet glaciers (Sulak and others, 2017; Crawford and others, 2018; Scheick and others, 2019) and icebergs calved from Columbia Glacier in Alaska (Neuhaus and others, 2019). The observed power law distributions have been attributed to brittle fragmentation (Kirkham and others, 2017; Benn and Åström, 2018) and variations in the power law exponent have been attributed to differences in calving rate, glacier geometry and environmental conditions (Sulak and others, 2017; Neuhaus and others, 2019).

The influence of calving rate and glacier geometry on iceberg size distributions is most apparent for the rapid ice shelf disintegration events observed on the Antarctic Peninsula (Doake and Vaughan, 1991; Glasser and Scambos, 2008). These disintegrations consist of a substantial retreat of the ice shelf calving front as a large portion of the ice shelf area rapidly disintegrates into a multitude of small icebergs. Disintegration events are characterised by a rapid and distinct change in the shape of icebergs produced by the ice shelf, shifting from a typical pattern of irregular releases of large tabular icebergs, to the formation of a large number of narrow ice blocks that topple, collide and grind against each other (Glasser and Scambos, 2008; Scambos and others, 2009). These observed changes in iceberg size distributions suggest that the controls on iceberg calving may evolve with changes in dynamic mass loss from glaciers and ice shelves. However, the precise nature of the relationship between the iceberg size distribution and rate of mass loss is not well understood.

In many cases, iceberg distributions evolve over time as the result of external influences like submarine melting of icebergs in relatively warm ocean waters (Enderlin and others, 2016) and subsequent fragmentation events or crushing (Bouhier and others, 2018; Crawford and others, 2018). These secondary processes can cause the power law exponent to change (Neuhaus and others, 2019) or the distribution to take on a new function (Kirkham and others, 2017). In the ice mélange observed near most Greenland glacier termini throughout at least a portion of the year (Amundson and others, 2010; Enderlin and others, 2014), interactions between icebergs can generate sufficient stress to inhibit new iceberg calving (Amundson and Burton, 2018; Burton and others, 2018). The crushing and grinding of icebergs in this high-stress environment should also alter the iceberg size distributions, such that analysis of iceberg size distributions may yield insights into the influence of ice mélange on calving.

Here we apply and expand on fragmentation theory to provide theoretical interpretations of the size distributions of a range of calving behaviours. The theory produces four possible types of FSD relating to the environment in which fragmentation occurs and the energy involved in the process. We use the term ‘energy’ to refer to the potential energy associated with a calving event. This could be derived from buoyancy for ice fragments below the sea surface, gravitational potential energy above the sea surface, and elastic stress in the ice. Little energy is released in the formation of a large tabular iceberg which has the same thickness as the unbroken shelf, while a considerable amount of energy is dissipated in the release of subaerial or subaqueous icebergs, or during iceberg toppling.

The theoretical FSDs are validated against outputs from Rink Isbræ and Totten Ice Shelf (TIS) numerical model simulations performed using the Helsinki Discrete Element Model (HiDEM) (Åström and others, 2013), as a means to assess climate controls on calving from Greenland and Antarctica, respectively. The theoretical FSDs are also examined with respect to satellite-observed FSDs from a range of Greenland outlet glaciers and Antarctic ice shelves to assess the applicability of the theoretical FSDs to iceberg distributions produced by glaciers and ice shelves with varying geometries and in diverse environmental regimes.

From our results a range of different calving styles emerge: (1) low-energy tabular calving events with an exponential distribution, (2) typical ‘stable’ iceberg calving with a universal power law for smaller icebergs and nonuniversal exponential functions in the larger iceberg size limits, and (3) a high energy ‘unstable’ power-law distribution relating to grinding and crushing. Using both observations and models we demonstrate that iceberg sizes can transition from one distribution to another, either as icebergs are advected down a fjord, or as ice shelf stability decreases under increasing ocean-driven melt and/or hydrofracture. Finally, we discuss the possible implications of changes in iceberg size distributions in a warming climate in which the low-energy calving of large tabular icebergs is gradually replaced by more frequent high-energy calving of smaller-sized icebergs.

2. Fragmentation theory

Our hypothesis is that icebergs that calve from glaciers and ice shelves can be described by the size-distribution of generic fragmentation theory for elastic-brittle materials and a grinding/crushing process. Such size-distribution functions can be defined as the density in the number of fragments of volume v formed during fragmentation events. An FSD function, n(v), that is valid for a large class of elastic-brittle fragmentations can be written for fragments with size between v and v + dv as (Åström, 2006):

(1)$$n( v ) = c_1v^{-\alpha }f_1( {v/c_2} ) + c_3f_2( {v/c_4} ) , \;$$

where α = (2D − 1)/D, and D is the dimension. D = 2 should be used for 2-D fragmentation when the linear dimension of fragments is larger than the thickness of, for example, an ice shelf. Then volume v can be replaced by area a multiplied by a, more or less, constant shelf thickness. D = 3 should be used when calved pieces are smaller than the thickness of the glacier from which they have emerged. Below we differentiate between α _a for D = 2 and α for D = 3. The functions f ₁(x) and f ₂(x) are almost always of exponential form (i.e., exp(− x)), but can, at least in principle, have other forms. That is, they are almost constant for small magnitudes, x < 1, and they decay rapidly for large magnitudes, x > 1. Equation 1 is the elastic-brittle characteristic fragment size distribution, hereafter referred to as the E-BC. A detailed derivation of Eqn 1 is given in Kekäläinen and others (2007).

The function, n(v), on the left side of Eqn (1) is the number density of fragments. This function, integrated over an iceberg size interval, describes the number of observed icebergs of that size in a dataset. This is slightly different from a probability density function, which when integrated over all sizes equals unity. The second term on the right-hand side of Eqn 1 arises from uncorrelated fracture, while the first term arises from correlated fracture that is scale-invariant in nature. For calving, the processes behind the second term can perhaps best be understood by imagining an ice shelf that is sporadically fractured by slowly propagating rifts that eventually delineate large tabular icebergs. Cracks that develop far from each other can be assumed to appear at random locations and independent of each other. This defines a Poisson process and leads to an exponential form for f₂(x), where x may be volume, area, or linear size of fragments (Grady and Kipp, 1985). If such a fragmentation process is very ‘gentle’ with a low rate of energy-release via slowly creeping rifts that do not form crack branches, the FSD will consist of only the second term in Eqn 1. In contrast, if the energy-release rate is high, cracks preferentially form a multitude of branches that readily merge to form additional fragments much smaller in size than those produced by the Poisson process. The FSD of these smaller fragments is described by the first term on the right-hand side in Eqn 1, and we refer to this as ‘correlated’ fracture (Åström, 2006). If the constants c₂ and c₄ in Eqn 1 are of the same order of magnitude, the first term will dominate the second term in Eqn 1. In that case, the FSD reduces to a power law (i.e., first term on the RHS of Eqn 1) with an exponential ‘cut-off’, that may be apparent in observational data depending on the amount and accuracy of observations.

There are three parameters in Eqn 1 that have physical relevance: c₁/c₃, c₂ and c₄. The ratio of c₁/c₃ sets the relative abundance of fragments represented by the two terms in Eqn 1, c₂ depends on the energy that drives the fragmentation and its dissipation, such that low-energy fragmentation and high damping induce small values for c₂. c₄ depends on initiation density of Poisson-type cracks, and in the case that cracks are very low density, the dimensions of the fragmenting system. A detailed description of the meaning of the c₁–c₄ parameters is given in Åström and others (2004).

The E-BC distribution described in Eqn 1 is the result of an initial fragmentation process. This represents the FSD that we expect to find close to the termini of glaciers and ice shelves unless the calving event itself is extremely energetic, such as an ice shelf disintegration. Fragmentation can, however, continue beyond the initial fragmentation as grinding and crushing processes once the icebergs enter the water. Such processes break existing fragments into smaller pieces and thus induce further evolution of the fragment distribution function. This type of process can occur in the densely packed ice mélange found in Greenlandic fjords, where icebergs interact and collide, causing large icebergs to be broken into smaller fragments. A similar size distribution can also be generated during ice shelf disintegration when the capsizing of icebergs releases gravitational potential energy and drives further fragmentation, and the high density of ice fragments causes icebergs to interact (MacAyeal and others, 2003).

The grinding and crushing processes typically retain an approximate scale-invariance, thereby preserving a power law FSD (Weiss, 2003). Thus, the FSD gradually transforms to:

(2)$$n( v ) = c_5v^{-\beta }, \;$$

where β > α, because such a process preferentially breaks up large fragments and thus produces more small fragments. In a grinding process, it is common that all large pieces are destroyed, which means that only a pure power law typically remains. This is the grinding and crushing FSD, hereafter referred to as the GC, which applies to both secondary fragmentation processes and extreme events such as ice shelf disintegration. There does not seem to be universal values for β in the same way as for α. As long as grinding and crushing continues the FSD changes and this manifests as a gradually increasing value of beta. This is best demonstrated when considering a dense mélange extruded from a fjord. Such a process can be viewed as a granular compaction and shear-deformation process and it can, at least partly, be linked to abstract structures called self-similar space-filling packings and their shear stiffness. There are a large number of such structures with different fractal dimensions, which correspond to the exponent β in Eqn 2 and different shear deformation properties (Stäger and Herrmann, 2018) such that there is no universal value for β.

A schematic representation of Eqn 1 and 2 is plotted in Figure 1. As is evident in Figure 1, the two terms on the right-hand side of Eqn 1 separate completely with the first term representing small icebergs (i.e., bergy bits) and the second term the large icebergs. The line representing Eqn 2 demonstrates the steeper slope (i.e., larger exponent) expected for size distributions subject to secondary fragmentation by grinding/crushing. In this particular case β = 3.0.

Fig. 1. Idealized iceberg size distributions. Schematic representations of fragment size distributions from E-BC (Eqn 1, dashed black line) and GC (Eqn 2, dashed green line). For E-BC, distributions from correlated (first term on right-hand side of Eqn 1, solid red line) and uncorrelated (second term on right-hand side of Eqn 1, solid blue line) fracture are plotted separately. The interpretations of the constants c₁ to c₅ are easily extracted from the Figure: c₁ = 100, is the value on the y-axis of the first term of Eqn (1) at v = 1. c₃ = 0.01, is the corresponding value for the second term. c₂ = 300, and c₄ =  10 000, are the cut-off values for the two terms of Eqn 1. c₅ = 500, is the value on the y-axis for Eqn 2 at v = 1. The y-axis represents fragment abundance, while the x-axis is a measure of volume (e.g. (10 m)³). This plot has volume v on the x-axis and therefore α = 5/3. Eqn 2 is plotted with β = 3.0.

In Antarctica, icebergs often calve from ice shelves as tabular icebergs of dimensions larger than the thickness of the shelf. In such a case, the fragmentation process is 2-D and volume (v) can be replaced by area (a) in Eqns 1 and 2 and α = (2D − 1)/D = 1.5 (formally, v could be replaced by ah, where h is a measure of shelf thickness). For smaller icebergs, we must use D = 3, yet there is no efficient direct way to estimate volumes for all icebergs in the ice mélange observed in satellite images of glacier fjords. If the ice mélange is so dense that smaller bergy bits are rafted on top of each other, then in a limiting case the FSD would be best represented by setting D = 2 as only a 2-D cross-section of the 3-D fragmentation process may be visible. If the ice mélange is more loosely packed such that all ice fragments float on the water surface, all fragments formed may be visible. Then the FSD can be obtained by a parameter transformation from volume v to area a, where a ∝ v ^2/3 and the power exponent in Eqn 1 becomes α _a = 2. However, the conversion exponent might be larger than 2/3, because an irregular piece of ice will float on water such that it tends to expose its largest possible surface area to the vertical. That is, an iceberg will float ‘lying down’ on the water. From Rink Glacier in west Greenland, this exponent has been estimated as ~0.77, which gives α _a ≈ 1.87 (Sulak and others, 2017).

For the Greenlandic cases investigated here, we expect that the power-exponent for marine-terminating glaciers is α _a ≈ 2, if measured as area, and α = 5/3 if measured as volume. If there is crushing/grinding of ice, the effective exponents will be larger. For ‘standard’ calving at Antarctic ice shelves, we expect α _a = 1.5. In contrast, for ice shelf disintegration events we accordingly expect β > 1.5. The cut-off functions, f _1,2(x), are far less universal than the power-exponents and we expect them to vary between different glaciers or shelves and also with time. They should, however, typically be of the form f _1,2(x) ~ exp(− x), where x can be the volume (v), area (a), or linear size (l) of icebergs, or some intermediate.

There is at least one more possible functional form for the FSD, the log-normal form, although we do not consider this form here because it is more appropriate for icebergs subject to considerable melt and fragmentation as they drift across the open ocean (Kirkham and others, 2017). This functional form of the FSD can be conceptualized by the following: assume that small icebergs preferentially melt in the open ocean due to their large surface area to volume ratio and that the maximum size of icebergs entering the open ocean is limited by the presence of bathymetric high points within the fjords, then the resulting iceberg volumes (v₀) will be more uniform than near the glacier margins. Assuming that these icebergs promote n random and independent breakup events and that each breakup, for simplicity, cuts an iceberg in half, then the remaining volumes of the disintegrating icebergs will become v₀ /2ⁿ. This means that the logarithm of the iceberg volumes will be proportional to n, and n is the number of random and independent breakup events for each iceberg, resulting in an FSD of log-normal form.

3. Methods

Fragmentation theory was applied to iceberg size distributions extracted from HiDEM (Åström and others, 2013) simulations performed for Rink Isbræ and TIS in Greenland and Antarctica, respectively and numerous satellite datasets. The HiDEM simulations are described below, followed by the observational datasets.

3.1 Discrete element model

We use HiDEM to numerically compute crevasse formation and calving on TIS and Rink Isbræ (Rink). Discrete element models (DEMs) are useful for simulating fracture on ice shelves and glaciers because they discretise the ice shelf geometry into separate blocks (or elements) joined by bonds and allow these bonds to break once a critical stress is exceeded. This framework allows fractures and calving events to form explicitly within the model domain (Åström and others, 2013, 2014; Riikilä and others, 2014). There are a few important limiting factors of DEM-simulations for calving: a lack of viscous deformation and a limitation in the longest times that can be simulated of ~10–100 min, depending on the number of DEM-blocks and available computational capacity. This means that DEM simulations are good at capturing elastic-brittle calving events, but not the accumulation of damage over viscous timescales (days and longer).

Our numerical simulations utilize observed glacier geometries that are discretised into equally sized spherical blocks connected by inelastic beams. Blocks are arranged in a densely-packed HCP-lattice. The shape of this lattice introduces a weak directional bias in the elastic and fracture properties of the ice, however this is easily broken by the propagating cracks, as evidenced by the results. At the beginning of a fracture simulation, the ice is assumed to contain a small density of randomly scattered pre-existing cracks. The movement of the ice blocks is computed using a discrete version of Newton's equation of motion and iteration of time steps. A buoyancy force is applied to each block of ice below the water surface and a gravitational force is applied to blocks above the water surface.

Since the model cannot cover the whole glacier or ice shelf due to computational cost, some of the boundaries of the domain do not float freely, but experience pressure from the land-based ice. This pressure is reproduced by scaling the hydrostatic pressure at the back of the simulated area. At each time step, the stress acting on each beam is calculated and if that stress exceeds a critical threshold of the order of 1 MPa the beam is broken. The model is then iterated forward with a time step of 0.0001 s.

We implement into HiDEM a geometry for Rink that is a specific and preliminary version of BedMachine Greenland (Morlighem and others, 2017). Rink is a 500–800 m thick tidewater glacier in western Greenland that has maintained a relatively stable terminus position in the last several decades. We compute fracture and calving at Rink terminus using ~5 km × 5 km domain. We use DEM blocks with a 30 m diameter which requires just under 700 000 particles with ~41 million connections to construct the glacier in Figure 2b. As observed, glacier geometries are almost by definition stable against immediate calving, we alter the geometry to induce calving and produce measurable FSDs. For Rink we perform two simulations with distinct modifications: (1) sea ice is introduced in front of the terminus which is fragmented by the forward moving glacier and at the same time the ice-cliff is partly crushed, and (2) melt-undercutting is applied at the base of the terminus.

Fig. 2. Rink Isbræ Fragment Size Distributions. (a) Computed and observed FSDs at Rink fitted to Eqn 1 (E-BC) and 2 (GC) with D = 3 (α = 5/3). For the observed data c₃ = 0, and c₂ ~ 10⁸ m³, which means the largest icebergs are ~1 km². For the sea ice case c₄ ~ 10⁷ m³ and c₂ ~ 10⁶ m³. The exact value of β depends on the precise method used for fitting, but in all cases, it is close to 3, which is used in the plot. (b) A snapshot of the Rink HiDEM simulation that produced the FSDs with sea ice in Figure 2a. Basal crevasses were formed during the simulation but cannot be observed at the ice surface. Dark-blue sea water is visible through cracks in the sea ice in front of the glacier terminus. The contours on the glacier surface represent 30 m (i.e., DEM block size) elevation increments. The dark-grey margins indicate bedrock and a shadow is cast from the fjord-side hills.

TIS is ~5000 km² floating part of the terminus region of Totten Glacier in East Antarctica. Totten drains a considerable part of the East Antarctic ice sheet, with a drainage basin that extends ~1000 km inland. We simulate ~50 km × 50 km of the terminus of TIS including parts of its grounded sides with ~5 million particles with 30 million connections and a DEM-size of 50 m. The geometry used for HiDEM model setup is derived from aerogeophysical data collected by the International Collaboration for Exploration of the Cryosphere through Aerogeophysical Profiling (ICECAP) project which has flown a dense network of flights over the TIS. The surface elevation (Blankenship and others, 2015) and ice thickness (Blankenship and others, 2012) data are interpolated using TELVIS (Thickness Estimation by a Lagrangian Validated Interpolation Scheme) (Roberts and others, 2011), to produce a geometry which is glaciologically self-consistent (Roberts and others, 2018). The surface area of the model is defined by a calving front manually digitised from a Landsat 7 image from November 2010.

As for Rink, the initial geometry is inherently stable against immediate calving, therefore we alter the geometry to induce calving. For TIS we perform three types of simulations: (1) one long simulation with enhanced driving pressure to force TIS advection and induce fragmentation, (2) shorter simulations with the ice shelf geometry altered by enhanced melting and (3) simulations that are the same as (2) except with unbalanced buoyancy. These simulations represent extreme versions of natural ice shelf processes which might be produced for example by increased rates of ocean melt and by changes in surface meltwater loading.

In the first experiment we apply forcing on TIS by using the maximum pressure that can be applied without crushing the glacier at the upstream face of the simulation domain. This induces forward motion and calving of TIS. In reality, this process is much too slow to be captured by HiDEM. We, therefore, apply as little dissipation as possible and simulate for as long as is feasible (~100 min) to obtain reasonable FSDs (Cook and others, 2018).

In the second set of numerical experiments, we test the effect of thinning of the ice shelf using basal melt rates, $\dot{m}( {x, \;y} )$, from regional ocean modelling (Gwyther and others, 2014) averaged over a period 1992–2007 (Fig. 5d). Since the ice shelf did not thin over this period (Paolo and others, 2015), we apply this melt-rate in increasing multiples to thin the model geometry, keeping the glacier at neutral buoyancy everywhere. Melt is stopped if the local ice thickness reaches a minimum of the initial subaerial height above sea level. Ideally, an updated model geometry would be derived from a coupled ice-ocean model, allowing the ice shelf geometry to synchronously adapt to basal melt; however, this simplistic approach has been chosen as a first-order approximation to investigate the fragmentation of a significantly thinning shelf geometry. These additional simulations are performed over the same domain as the initial simulation (1) and for a shorter time.

A final set of experiments investigate the role of unbalanced buoyancy in iceberg calving. This process has been hypothesised to lead to increased iceberg calving, as draining surface ponds can leave areas of an ice shelf out of hydrostatic balance (Banwell and others, 2013; Banwell and Macayeal, 2015). In order to create an unbalanced modelled ice shelf, we repeat the thinning experiments but remove ice from the base without any hydrostatic adjustment.

3.2 Observed iceberg size distributions

Iceberg size data are extracted from satellite images for several locations around Greenland and Antarctica.

3.2.1 Amundsen sea

Iceberg surface areas are automatically digitized using a segmentation algorithm (Mazur, 2017) from 74 ENVISAT ASAR Wide Swath Mode (WSM) Level 1b images, one per month from February 2006–December 2012. Although the ENVISAT imagery is relatively low resolution (75 m) the high detectability of icebergs in SAR data and the efficient automated detection method mean that icebergs of only two pixels can be detected (Mazur, 2017; Mazur and others, 2017), producing a large dataset of 7397 icebergs with areas between 0.09 and 656 km².

3.2.2 Larsen B

Iceberg surface areas are manually digitised from a Landsat 7 ETM + (Enhanced Thematic Mapper Plus) image from April 2002, which provides a spatial resolution of 30 m. This is shortly after the disintegration of the ice shelf in February–March 2002, and many of the calved ice blocks remain in the region (Glasser and Scambos, 2008). The Larsen B dataset consists of 1242 icebergs with areas between 0.002 and 22.2 km².

3.2.3 Totten ice shelf

Iceberg surface areas are manually digitised from multiple satellite images, since the release of icebergs is much less frequent than during ice shelf disintegration events. For this dataset icebergs were digitised from three Landsat 7 images from 23 November 2009, 26 November 2010 and 29 November 2011. The dataset consists of 370 icebergs with areas between 0.001 and 5.8 km².

3.2.4 Wilkins ice shelf

Iceberg surface areas are manually digitised from a Formosat-2 image acquired on 8 March 2008, immediately after the break-up event which lasted from 28 February to 6 March. The high spatial resolution of the Formosat image (2 m) allows for much smaller icebergs to be digitised. The Wilkins dataset consists of 2451 icebergs with areas between 2.5 × 10⁻⁴ and 1.2 km².

3.2.5 Greenland fjords

Landsat 8 image-derived iceberg size distributions from Rink, Kangerlussuup Sermia and Sermilik fjords are from Sulak and others (2017). These size distributions span tens of kilometres from the glacier termini and, therefore, include icebergs that have been considerably altered by ablation processes. To extract iceberg size distributions from the dense ice mélange, we divide high-resolution digital elevation models of the ice mélange into along-fjord bins and extract ice elevation distributions from each bin. Digital elevation models are created from DigitalGlobe, Inc. imagery using the NASA Ames Stereo Pipeline (Shean and others, 2016). The elevation distributions are converted to iceberg size distributions following the approach of Enderlin and others (2016). The along-fjord bin size is restricted to the maximum width of icebergs observed in each fjord (~1 km) and the size classes are limited by the vertical uncertainty (~3 m) of the elevations.

4. Results

The FSDs produced by HiDEM simulations of Rink Isbræ and TIS are compared to iceberg size data extracted from satellite images for several locations around Greenland and Antarctica, respectively. Specifically, size distributions were extracted for icebergs near Rink Isbræ and Kangerlussuup Sermia and in Sermilik and Ilulissat fjords in Greenland as well as Totten, Larsen B and Wilkins ice shelves and the Amundsen Sea embayment in Antarctica.

4.1 Greenland glaciers

In our Rink geometry for HiDEM, a large part of the terminus is slightly below flotation and thus experiences buoyant uplift during computation. This induces a set of distinct crevasses forming at the base and propagating upwards (Benn and others, 2017b). In the simulations, these cracks are almost enough to form a tabular iceberg in the terminus region but propagation stops before a freely floating iceberg can be fully formed. Tabular icebergs of a similar size to those that nearly form in our simulations have been observed to form at the same location, albeit rarely, at Rink (Medrzycka and others, 2016).

A straight-forward implementation of the observed Rink geometry in HiDEM produces, as expected, a FSD with only a few minor icebergs (<10 DEM blocks, or <0.125 km³, in volume). Two further cases which produce more typical FSDs are displayed in Figure 2a. In the first case we introduce sea ice in front of the terminus (Fig. 2b). In the second case there is no sea ice but instead we introduce a melt-undercut (~100 m) near the base at the front of the glacier. In the sea ice scenario, a gentle spontaneous forward sliding motion of the glacier during the simulation induces break-up of the sea ice and small ice fragments to calve off the edge of the glacier. The simulated FSD is fit with Eqn 1 (E-BC) with D = 3 (Fig. 2a red markers). In this simulation, the rather noticeable contribution of the exponential term, f ₂(v), originates from uncorrelated fracture in the broken sea ice, while the power law component of Eqn (1) is dictated primarily by icebergs. Figure 2a also displays the FSD and best-fit for the under-cut geometry HiDEM simulation. In contrast with the sea ice configuration, the FSD can be modelled with a steep power law alone, representative of GC (i.e., Eqn 2). In this case, the large undercut makes the computed ice-cliff unstable and a high-energy calving event produce a GC distribution as the calving ice is crushed during a rapid disintegration event.

Together with an observed FSD from Rink, the distributions displayed in Figure 2a have the same general structures as the schematic FSDs in Figure 1. The modelling exercise for Rink Isbræ demonstrates three fragmentation processes presented above: the observed iceberg calving at Rink can be described by correlated fracture and the first term in Eqn 1, the rather low-energy fragmentation of the sea ice produces uncorrelated fracture described by the second term in Eqn 1, and the high-energy calving of the collapsing front in the under-cut case is described by Eqn 2 with β ≈ 3.

Observational data for icebergs calved from several other Greenlandic glaciers show the same type of FSDs. Figure 3a includes optical image-derived iceberg volumes from summer (July–September) 2013–2015 for Sermilik and Kangerlussuup Sermia fjords and FSDs from Eqn 1 with D = 3 that were fit to the observations. Sermilik observations are separated into open water distributions and distributions within the ice mélange adjacent to the terminus, hereafter referred to as Sermilik Fjord and Sermilik Fjord Mélange, respectively. In previous work (Sulak and others, 2017), the data in Figure 3a were described by pure power laws, resulting in slightly different power exponents for the datasets. By fitting Eqn 1 instead of pure power laws, we find that the power-exponents remain the same for the datasets and that the differences instead appear in the f₁ and f₂ cut-off functions. Effectively, however, we observe the same phenomenon as in the initial data analysis: Kangerlussuup Sermia fjord has the lowest density of small icebergs (~10⁴–10⁵ m³) and is the only dataset that needs to be fitted with both the correlated and uncorrelated terms in Eqn 1 because of its relatively high fraction of large icebergs (~10⁷ m³). In contrast, the Sermilik Fjord FSD has the highest density of small icebergs and an exponential cut-off that steepens the slope of the size distribution more dramatically than for the nearby Sermilik Fjord Mélange, indicating a relative lack of large icebergs in the open portion of the fjord.

Fig. 3. FSDs for Greenland tidewater glaciers. (a) Eqn 1 with D = 3 (α = 5/3) fitted to FSD observation data for Kangerlussuup Sermia and open water and ice mélange in Sermilik Fjord. c₂, c₄ ~ 10⁷ m³ (b) Sermilik Fjord Mélange FSDs in October 2014 at 1 km and 4 km from terminus fitted to the E-BC (Eqn 1) distribution at 1 km, and the GC (Eqn 2) distribution at 4 km. (c) Ilulissat Isfjord FSDs in April 2014 fitted to the E-BC (Eqn 1) distribution with the power law term for correlated fracture absent and August 2015 fitted to the E-BC distribution with the Poisson term for uncorrelated fracture absent. (d) Ilulissat Isfjord FSDs on 8 August 2015 at 1500–12 500 m from terminus compared to power laws with α _a = 2.0 (E-BC), and β = 2.4 (GC).

Figure 3b shows FSDs inferred from elevation distributions in high-resolution digital elevation models for the Sermilik Fjord Mélange in October 2014 (Enderlin and others, 2016). Since the ice mélange is consistently kilometres-long, the iceberg size distributions were parsed into ~1 km bins from the glacier terminus to the end of the observational domain. Here we show FSDs computed ~1–2 km and ~4–5 km from the terminus, each with an area of ~3 km². There is a noticeable along-fjord change in FSDs: near the terminus, the FSD can be described using the E-BC distribution but farther from the calving front the FSD is a very steep power-law with β ≈ 3.2, consistent with the GC distribution.

Figure 3c shows FSDs for Ilulissat Isfjord on 19 April 2014 and on 13 August 2015, both obtained near the terminus of Sermeq Kujalleq (Jakobshavn Isbræ) using the same methods as Enderlin and others (2016). In April, at the end of the winter, the glacier typically has a floating tongue that calves in a rather gentle fashion with little potential energy released in each calving event and we observe the E-BC distribution without the power law term for correlated fracture. In contrast, the iceberg FSD in August is produced by a grounded glacier that calves smaller icebergs in a much more energetic fashion and we observe an E-BC distribution in which the power law term for correlated fracture dominates.

Figure 3d shows FSDs obtained on 8 August 2015 from Ilulissat Isfjord parsed in 1 km bins from ~1 to 13 km from the terminus, with areas ranging from ~5 to 10 km². For these data, the power law exponent gradually increases in magnitude from 2.0 at 500 m to 2.4 at 5500 m. This suggests gradual grinding and crushing of icebergs within the ice mélange, as larger icebergs are continuously broken into smaller pieces. Although the ice mélange persists for another ~7 km, the exponents remain the same throughout the remainder of the ice mélange to the outside of the fjord (Scheick and others, 2019) as the interaction between icebergs diminishes then largely ceases.

4.2 Antarctic ice shelves

The TIS geometry we used in HiDEM is close to neutral buoyancy and as expected, not much calving takes place without altering the geometry or introducing an external forcing. To induce calving, we implemented as large a pressure at the back of TIS we could without crushing it (~4 times the depth-dependent hydrostatic pressure) and computed the fracture of TIS as the ice shelf advected. Eventually a fracture pattern emerged that resembles that of TIS and an FSD that is almost identical to the observed FSD for large icebergs (Cook and others, 2018).

Figure 4a displays the HiDEM – and remote sensing –derived TIS iceberg size distributions fitted with E-BC distributions. For the largest icebergs (>0.3 km²), the observed and computed FSDs coincide, and can be fitted with a 1-D Poisson term: exp(l), where l = √a. The FSDs for the largest icebergs also can be reasonably well described by a power law with α _a = 1.5 (Cook and others, 2018). The HiDEM-computed and the observed FSDs deviate significantly in the small-area range (10⁴–10⁵ m²) in Figure 4a, indicating that the computed calving is much more energetic than observed. This is not very surprising since the short simulation times attainable with HiDEM forces us to simulate with much higher advection velocities than TIS has in real life.

Fig. 4. Iceberg area distributions from Antarctica. (a) Observed and computed FSDs at TIS. The observed FSD is fitted with a 1-D Poisson term: exp(l), and the computed FSD is fitted with the E-BC function and α _a = 2.0. (b) Abundance of icebergs as functions of iceberg area from the Amundsen Sea. Data are shown separately for different months of the year and fitted with the power law of Eqn 1 with D = 2 (α _a = 1.5, red line).

Previous work measuring the size of icebergs observed in Antarctica indicates that they follow a power law distribution, with an exponent of 1.5 as expected for 2-D E-BC fragmentation (Tournadre and others, 2016). This is confirmed by examining iceberg FSDs in the Amundsen Sea region from 2006 to 2012, using automated detection methods from synthetic aperture radar imagery (Mazur, 2017; Mazur and others, 2017), shown in Figure 4b. There are no significant seasonal variations in FSDs (Fig. 4b) and a power law distribution with exponent 1.5 and cut-off size of ~5 km², with no significant contribution from the Poisson term for uncorrelated fracture, reasonably matches the data for large iceberg sizes (>0.3 km²). Together Figures 4a, b indicate that Antarctic icebergs may be fragile and easily break apart in the scale-invariant fashion characteristic for elastic-brittle materials, but that ice-shelf rift propagation also may be so slow and gentle that a Poisson process dominates.

Given the over-abundance of small icebergs in the initial TIS HiDEM simulation, the second numerical experiment for TIS included a set of shorter computations without any extra pressure pushing the ice shelf. Instead, the ice shelf was thinned via basal melting. Figures 5a–c show calving and disintegration of TIS at 0 (i.e., Control), 50 and 100 times the current annual basal melt rate ($\dot{m}( {x, \;y} )$, average of ~5 ma⁻¹). At 100 times basal melt the fastest melting regions of the domain have reached the minimum allowed thickness. In the slow-melting regions the ice is still 100–200 m thick. As thinning by basal melt increases, the rate of fracture and iceberg production increases. For the 100$\dot{m}( {x, \;y} )$ simulation, the modelled ice shelf area is reduced by 8% due to enhanced calving.

Fig. 5. Modelled fragmentation of Totten Ice Shelf under basal melt. Red represents large icebergs, followed by black (~0.1 km²) and violet (~10⁴ m²). Single element ice blocks (~10³ m²) are coloured in light blue. Grey represents the largest unbroken ice block. (a–c) Model runs under increasing melt, with neutral buoyancy. (d) modelled average annual melt-rate for TIS 1992–2007. (e and f) Model runs under increasing melt with unbalanced buoyancy.

Figure 6a shows the FSD for the different melt rate simulations, with the power law distributions of exponent 1.5 (E-BC) and 2.5 (GC) included for reference. The 1.5 exponent is consistent with calving observation in Antarctica (Stern and others, 2016; Tournadre and others, 2016) and the 2.5 exponent seems to be the approximate value for shelf collapse. As the melt-induced thinning increases, the cross-over point between E-BC and GC distributions (square markers in Fig. 6) shifts towards larger iceberg sizes and the proportion of small icebergs increases, suggesting that grinding and crushing become important for progressively larger iceberg sizes.

Fig. 6. Fragment size distributions (FSD) under increasing basal melt. (a) Modelled iceberg sizes under increasing thinning by basal melt with model in neutral buoyancy, compared to observed size distribution for the disintegration of the Wilkins Ice Shelf. (b) Modelled iceberg sizes under increasing thinning by basal melt with model in unbalanced buoyancy, compared to observed size distribution for the disintegration of the Larsen B Ice Shelf. Dashed lines show the theoretical slopes for stable calving (area^−1.5) and disintegration (area^−2.5) FSDs. Square boxes indicate the cross-over point between stable-type (E-BC) and disintegration-type (GC) FSD.

Given observations of surface meltwater build-up on Antarctic ice shelves prior to their disintegration, we also modified the enhanced basal melt simulations to account for surface loading effects (i.e., unbalanced buoyancy). In these simulations, as thinning increases, the ice shelf becomes unstable and a large portion disintegrates into icebergs (coloured shading in Figs 5e, f). The incidence of small icebergs (shaded light blue in Fig. 6) is much higher than in experiments with neutral buoyancy. In experiments with both 50$\dot{m}( {x, \;y} )$ and 100$\dot{m}( x, \;y$) applied to the model geometry, the front of TIS entirely decouples from the rest of the shelf upstream, roughly simulating a rapid ice shelf disintegration event (Figs 5e, f). As with the buoyancy-balanced thinning experiments, the cross-over point between E-BC and GC distributions (square markers in Fig. 6b) shifts towards larger iceberg sizes and the proportion of small icebergs increases, suggesting that grinding and crushing become important for progressively larger iceberg sizes. This modelled change in the iceberg size distribution is supported by satellite observations of iceberg sizes during rapid ice-shelf disintegration. Over the past three decades, there have been a number of ice shelf disintegration events along the Antarctic Peninsula, including the Wordie (Doake and Vaughan, 1991), Wilkins (Scambos and others, 2009) and Larsen Ice Shelves (Rott and others, 1996; Glasser and Scambos, 2008). Figures 6a and b show FSD measured during the disintegrations of both Wilkins and Larsen B ice shelves, respectively. The cross-over between E-BC and GC is near the upper limit of the Wilkins and Larsen B FSDs, supporting previous observations of abundant small ice fragments formed during ice shelf disintegration events (Glasser and Scambos, 2008; Scambos and others, 2009), consistent with capsizing ice shelf fragments after ice shelf disintegration (MacAyeal and others, 2003; Macayeal and others, 2011; Burton and others, 2013).

5. Implications for iceberg calving in a changing climate

Using satellite observations from a range of Greenland glaciers, we find three distinct iceberg size distributions that are compatible with fragmentation theory. Icebergs observed near calving fronts can all be described by a FSD characteristic of elastic-brittle processes. Nearly all elastic-brittle distributions included a power law component indicative of correlated fracture branching, implying that fracturing of the glacier during calving of relatively small icebergs is pervasive. When the terminus geometry is favourable for the formation of large tabular icebergs, as observed near the terminus of Sermeq Kujalleq in April 2019 and Helheim Glacier in October 2014, the exponential part of the E-BC distribution becomes important, implying that fractures governing calving grew independently of each other (uncorrelated fracture).

While the elastic-brittle distribution dominates near the calving front, our observations also suggest that iceberg FSDs often change in character down-fjord as the result of grinding and crushing within the ice mélange that occupies many fjords. The HiDEM simulations of Rink Glacier indicate that a grinding-crushing distribution can also be produced by a grounded marine-terminating glacier subject to strong terminus undercutting. In our simulations, deep undercutting led to high-energy calving of the collapsing front, promoting fragmentation and the generation of smaller icebergs. Thus, our model results suggest that the increase in undercutting estimated for several glaciers in west Greenland over the last two decades (Rignot and others, 2016) has likely altered the iceberg size distributions near Greenland termini.

For Antarctic ice shelves with stable calving (TIS and Amundsen Sea) we find agreement between observed, computed and theoretical elastic-brittle FSDs for icebergs larger than ~0.3 km². For smaller icebergs, there is a significant discrepancy between the observed and the computed FSDs for TIS. There are three identifiable reasons for this: (1) in the computations, the glacier must be forced to move much faster than in nature to speed up calving simulations, which causes more energetic fragmentation; (2) in nature much energy is dissipated on smaller scales than the DEM size (50 m). If simulation could be run with smaller particles, a similar effect would appear; and (3), the manually extracted observational datasets for TIS are missing many of the smallest icebergs (Cook and others, 2018). Consequently, all forms of FSDs (computed, observed and theoretical) for Antarctic ice shelves are consistent with an elastic-brittle distribution for large icebergs and there is probably also a grinding-crushing component to these FSDs on scales smaller than those contained within current observational data.

In line with the Greenlandic results, the HiDEM simulations of TIS also indicate a transition in the FSDs. In this case, the transition from elastic-brittle to grinding-crushing behaviour does not appear as a result of iceberg interactions in a dense mélange, but instead as a result of a transition from stable calving to ice shelf disintegration events induced by thinning. The simulated changes in iceberg sizes match closely the FSDs produced during the rapid disintegration of the Larsen B and Wilkins ice shelves. It is worth noticing that HiDEM shelf disintegration for balanced buoyancy is only partial as was observed for Wilkins (Scambos and others, 2009). If, in contrast, an ice shelf suddenly switches from balanced to unbalanced buoyancy, a catastrophic disintegration may appear, as in the case of Larsen B (Glasser and Scambos, 2008).

Our results provide a theoretical underpinning to the variety of iceberg size distributions previously reported. At the point of iceberg calving, size distributions typically follow a power law function with a power exponent of –(2D – 1)/D, where D = 2 for tabular icebergs and D = 3 for smaller ice debris, with an exponential ‘cut-off’ for the largest iceberg sizes. However, this function can be modified as glacier and ice shelf geometry evolves in response to climate change, as for the rapid ice shelf disintegration events observed on the Antarctic Peninsula.

The application of the theoretical elastic-brittle and grinding-crushing distributions to observed iceberg size distributions opens the possibility of using iceberg size distributions to investigate the physical processes controlling calving and diagnosing the vulnerability of glaciers and ice shelves to undergo rapid retreat. Additionally, although limited in spatial and temporal scope, our modelled and observed size distributions suggest that climate change has and will continue to alter the size distributions of icebergs. Alteration of iceberg size distributions has the potential to influence the timing, magnitude and spatial distribution of fresh water released from icebergs into ocean basins near glaciers and ice shelves (Enderlin and others, 2016) as well as the trajectory of iceberg drift. Thus, additional application of the theoretical elastic-brittle and grinding-crushing fragment size distributions to observed iceberg size distributions has the potential to yield insights into feedbacks between glaciers and oceans in a changing climate.