1. Introduction
Iceberg production is a key mass loss process for glaciers and ice sheets around the world (Benn and others, Reference Benn, Warren and Mottram2007). 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, Reference Pollard, Deconto and Alley2015; Choi and others, Reference Choi, Morlighem, Wood and Bondzio2018). One complicating factor is that iceberg calving can occur as the result of several mechanisms (Benn and others, Reference Benn, Warren and Mottram2007), from the release of large tabular icebergs from floating ice shelves (Joughin and MacAyeal, Reference Joughin and MacAyeal2005), to the toppling of smaller icebergs from fast-moving outlet glaciers (Benn and others, Reference Benn, Cowton, Todd and Luckman2017a). 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, Reference Tournadre, Bouhier, Girard-Ardhuin and Remy2016; Sulak and others, Reference Sulak, Sutherland, Enderlin, Stearns and Hamilton2017; Barbat and others, Reference Barbat, Rackow, Hellmer, Wesche and Mata2019; Neuhaus and others, Reference Neuhaus, Tulaczyk and Begeman2019). Iceberg size distributions can yield a range of useful information: they provide insight on the fundamental processes driving iceberg calving (Benn and Åström, Reference Benn and Åström2018; Cook and others, Reference Cook2018); they can be useful for assessing iceberg risks to shipping and hydrocarbon infrastructure (Fuglem and others, Reference Fuglem, Jordaan, Crocker, Cammaert and Berry1996; Habib and others, Reference Habib, Hicks, Stuckey and King2016); iceberg size distributions have been shown to affect fresh water transport in global ocean models (Stern and others, Reference Stern, Adcroft and Sergienko2016); 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, Reference Enderlin, Hamilton, Straneo and Sutherland2016; Rezvanbehbahani and others, Reference Rezvanbehbahani, Stearns, Keramati, Shankar and van der Veen2020).
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, Reference Tournadre, Bouhier, Girard-Ardhuin and Remy2016; Barbat and others, Reference Barbat, Rackow, Hellmer, Wesche and Mata2019), datasets from Greenland outlet glaciers (Sulak and others, Reference Sulak, Sutherland, Enderlin, Stearns and Hamilton2017; Crawford and others, Reference Crawford, Mueller, Desjardins and Myers2018; Scheick and others, Reference Scheick, Enderlin and Hamilton2019) and icebergs calved from Columbia Glacier in Alaska (Neuhaus and others, Reference Neuhaus, Tulaczyk and Begeman2019). The observed power law distributions have been attributed to brittle fragmentation (Kirkham and others, Reference Kirkham2017; Benn and Åström, Reference Benn and Åström2018) and variations in the power law exponent have been attributed to differences in calving rate, glacier geometry and environmental conditions (Sulak and others, Reference Sulak, Sutherland, Enderlin, Stearns and Hamilton2017; Neuhaus and others, Reference Neuhaus, Tulaczyk and Begeman2019).
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, Reference Doake and Vaughan1991; Glasser and Scambos, Reference Glasser and Scambos2008). 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, Reference Glasser and Scambos2008; Scambos and others, Reference Scambos2009). 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, Reference Enderlin, Hamilton, Straneo and Sutherland2016) and subsequent fragmentation events or crushing (Bouhier and others, Reference Bouhier, Tournadre, Rémy and Gourves-Cousin2018; Crawford and others, Reference Crawford, Mueller, Desjardins and Myers2018). These secondary processes can cause the power law exponent to change (Neuhaus and others, Reference Neuhaus, Tulaczyk and Begeman2019) or the distribution to take on a new function (Kirkham and others, Reference Kirkham2017). In the ice mélange observed near most Greenland glacier termini throughout at least a portion of the year (Amundson and others, Reference Amundson2010; Enderlin and others, Reference Enderlin2014), interactions between icebergs can generate sufficient stress to inhibit new iceberg calving (Amundson and Burton, Reference Amundson and Burton2018; Burton and others, Reference Burton, Amundson, Cassotto, Kuo and Dennin2018). 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, Reference Åström2013), 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, Reference Åström2006):
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 1(x) and f 2(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 (Reference Kekäläinen, Åström and Timonen2007).
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 f2(x), where x may be volume, area, or linear size of fragments (Grady and Kipp, Reference Grady and Kipp1985). 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, Reference Åström2006). If the constants c2 and c4 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: c1/c3, c2 and c4. The ratio of c1/c3 sets the relative abundance of fragments represented by the two terms in Eqn 1, c2 depends on the energy that drives the fragmentation and its dissipation, such that low-energy fragmentation and high damping induce small values for c2. c4 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 c1–c4 parameters is given in Åström and others (Reference Åström, Linna, Timonen, Møller and Oddershede2004).
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, Reference MacAyeal, Scambos, Hulbe and Fahnestock2003).
The grinding and crushing processes typically retain an approximate scale-invariance, thereby preserving a power law FSD (Weiss, Reference Weiss2003). Thus, the FSD gradually transforms to:
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, Reference Stäger and Herrmann2018) 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 c1 to c5 are easily extracted from the Figure: c1 = 100, is the value on the y-axis of the first term of Eqn (1) at v = 1. c3 = 0.01, is the corresponding value for the second term. c2 = 300, and c4 = 10 000, are the cut-off values for the two terms of Eqn 1. c5 = 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)3). 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, Reference Sulak, Sutherland, Enderlin, Stearns and Hamilton2017).
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, Reference Kirkham2017). 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 (v0) 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 v0 /2n. 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, Reference Åström2013) 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, Reference Åström2013, Reference Åström2014; Riikilä and others, Reference Riikilä, Tallinen, Åström and Timonen2014). 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, Reference Morlighem2017). 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 c3 = 0, and c2 ~ 108 m3, which means the largest icebergs are ~1 km2. For the sea ice case c4 ~ 107 m3 and c2 ~ 106 m3. 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 km2 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, Reference Blankenship, Kempf, Young and Lindzey2015) and ice thickness (Blankenship and others, Reference Blankenship, Kempf and Young2012) data are interpolated using TELVIS (Thickness Estimation by a Lagrangian Validated Interpolation Scheme) (Roberts and others, Reference Roberts2011), to produce a geometry which is glaciologically self-consistent (Roberts and others, Reference Roberts2018). 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, Reference Cook2018).
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, Reference Gwyther, Galton-Fenzi, Hunter and Roberts2014) averaged over a period 1992–2007 (Fig. 5d). Since the ice shelf did not thin over this period (Paolo and others, Reference Paolo, Fricker and Padman2015), 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, Reference Banwell, MacAyeal and Sergienko2013; Banwell and Macayeal, Reference Banwell and Macayeal2015). 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, Reference Mazur2017) 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, Reference Mazur2017; Mazur and others, Reference Mazur, Wåhlin and Krężel2017), producing a large dataset of 7397 icebergs with areas between 0.09 and 656 km2.
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, Reference Glasser and Scambos2008). The Larsen B dataset consists of 1242 icebergs with areas between 0.002 and 22.2 km2.
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 km2.
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−4 and 1.2 km2.
3.2.5 Greenland fjords
Landsat 8 image-derived iceberg size distributions from Rink, Kangerlussuup Sermia and Sermilik fjords are from Sulak and others (Reference Sulak, Sutherland, Enderlin, Stearns and Hamilton2017). 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, Reference Shean2016). The elevation distributions are converted to iceberg size distributions following the approach of Enderlin and others (Reference Enderlin, Hamilton, Straneo and Sutherland2016). 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, Reference Benn2017b). 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, Reference Medrzycka, Benn, Box, Copland and Balog2016).
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 km3, 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 2(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, Reference Sulak, Sutherland, Enderlin, Stearns and Hamilton2017), 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 f1 and f2 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 (~104–105 m3) 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 (~107 m3). 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. c2, c4 ~ 107 m3 (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, Reference Enderlin, Hamilton, Straneo and Sutherland2016). 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 km2. 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 (Reference Enderlin, Hamilton, Straneo and Sutherland2016). 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 km2. 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, Reference Scheick, Enderlin and Hamilton2019) 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, Reference Cook2018).
Figure 4a displays the HiDEM – and remote sensing –derived TIS iceberg size distributions fitted with E-BC distributions. For the largest icebergs (>0.3 km2), 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, Reference Cook2018). The HiDEM-computed and the observed FSDs deviate significantly in the small-area range (104–105 m2) 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, Reference Tournadre, Bouhier, Girard-Ardhuin and Remy2016). 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, Reference Mazur2017; Mazur and others, Reference Mazur, Wåhlin and Krężel2017), 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 km2, with no significant contribution from the Poisson term for uncorrelated fracture, reasonably matches the data for large iceberg sizes (>0.3 km2). 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−1). 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 km2) and violet (~104 m2). Single element ice blocks (~103 m2) 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, Reference Stern, Adcroft and Sergienko2016; Tournadre and others, Reference Tournadre, Bouhier, Girard-Ardhuin and Remy2016) 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, Reference Doake and Vaughan1991), Wilkins (Scambos and others, Reference Scambos2009) and Larsen Ice Shelves (Rott and others, Reference Rott, Skvarca and Nagler1996; Glasser and Scambos, Reference Glasser and Scambos2008). 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, Reference Glasser and Scambos2008; Scambos and others, Reference Scambos2009), consistent with capsizing ice shelf fragments after ice shelf disintegration (MacAyeal and others, Reference MacAyeal, Scambos, Hulbe and Fahnestock2003; Macayeal and others, Reference Macayeal, Abbot and Sergienko2011; Burton and others, Reference Burton, Mac Cathles and Wilder2013).
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, Reference Rignot2016) 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 km2. 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, Reference Cook2018). 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, Reference Scambos2009). 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, Reference Glasser and Scambos2008).
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, Reference Enderlin, Hamilton, Straneo and Sutherland2016) 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.
Acknowledgements
Computing resources were provided by CSC – IT Center for Science Ltd. and National Computational Infrastructure grant m68. Imagery from Formosat-2 was provided through the NATO Support and Procurement Organisation in consultation with Dr Chen-Chien Liu. This work was supported by the Australian Government's Business Cooperative Research Centres Programme through the Antarctic Climate and Ecosystems Cooperative Research Centre (ACE CRC) and the Australian Research Council's Special Research Initiative for the Antarctic Gateway Partnership (project ID SR140300001). During this publication, Sue Cook received support from the Australian Government as part of the Antarctic Science Collaboration Initiative program. Antarctic iceberg data are available at: https://doi.org/10.25959/5e41d14444c88. Greenland iceberg data available on request from the authors.
Open access
