Calving processes

Calving and related dynamic processes are poorly represented in the current generation of ice-sheet models (Benn and others, Reference Benn2017a). Owing in large part to a lack of observational data (e.g. James and others, Reference James, Murray, Selmes, Scharrer and OLeary2014), the true mechanisms of calving are still largely unknown and thus the development of a universal calving law remains unresolved. Ice front disintegration becomes more likely as the glacier approaches flotation due to thinning/stretching from steepening stress gradients caused by front acceleration. Once the glacier reaches flotation, full-depth fracture penetration is allowed because high water pressures in basal crevasses can offset the stabilizing effect of ice overburden pressure (Van der Veen, Reference Van der Veen1998; Benn and others, Reference Benn, Warren and Mottram2007; Murray and others, Reference Murray2015c; Benn and Åström, Reference Benn and Åström2018). Although smaller calving events may happen and still contribute to the discharge (e.g. Rignot and others, Reference Rignot, Fenty, Xu, Cai and Kemp2015), full-glacier thickness iceberg calving at near-floating termini provides a mechanism for rapid ice front disintegration. Tidewater glacier calving behavior exhibits a complex seasonality and includes a wide variety of calving (calving style) which is sensitive to the glacier geometry (Meier and Post, Reference Meier and Post1987; Benn and others, Reference Benn, Warren and Mottram2007). Calving style varies with changing boundary conditions. Two major calving scenarios can be identified: (1) calving of large tabular icebergs caused by longitudinal stretching at long floating ice-tongues or near the grounding line (Figs 1a and b), and (2) buoyancy-driven calving which leads to nontabular icebergs capsizing due to their small longitudinal width-to-height ratio (Fig. 1c). Capsizing calving is facilitated when the glacier terminus is lightly grounded and approaches flotation (Amundson and others, Reference Amundson2010; James and others, Reference James, Murray, Selmes, Scharrer and OLeary2014). At a floating ice-tongue, deviatoric stresses and strain rates tend to be smaller than for a grounded terminus, which likely causes slower rates of rift propagation and thereby promotes large tabular calving (Reeh, Reference Reeh1968).

Fig. 1.

Scenarios for km-scale iceberg calving at fast tidewater glaciers. At long floating ice-tongues, calving of large tabular icebergs is preferred due to smaller stress and strain rates. At grounded or near-grounded glaciers, calving of tabular (b) and nontabular (c) full-glacier thickness icebergs is enhanced when the terminus approaches flotation. Buoyancy-driven calving (c) is likely to produce icebergs with small width-to-height ratios which will capsize against the terminus front. The generated iceberg-to-terminus contact force is responsible for the production of glacial earthquakes.

Buoyancy-driven or capsizing calving appears to be a major component of current glacier dynamics in Greenland (e.g. Murray and others, Reference Murray2015c). Rapid ice flow into deepening water can create ‘super-buoyant’ conditions, in which the ice is held below buoyant equilibrium and is thus subjected to upward-directed buoyant forces (Fig. 1c). These forces can be relaxed by rapid ice viscous flow, and/or can lead to upward propagation of basal fractures and calving (James and others, Reference James, Murray, Selmes, Scharrer and OLeary2014; Wagner and others, Reference Wagner, James, Murray and Vella2016; Benn and others, Reference Benn2017a). In the meantime, the nascent iceberg starts rising upward promoting crevasse enlargement. Nontabular calving has strong implications for glacier dynamic response and stability as such events cause glacier speed-up over 4–5 days (Joughin and others, Reference Joughin2008a; Nettles and others, Reference Nettles2008; Rosenau and others, Reference Rosenau, Schwalbe, Maas, Baessler and Dietrich2013; Holland and others, Reference Holland2016).

Glacial earthquakes

GEs are magnitude M _SW ~ 5 events that are recorded globally (Ekström and others, Reference Ekström, Nettles and Abers2003). They typically are characterized by long-period (up to 150 s) surface waves propagating in the solid Earth which are produced by capsizing icebergs colliding with the glacier terminus (Amundson and others, Reference Amundson2008; Murray and others, Reference Murray2015a). As for calving behavior, there is a strong seasonality in the occurrence rate of GEs with generally higher activity in summer months when glaciers have retreated to near-grounded positions and accelerated (Joughin and others, Reference Joughin2008a; Veitch and Nettles, Reference Veitch and Nettles2012) and/or in the absence of a stiff ice-mélange layer in the proglacial fjord (Amundson and others, Reference Amundson2010; Olsen and Nettles, Reference Olsen and Nettles2017). From 1993 to 2013, 15 Greenland tidewater glaciers were observed to generate a total of 450 GEs. Comprehensive descriptions of GE patterns are addressed in Veitch and Nettles (Reference Veitch and Nettles2012) and Olsen and Nettles (Reference Olsen and Nettles2017). Since the 2000s, annual occurrence has more than quadrupled. GE productivity has expanded northward on the West coast to previously inactive glaciers. The evolution of GE production is consistent with the global acceleration of Greenland glacier flow and increasing numbers of glaciers which have started to retreat and calve at enhanced rates (Howat and others, Reference Howat, Joughin and Scambos2007; Howat and Eddy, Reference Howat and Eddy2011; Moon and others, Reference Moon, Joughin, Smith and Howat2012; Carr and others, Reference Carr, Vieli and Stokes2013; Velicogna and others, Reference Velicogna, Sutterley and Van Den Broeke2014; Murray and others, Reference Murray2015b).

GEs constitute a promising tool to track changes in individual glacier dynamics, terminus retreat, shape and buoyant conditions (Veitch and Nettles, Reference Veitch and Nettles2017; Olsen and Nettles, Reference Olsen and Nettles2017). Released seismic energy could provide insight into the amount of ice discharge. Creating a database for calved iceberg volumes and geometries from seismic records would significantly improve our understanding of calving processes and help in investigating glacier dynamic response to mass loss. GEs could contribute in an effective way to monitoring GrIS calving discharge and freshwater flux around Greenland, when image processing methods based on satellite records (e.g. Sulak and others, Reference Sulak, Sutherland, Enderlin, Stearns and Hamilton2017) are limited and labor-intensive on a routine basis.

Source mechanisms of GEs have been widely discussed as to whether they were due to glacier acceleration and slipping as observed in Antarctica (e.g. Wiens and others, Reference Wiens, Anandakrishnan, Winberry and King2008; Pratt and others, Reference Pratt, Winberry, Wiens, Anandakrishnan and Alley2014) or iceberg capsizing processes (Tsai and others, Reference Tsai, Rice and Fahnestock2008). Once fully detached, the iceberg slowly overturns and applies a horizontal contact force normal to the terminus, compatible with seismic source inversion solutions (Tsai and others, Reference Tsai, Rice and Fahnestock2008; Sergeant and others, Reference Sergeant2016; Veitch and Nettles, Reference Veitch and Nettles2017). During the capsize, temporary reverse motion of the glacier surface of cm-scale has been recorded and attributed to ice elastic deformation due to the iceberg contact forcing (Murray and others, Reference Murray2015a). The iceberg capsize perturbates deep water and lowers the water pressure beneath the floating tongue according to Murray and others (Reference Murray2015a). As a result, the floating terminus bends downward a few cms and this could be responsible for the small vertical component of the GE recorded force.

Since the work of Ekström and others (Reference Ekström, Nettles and Abers2003), GE sources have been commonly represented with a CSF model (centroid single force – a force time series described by two symmetric boxcars with opposite polarities, Fig. S1). Whereas the surface-wave magnitude M _SW is based on the amplitude of seismograms (Table 1) and constitutes a common representation of the size of ‘classic’ earthquakes, GEs are used to be characterized by the CSF source magnitude A _CSF which is found by twice-integrating in time the CSF models. Distributions of A _CSF vary from one glacier to another in Greenland (Veitch and Nettles, Reference Veitch and Nettles2012) but their relation to iceberg sizes is not straightforward. Observations (Walter and others, Reference Walter2012; Sergeant and others, Reference Sergeant2016) and modeling results (Tsai and others, Reference Tsai, Rice and Fahnestock2008; Amundson and others, Reference Amundson, Burton and Correa-Legisos2012a; Sergeant and others, Reference Sergeant2018) show that the magnitude A of the capsize-generated contact force is not an unambiguous function of the iceberg volume, but results from a complex combination of the ice-block dimensions and non-linear hydrodynamics effects. Quantification of iceberg sizes and associated discharge is then almost impossible based on the seismic magnitude A _CSF only.

Table 1.

List of quantities commonly used to characterize glacial earthquakes. Range values except for the contact force magnitude A are computed from seismic inversions and are based on the Nettles earthquake catalog (Tsai and Ekström, Reference Tsai and Ekström2007; Nettles and Ekström, Reference Nettles and Ekström2010; Veitch and Nettles, Reference Veitch and Nettles2012; Olsen and Nettles, Reference Olsen and Nettles2017). The values of A are from mechanical modeling of capsizing iceberg contact forces (Sergeant and others, Reference Sergeant2018) and are dependent on iceberg size. The filter values indicate the period corners of bandpass causal filters applied to glacial earthquake seismograms when converted to ground displacement.

Proposed method

Broadband seismic inversion accurately retrieves not only the force direction (Olsen and Nettles, Reference Olsen and Nettles2017) but also its time-evolution (force history) which captures the capsize dynamics (Sergeant and others, Reference Sergeant2016, Reference Sergeant2018). We use GE force history to recover associated iceberg volumes. To interpret the forces in terms of iceberg properties, we use horizontal contact forces computed with a two-dimensional (2-D) numerical model of iceberg capsize (Sergeant and others, Reference Sergeant2018). Numerous simulations were carried out and cataloged for different iceberg dimensions to reproduce possible contact forces. As the first step, we assume in our analysis that events are due to bottom-out capsizes only, which occur more frequently than top-out events (Amundson and others, Reference Amundson2010). GE forces are then compared with modeling results. By looking at the similarity between data and modeled force time-series, iceberg attributes are inverted when minimum misfit values are achieved. This study presents an innovative procedure for characterizing and monitoring calving volumes from GEs providing the first-of-its-kind calculation of the GrIS iceberg discharge from a seismo-mechanical coupling approach. Calculated calving volumes and related discharge are then analyzed in light of Greenland glacier dynamics over the investigated period 1993–2013.

Seismic waveform inversion and glacial earthquake forces

We use GE detections from 444 events occurring in 1993–2013 as the basis of our analysis. We combine all seismic events reported in Tsai and Ekström (Reference Tsai and Ekström2007); Nettles and Ekström (Reference Nettles and Ekström2010); Veitch and Nettles (Reference Veitch and Nettles2012) and Olsen and Nettles (Reference Olsen and Nettles2017) into one GE database that we call the ‘Nettles catalogue’. Previous studies use intermediate-period surface waves obtained from globally distributed seismometers. They invert for earthquake location (with a mean error of 15 km) and CSF source parameters (force direction and magnitude A _CSF) using a methodology similar to that routinely employed for tectonic events of similar magnitudes M _SW (Ekström, Reference Ekström2006a). Despite the good agreement between GE records and synthetic seismograms computed with the CSF solutions, GE sources may not be uniquely defined as shown by Tsai and others (Reference Tsai, Rice and Fahnestock2008) and Walter and others (Reference Walter2012) who obtain equivalent fits of GE waveforms when using different simplified source models (e.g. CSF with various source durations, simple boxcar function, Dirac delta function, etc). In order to obtain as precise GE forces as possible, we use the broadband seismic inversion method of Sergeant and others (Reference Sergeant2016) which has the advantage of requiring no restrictive a priori constraints on the force history.

We invert the forces of the Nettles catalog using seismic records from GLISN (Greenland ice-sheet Monitoring Network) broadband stations in Greenland. Because GE seismograms are polluted with a wide spectrum of signals (e.g. seismic ambient noise, high-frequency seismicity related to iceberg fracture, acceleration of the ice-mélange layer and ice-block collisions, low-frequency ground tilt and seiche phenomena, Amundson and others, Reference Amundson2010, Reference Amundson2012b; Walter and others, Reference Walter, Olivieri and Clinton2013; Sergeant and others, Reference Sergeant2016), seismic waveforms are inverted for periods of 20–100 s. Seismic inversion retrieves the force history in the North, East and vertical directions. With a polarization analysis of the force components, we can accurately compute the force direction given by its azimuth ψ (positive clockwise from the North) and vertical incidence angle when the force time-series are constantly linearly polarized. In this study, only horizontal forces are scanned and compared with mechanical simulations of iceberg-to-terminus contact forces. Azimuth angles are retrieved with a 180^° ambiguity which leads to a force that points either up-glacier or down-glacier. To ensure that inverted forces are rotated toward the glacier terminus to properly describe applied contact forces, we visually inspect individual GE locations with respect to the average orientation of local calving fronts (Veitch and Nettles, Reference Veitch and Nettles2017; Olsen and Nettles, Reference Olsen and Nettles2017) and we force the azimuth angles ψ, when aligned with the glacier flow, to be consistent with up-glacier directions (i.e. in general, 0° ≤ ψ ≤ 180° for GEs located at Western Greenland glaciers, and −180° ≤ ψ ≤ 0° for Eastern glaciers).

Out of 444 events, the inversion of only 203 GEs yield force sources with horizontal directions normal to the calving front (Fig. 2a) and vertical incidences generally <15° above the horizontal. These characteristics coincide with CSF solutions of other studies and with our understanding of iceberg capsize seismogenic source (Murray and others, Reference Murray2015a). For 241 events, we were unable to confidently calculate the force history and direction due to complex polarization patterns. A 1-D Earth's model and only Greenland stations were used here. Including other global seismometers and 3-D Earth models in the inversion scheme should give a better characterization of the source of these events.

Fig. 2.

Calved icebergs at Helheim glacier, GLISN seismic stations and geometry of the capsize model. (a) Locations of icebergs that calved on 10 July 2010, 17:49 UTC (green) and 25 July 2013, 03:13 UTC (blue), superimposed on Landsat 7 image from 10 July 2010. The iceberg surface projections are mapped from Murray and others (Reference Murray2015a) and James and others (Reference James, Murray, Selmes, Scharrer and OLeary2014). Black arrows indicate the horizontal direction of the forces inverted from glacial earthquake records at seismic stations shown in yellow triangles. (b) Schematic force balance (excluding drag forces and torques) acting on the 2-D iceberg (height H, along-glacier width εH) which tips over in bottom-out style in water with surface elevation z _B + Δz, where z _B is the water level for hydrostatic balance of the ice-block at the initiation of its capsize. Here the iceberg is super-buoyant as it experiences a positive Δz. The generated contact force F _c is computed by integrating the normal stress on the calving face and is upscaled by the iceberg across-glacier length L.

Iceberg capsize model and contact forces

The iceberg capsize model of Sergeant and others (Reference Sergeant2018) is used to compute a catalog of contact forces. It solves the 2-D dynamic equations for solid mechanics and uses the finite-element description. The model accounts for various rheologies and can be easily expanded to glacier mechanics. In this study we restrict the model set-up to a simple configuration of icebergs capsizing bottom-out (i.e. the iceberg tips over counterclockwise in Fig. 2b) against a vertical wall which represents the fixed glacier calving front. It implies that basal sliding and elastic deformation of the glacier are neglected as well as its viscous flow, which should be negligible at calving timescales.

The model is 2-D and uses box-shaped icebergs with rectangular sections of height H and along-glacier width W. To be able to capsize, the iceberg aspect ratio ε = W/H has to be lower than a critical value which depends on the ice-to-water density ratio (MacAyeal and others, Reference MacAyeal, Scambos, Hulbe and Fahnestock2003). The iceberg is partly submerged in a ‘virtual layer’ that represents the surrounding sea water defined by the water level z _w. The ice block then experiences the forces and moments exerted by gravity, buoyancy, water drag and contact reaction from the wall.

As nascent icebergs can be out of their hydrostatic equilibrium (e.g. James and others, Reference James, Murray, Selmes, Scharrer and OLeary2014; Murray and others, Reference Murray2015c), we control the water level to mimic various initial buoyancy conditions. The difference in water surface elevation experienced by the iceberg (z _w) with respect to the theoretical one for an initially neutrally buoyant iceberg of the same geometry (z _B) is denoted Δz. Icebergs which experience a negative Δz are then uplifted by buoyancy compared with their buoyant position at the initiation of the capsize. Conversely, for positive Δz, icebergs sink deeper.

The contact force generated on the terminus by the rotating iceberg is recovered by integrating the solid mechanics equation in time using the method of Hilber and others (Reference Hilber, Hughes and Taylor1977). A catalog of contact force historiesFootnote ¹ is then computed with varying sets of (ε, H, Δz) values. The force generated by a 3-D ice block is computed using a scaling factor equal to the iceberg across-glacier length L. Possible contact forces are explored for iceberg geometries which describe typical Greenland glaciers. Iceberg volumes covered by the computed contact force catalog range between 0.018 and 4 km³ corresponding to parameters: 0.1 ≤ ε ≤ 0.6, 600 m ≤ H ≤ 1100 m, 500 m ≤ L ≤ 6000 m and − 20 m ≤ Δz ≤ 20 m.

Synthetic inversion of iceberg volumes and induced errors

We first run synthetic inversions to estimate the errors we could obtain on the iceberg volume with the first step of this procedure. We run several tests taking as input 20–100 s filtered contact forces which correspond to icebergs of various dimensions ε and H but one given length L and at hydrostatic balance (i.e. Δz = 0). We also investigate the dependency of the inversion results on the data force's signal-to-noise ratio (SNR). For one SNR, we average the errors on calculated iceberg volumes over 100 tests, each one taking as input a force that has been polluted by random noise which has the same dominant frequency.

For noise-free inversions, maximum VR is always reached for the exact input aspect ratio and height which then give correct volumes V ₀ (Fig. 3a). When using noisy forces, maximum volume errors (15% for SNR=0.5 and 4% for SNR=5) are reached for the cases of intermediate aspect ratios. Depending on the tested geometry and notably the input aspect ratio, VR values decrease more or less fast around the maximum peak (Figs S2c and d). For each test, the volume uncertainty interval δV ₀ plotted in Fig. 3b is constructed by taking all (ε, H) combinations that yield VR values >98% of maximum VR. While δV ₀, as defined here, is always <15% for icebergs with ε ≤ 0.2, δV ₀ increases up to 40% for intermediate ε due to wider distributions of high VR values. The variation of the contact force with iceberg dimensions is approximately symmetric about ε ≈ 0.4 (see Sergeant and others, Reference Sergeant2018, Fig. 6). The confidence interval on ε and H and therefore volume will then be higher and more difficult to reduce when the iceberg has an intermediate aspect ratio close to 0.4, as rather similar contact forces are observed around this value. The V ₀ error or uncertainty δV ₀ are then linked to the ability to discriminate contact forces with the aspect ratio, which can be even more difficult for poorly inverted data forces. Given the synthetic tests with the force SNR, in the worst case scenario, the maximum volume error would be 0.3 km³ for the biggest icebergs that a Greenland glacier is able to produce.

Fig. 3.

Error on (a) estimated iceberg volumes and (b) associated confidence intervals from synthetic inversions run with 20–100 s filtered contact force histories for any iceberg of across-glacier length L=1. Synthetic tests are for noise-free input forces, and forces that have been polluted by random noise given a signal-to-noise ratio (SNR).

Validation study on two calving events at Helheim glacier

We illustrate the inversion process for two GEs that occurred at Helheim glacier. The associated calving events were captured on cameras, making it possible to estimate the iceberg dimensions (Fig. 2a). This section serves as a benchmark for our iceberg geometry calculations. One event (25 July 2013, 03:13 UTC) results from the capsize of an iceberg with a height equal to the full-glacier thickness H _d ~ 800 m, an across-glacier length L _d ~ 2500 m, and an aspect ratio ε_d ~ 0.23 (Murray and others, Reference Murray2015a). The iceberg volume was then estimated as 0.37 km³. The second event (12 July 2010, 17:49 UTC) results from the capsize of a 1.18 km³ iceberg with H _d ~ 800 m, ε_d ~ 0.37 and 4000 m < L _d < 5000 m (James and others, Reference James, Murray, Selmes, Scharrer and OLeary2014). The GE horizontal force components shown in Fig. 4a (blue and green lines) have maximum amplitudes of 3 × 10¹⁰ N and 4.6 × 10¹⁰ N for the 2013 and 2010 events, respectively.

Fig. 4.

Iceberg volume inversion for two well documented calving events of dimensions (ε_d, H _d, L _d). (a) Horizontal component of the force inverted from seismic data (blue or green for the 2013 and 2010 events, respectively) and associated best-fitting force models (red) bandpass filtered in the GE seismic band (20–100 s). Best force models correspond to parameters (ε₀, H ₀, Δz ₀, L _d) with the iceberg across-glacier length L _d derived from field observations. Red dashed lines show the corresponding actual (i.e. unfiltered) contact force F _c generated by the iceberg on the terminus. Grey-to-black shaded lines represent all force models that equivalently fit the seismic force with high variance reduction (VR) values. Yellow box indicates the iceberg-to-terminus contact duration. Snapshots of numerical simulations that are illustrated in (b). (c) Misfit function with iceberg dimensions ε and H for the catalog of contact forces scaled by the actual iceberg length (L = L _d). For this representation, the third model parameter Δz has been fixed to the value reached for maximum VR, i.e. Δz ₀ indicated above each panel. The real and calculated iceberg dimensions are (ε_d, H _d) and (ε₀, H ₀), indicated by white and yellow-filled circles, respectively. The white lines show isovolumetric contours, when L is kept constant. (d) Computed iceberg volumes when varying L-values. White dashed lines indicate the actual iceberg length L _d and volume V _d. Uncertainties on estimated volumes are calculated from all (ε, H) combinations that yield VR values>98% of maximum VR, i.e. within the domain indicated by black dashed contour lines in (c). Derived iceberg dimensions (ε₀, H ₀) are indicated for some L-values. Beyond a critical iceberg length L, inverted aspect ratios and volumes are stable around correct values. The red thick line indicates the average calculated volume for the indicated L-range, i.e. when ε₀ is stable and H ₀ does not exceed ± 15% of the glacier thickness H _d.

We first show in Fig. 4c the VR-misfit function with ε and H, obtained for L fixed to the actual iceberg length L _d. Yellow-filled and white circles indicate best-fitting parameters (ε₀, H ₀) and true iceberg dimensions (ε_d, H _d), respectively. For the 2013 event, the best VR value (0.95) is achieved for ε₀ = 0.2, H ₀ = 850 m and Δz ₀ = −9 m. A negative Δz corresponds to an uplift of the ice block prior to its capsize which is well justified by the observations of Murray and others (Reference Murray2015c). For the 2010 event, maximum VR (0.98) corresponds to a 11 m-uplifted iceberg with ε₀ = 0.33 and H ₀ = 825 m. For both events, minimization-determined parameters are in good agreement with the true iceberg geometry and yield iceberg volumes V ₀ of ~ 0.36 km³ and ~ 1.12 km³. Computed volumes are underestimated by <5%. While the volume uncertainty (illustrated by the black dashed line in Fig. 4c) is quite small (δV ₀ < 7%) for the thinner iceberg, it increases up to 40% of V ₀ for the larger calving event. Retrieved δV ₀ intervals follow the same trend as observed with synthetic inversions with higher uncertainties obtained for intermediate aspect ratios. Finally, best-determined force models (red in Fig. 4a) reproduce the force data components well during the time period coincident to the capsize and indicated by the yellow box (Murray and others, Reference Murray2015a). Force models determined within the confidence interval (black-shaded lines) are very similar to the best-fitting force. They still show different force amplitudes which can also encompass amplitude uncertainties on the data forces that were obtained from GE records (Sergeant and others, Reference Sergeant2016).

The results of all inversion run with other L values are reported in Fig. 4d. Best retrieved parameters (ε₀, H ₀) and volumes V ₀ vary with L, as L impacts the force amplitude to be matched. In a certain L range (indicated by red thick lines), the calculated aspect ratios and volumes are consistent and match the field observations, independently of the iceberg length. While the aspect ratio essentially controls the force time evolution and then gives consistent results over the L-fixed inversions, the force amplitude is more affected by the other iceberg dimensions as it is approximately proportional to H ³ and L (Tsai and others, Reference Tsai, Rice and Fahnestock2008; Amundson and others, Reference Amundson, Burton and Correa-Legisos2012a; Sergeant and others, Reference Sergeant2018). This means that for a given force, considering a well resolved ε as it is, a specific combination of decreased H and increased L values will result in the same iceberg volumes.

This analysis shows that, unless at least one dimension of the iceberg is known (as here L), it is difficult to invert all three individual iceberg dimensions (and then the iceberg geometry) from GE horizontal forces. However a volume can be robustly estimated when consistent results are obtained on a loop of inversions.

General inversion procedure for glacial earthquakes

The above procedure is applied to the forces computed for 203 GEs. Final calving volumes V are calculated as the average of all $V_0=\epsilon _0H_0^2L_{\rm i}$ that are estimated for a given iceberg length L _i, when accepted. This means that V ₀, to be included in the final calculation of V, must (1) be stable around a certain value (i.e. V ₀(L) function must have reached a plateau), while (2) ε₀ is also stabilized and (3) H ₀ does not exceed the actual glacier thickness ± 15%. In this study, we consider glacier thicknesses estimated in previous studies (e.g. Joughin and others, Reference Joughin2008a,Reference Joughinb, Reference Joughin, Smith, Shean and Floricioiu2014; Stearns and Hamilton, Reference Stearns and Hamilton2007; Howat and others, Reference Howat2011; Murray and others, Reference Murray2015c) or fixed to 800 m for most events. The latter condition (3) on the iceberg height is, in the end, not too restrictive as a maximum deviation of H ₀ values of 15% of 800 m would nearly cover the entire range of investigated H given considered L range. The constraint on H and then V could be improved in future work by using bedmap data for Greenland (Griggs and others, Reference Griggs2012) or recent ice thickness measurements from IceBridge surveys (Rignot and others, Reference Rignot, Mouginot, Larsen, Gim and Kirchner2013). However the inversion results are consistent both with or without constraining H (Fig. S3). The final volume uncertainty δV for one event is computed as the quadrature propagated uncertainty

(1)$$\delta V = \displaystyle{1 \over N}\sqrt {\sum\limits_i^N {\delta V_0{(L_{\rm i})}^2} } $$

where δV ₀(L _i) is the uncertainty related to each of the N accepted models with given L _i. We find δV uncertainties that vary between 0.01 and 0.33 km³ with 70% of the events reaching δV < 0.1 km³ (Fig. S4). The median value of δV is 0.07 km³ across all inverted GEs. This represents an average volume uncertainty of 6%.

Earthquake sizes and calving volumes

As noticed by Veitch and Nettles (Reference Veitch and Nettles2012), Greenland glaciers have different shapes of GE size distribution as measured by CSF. While every glacier's seismic activity shows a characteristic peak in magnitudes with a rapid decline at larger and smaller sizes of GEs, the peak occurs at different sizes depending on the source region (Fig. 5). The distribution width also varies from one place to another. Previous studies argue that the shape of GE size-frequency distributions is controlled by the glacier geometry with the upper bound of the magnitude being related to the biggest icebergs whose dimensions cannot exceed actual glacier thickness and width (Nettles and Ekström, Reference Nettles and Ekström2010; Veitch and Nettles, Reference Veitch and Nettles2012). As shown in Figs 5 and 6a–h, the biggest icebergs are observed at the largest glaciers Jakobshavn Isbrae, Kangedluqssuaq and Helheim where highest GE magnitudes A _CSF are also recorded. The iceberg size-frequency distributions we obtain vary from one glacier to another and somehow follow the same evolution as for GE magnitudes, but not systematically at every calving site (Fig. 5). For example, the calving volume distribution at Kangerdlussuaq is wider than at Helheim as is the case for GE sizes, with a volume peak that is larger (1.5 km³) than Helheim's peak (0.7 km³). On the contrary, calving size distribution at Jakobshavn Isbrae shows a peak at 0.9 km³, which is larger than at Helheim, despite a smaller seismic magnitude peak. Even if all inverted calving volumes are actually below the maximum size of icebergs that each glacier is able to produce given the glacier width, the largest icebergs do not systematically correspond to the highest GE magnitudes as illustrated below.

Fig. 6.

Earthquake force magnitudes A _CSF (squares) and associated contact force magnitudes A (dots) as a function of inverted calving volumes V at 8 Greenland glaciers (a-h). A _CSF are obtained by seismic waveform inversion (Tsai and Ekström, Reference Tsai and Ekström2007; Veitch and Nettles, Reference Veitch and Nettles2012; Olsen and Nettles, Reference Olsen and Nettles2017). A are the magnitudes of actual contact forces (i.e. non-filtered) obtained when with our volume inversion method. Color indicates the inverted iceberg aspect ratio ε. Gray-shaded boxes denote an estimate of the range of possible iceberg size and corresponding force magnitudes A based on the geometry of the glacier terminus. For a comparison, synthetic evolutions of A _CSF and A with iceberg volume and ε are presented in (i), as well as the evolution of the magnitude of the contact force A _c ∝ H (1 − ε) V (Eqn (3), color lines) when hydrodynamic effects are neglected and H is fixed to 900 m.

Figure 6 shows the evolution of GE magnitudes A _CSF (squares) with inverted volumes V , and associated actual contact force magnitudes A (dots, actually refer to non-filtered contact force models). Color indicates the iceberg aspect ratio ε calculated for each event. Comparing the dependence of A and A _CSF on computed V at individual glaciers (Figs 6a–h) with the synthetic tests run with different iceberg sizes (Fig. 6i), the calving volumes we obtain for GEs yield reasonably good expected GE magnitudes. From this representation, we see that equivalent GE magnitudes A _CSF can be reached for different sizes of icebergs and aspect ratios. These results tell us that the GE size-frequency distribution should not be interpreted in terms of exact calving volume distribution as GE magnitudes are dependent on the iceberg size, aspect ratio and height, collectively.

Mass loss related to glacial earthquakes and contribution to dynamic discharge

We now investigate the temporal evolution of the discharge associated with GEs over 1993–2013. To extrapolate the calving volume inversion results to the entire dataset of GEs which could not be characterized, we simply add to the iceberg catalog the most common ice volume that was calculated at each glacier (Table 2, column 5), when one seismic event was missed. This enables us to build the projected iceberg discharge we could expect if we had inverted the calving volumes for all GEs referenced in the Nettles catalog. Lower and higher expectations are built based on the smallest and biggest iceberg volumes that were inverted at each glacier. Results for individual glaciers and averaged over Greenland regions are shown in Figs 12–14 and Fig. 7, respectively. We will only discuss here the projected GE-associated discharge (called GED in the following). Annual and cumulated GED over the two decades are shown by blue bars and the blue line in Figs 7d–f.

Fig. 7.

Glacial earthquake (GE) and iceberg production across Greenland in 1993–2013. Evolution of (a–c) number of GEs, (d–f) associated calving volumes, (g–i) GE-associated discharge (GED, blue line) and dynamic discharge (DD, black line) at investigated glaciers and (j–l) GED contribution to DD in 2000--12. In (a–c) and (d–f), orange bars represent the number of GEs that were inverted in this study, and associated volumes; blue bars represent all GEs from the Nettles catalog and associated volumes that are expected based on the iceberg size-frequency distribution inverted at each glacier for the dataset in orange (Fig. 6). Right y-axis correspond to (a–c) the proportion of inverted GEs with respect to the Nettles catalog in the region and updated each year, and (d–f) cumulated GED based on inverted events only (yellow line), and based on expectations for all GEs (blue line) with their uncertainty (gray area). In (g–i), gray line (right y-axis) show dynamic discharge DD in the region when including also tidewater glaciers which did not produce GEs. Only eight GE-producing glaciers are accounted for in this study (Table 2). Details about GE production at individual glaciers are provided in Figs 12–14.

From our calculations, the eight glaciers included in this study produced between 270 km³ and 700 km³ with a median value of 398 km³ of ice (367 Gt equivalent) via GEs over the two decades 1993–2013. The biggest producers are Jakobshavn (71 km³), Helheim (86 km³) and Kangerdluqssuaq (117 km³). Altogether they account for half of the total GED we extrapolated for Greenland. While Eastern glaciers (Kangerdluqssuaq and Helheim) cumulated GED appears to have increased constantly over the whole time period, Western glacier mass loss has doubled since 2010 and thus contributes to the acceleration of Greenland GED since then. Before 2009, Western glacier annual GED never exceeded the average annual GED in East Greenland (i.e. 10 km³) and contributed to only 25% of the overall Greenland GED. In 2013, West Greenland GED exceeded the calving volumes produced by both Kangerdluqssuaq and Helheim cumulated over 2004/05 (i.e. 32 km³) when the two Eastern glaciers experienced a synchronous and fast retreat of more than 5 km, reaching their most retreated position ever recorded (Howat and others, Reference Howat, Joughin, Tulaczyk and Gogineni2005, Reference Howat, Joughin, Fahnestock, Smith and Scambos2008; Joughin and others, Reference Joughin2008a; Bevan and others, Reference Bevan, Luckman and Murray2012). Presently, the Greenland GED cumulated over 1993–2013 seems to have been equally partitioned between East and West regions (Table 2, column 6). If trends continue, the early dominance of Eastern glaciers is about to change notably due to a decrease in GE activity at Kangerdluqssuaq, which was the most productive Eastern glacier before its retreat, and most importantly to increasing GED contributions from Western glaciers including Jakobshavn and North West glaciers.

One concern is the contribution of buoyancy-driven calving to Greenland dynamic discharge (called DD in the following) which also includes submarine melting. Measurements from satellite imagery allow linear or areal frontal ablation rates to be determined, but it is rarely possible to measure the relative contribution of calving and melting below the surface (Benn and others, Reference Benn, Cowton, Todd and Luckman2017b). Because of these difficulties and the importance of glacier flow in governing rates of ice delivery to the terminus, DD is commonly quantified using measurements of the ice volume flux passing through a glacier cross-section. To estimate the proportion of GED to DD, we use DD measurements of Enderlin and others (Reference Enderlin2014) for years 2000–12. The GED generally follows the trend of DD (blue line versus black line in Figs 7g–i), except after 2010 when an acceleration of GED occurred at Western glaciers while the dynamic discharge observed in this area continued to increase constantly. We observe different levels of GE contribution to DD through time but also depending on the glacier and on the region (Table 2, column 7; Figs 7j–l). In 2012, Western glaciers experienced up to 25% of mass loss via GEs while the amount of GED tended to stabilize or even decrease when averaged over both Eastern glaciers, ~18%. In total, capsizing calving which produced earthquakes could have contributed between 8 and 21% of the discharge averaged over the Greenland glaciers included in our analysis, given the uncertainties on the projected calving volume time series. Including in our calculation every tidewater glacier in Greenland which is not observed to produce GEs (gray line in Fig. 7g), GED contribution to the GrIS discharge should be 3--4 times smaller, i.e. ~3--6% of Greenland mass loss (Table 3), although it is expected to account for more if trends continue.

Table 3.

Contribution of glacial earthquake associated discharge (GED) to Greenland glacier dynamic discharge (DD) in 2000, 2005, 2012 and over 2000–12. Results (lower-upper bounds, median value is within brackets) in lines 1–2 are when (1) considering only the eight earthquake-producing glaciers for the calculation of DD, (2) when DD measurements account for every tidewater glacier in Greenland. Dynamic discharge data are from Enderlin and others (Reference Enderlin2014).

Spatio-temporal variability of calving activity

Results are summarized in Fig. 8 and averaged over four time periods. Cumulated calving volumes and their uncertainty are represented by the size of the black circles. Green dashed circles show GED contribution to the dynamic discharge DD at each glacier. Since the evolution of annual GE associated discharge is related to earthquake occurrences at one glacier, it is not possible with the former representation to assess whether the associated mass loss is accelerating through more frequent events or through larger volumes per event, or both. We therefore introduce the cumulative GE-associated discharge anomaly C defined as

(2)$$C=\sum_k^N (V_k - \lt V\gt)$$

where V _k denotes iceberg volumes associated with N GEs within a given time period (e.g., 1 year) at one glacier, and < V > is the mean iceberg volume, averaged over all events produced at the glacier between 1993 and 2013. For constant discharge, C ≈ 0. C>0 (red colors in Fig. 8) when bigger icebergs are produced over one time period, on average. Conversely, C<0 (blue colors) when produced icebergs were smaller than usual. C results are similar whether they are computed on projected iceberg volume time series or on volumes which were actually inverted from GE records. This supports the argument that our calculations of discharge expected for the complete Nettles catalog did not introduce any bias in the analysis that follows.

Fig. 8.

Evolution of calving volumes with lower and upper uncertainties (black circles) that were produced at individual glaciers in Greenland through glacial earthquakes (GEs) for four time periods. Inner color of circles indicate value of the cumulative GE-associated discharge anomaly C at the considered time period. Red colors (C > 0) indicate that single glaciers discharge bigger icebergs with respect to their average size calculated over 1993–2013, and blue colors (C < 0) indicate that GE-associated discharge happens through more but smaller individual events. Conversely, white color (C ~ 0) indicates when glaciers produce icebergs of approximately constant sizes, compared to the average size over the whole time period. GE contributions to the glacier dynamic discharge DD in each period of time are indicated by the size of green circles. Iceberg production is only represented at glaciers where enough GEs (> 15% of the GE production) could be inverted.

Eastern glaciers

On the contrary, GE production globally decreased in South Eastern Greenland after 2005 as Helheim and Kangerdluqssuaq stabilized after retreat (Bevan and others, Reference Bevan, Luckman and Murray2012). Nevertheless, differences occur between the two glaciers. While GE-associated discharge has been sustained at Helheim with successive periods of bigger and smaller iceberg release (Fig. 9b), Kangerdluqssuaq production in terms of events, iceberg sizes and GED contribution to discharge DD has been decreasing since 2006 and particularly 2011 (Fig. 9d). Since 2005, GEs have been active at Kangerdluqssuaq only during retreat phases when the glacier retreated toward specific positions in the fjord and accelerated (Fig. 9c), suggesting that terminus position and glacier speed play an essential role in capsizing calving. By 2011 as a consequence of sustained thinning (Howat and Eddy, Reference Howat and Eddy2011), Kangerlussuaq's grounding line had retreated into shallower water and remained stable as the glacier developed a 5 km-long floating ice-tongue, affecting its calving behavior (Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017). Since 2011, we observe smaller calving events at Kangerdluqssuaq which contrasts with the 2005–13 calving behavior at Helheim which had maintained a lightly grounded terminus (James and others, Reference James, Murray, Selmes, Scharrer and OLeary2014; Kehrl and others, Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017).

Fig. 9.

1993–2013 front positions (orange) and speed (black) at (a) Helheim, (c) Kangerdluqssuaq and (e) Jakobshavn. Vertical bars indicate timing of glacial earthquakes (GEs). (b, d, f) Evolution of the cumulative GE-associated discharge anomaly C at corresponding glacier. 1993–2011 Eastern glacier positions are from Bevan and others (Reference Bevan, Luckman and Murray2012) and updated to 2013. 2009–13 Jakobshavn pos: itions are from Joughin and others (Reference Joughin, Smith, Shean and Floricioiu2014).

Links to glacier dynamics

To link calving patterns to ice flow dynamics, we analyze the evolution of calculated calving volumes with glacier speed and front position at Jakobshavn and Helheim. Jakobshavn Isbrae, while retreating, has continuously produced GEs since 2005 (Figs 9e and f) after the disintegration of a 15 km-long floating tongue (Joughin and others, Reference Joughin, Abdalati and Fahnestock2004). Between 2009 and 2013, the annual number of seismic events and related discharge has almost tripled and as previously reported, the average trend has produced bigger calving events. If we look at the dependency of GE production and calving volumes on ice flow speed and glacier geometry (Fig. 10), we observe that the largest icebergs calve at maximum speeds (i.e. 15 km a⁻¹) and at specific terminus positions (i.e. 10–11 km relative to the reference point in the along-flow profile in Joughin and others, Reference Joughin, Smith, Shean and Floricioiu2014). In 2009–13, Jakobshavn remained near-grounded as the calving front was resting over a bed depression (Joughin and others, Reference Joughin, Smith, Shean and Floricioiu2014). In July--September 2012/13 (orange and dark red data points), the glacier reached its most retreated position at greatest water depths coinciding with highest flow speeds, largest icebergs and also fewest calving events. As a terminus with little or no floating extension advances into deep water and approaches buoyancy, hydrostatic imbalance at the front causes downward bending (Wagner and others, Reference Wagner, James, Murray and Vella2016). This initially prevents the ice from attaining buoyant equilibrium and buoyant flexure of the terminus must be balanced upstream by enhanced longitudinal stress. This can lead to glacier stretching and acceleration (Nick and others, Reference Nick, Vieli, Howat and Joughin2009; Joughin and others, Reference Joughin2012, Reference Joughin, Smith, Shean and Floricioiu2014) and the development of a super-buoyant tongue whose length increases as ice advances (Benn and others, Reference Benn2017a). It can become several hundreds of meters long before buoyant forces are high enough to initiate calving. Our observations of increased calving sizes with deeper water and increased flow velocity agree with modeling results for glacier retreat behavior over a reverse bed slope (Nick and others, Reference Nick, Vieli, Howat and Joughin2009) or glacier advance in deepening water (Wagner and others, Reference Wagner, James, Murray and Vella2016). A similar analysis can be done for Helheim calving production (Fig. 11): as the glacier was retreating over a reverse bed slope in 2004/05, it reached its most retreated position in deeper water in 2005 (green data points) producing bigger iceberg sizes; since that time it has alternated advance and retreat phases over the bed overdeepening producing a wide range of calving events.

Fig. 10.

Evolution of 2009–13 GE and calving production with Jakobshavn's geometry (i.e. along the southern ice-stream). Number of detected GEs (bars) and associated calving volumes (color squares) as a function of terminus (a) speed and (b) front position along the profile described in Joughin and others (Reference Joughin, Smith, Shean and Floricioiu2014). Color indicates relative time since1 January 2009. In (c) color circles indicate front positions relatively to time. Filled circles indicate GE-active periods (i.e. mostly in summer). Bed topography is in black. Blue dashed line indicates the flotation threshold, above which the ice should be grounded and below which the glacier is super-buoyant (i.e. deeper than buoyant equilibrium). Glacier surface elevation profiles are plotted for February 2009 and March 2012. Gray box indicates the magnified region in (b). Glacier geometry-related data are from Joughin and others (Reference Joughin, Smith, Shean and Floricioiu2014).

Fig. 11.

Same as in Fig. 10 for 1996–2013 calving production at Helheim glacier. 1996–2011 calving front positions are from Bevan and others (Reference Bevan, Luckman and Murray2012) and updated to 2013. Glacier geometry in (c) is from Kehrl and others (Reference Kehrl, Joughin, Shean, Floricioiu and Krieger2017).

Fig. 12.

GE and iceberg production at Helheim and Kangerdluqssuaq (South East Greenland) in 1993–2013. Evolution of (a) number of GEs, (b) associated calving volumes, (c) cumulative GE-associated discharge anomaly C, (d) GE-associated discharge (GED, blue line) and dynamic discharge (DD, black line) from Enderlin and others (Reference Enderlin2014), (e) GED contribution to DD. In (a) and (b), orange bars represent the number of GEs that were inverted in this study, and associated volumes; blue bars represent all GEs from the Nettles catalog and associated volumes that are expected based on the iceberg size-frequency distribution inverted at each glacier for the dataset in orange (Fig. 6). Right y-axis correspond to (a) the proportion of inverted GEs with respect to the Nettles catalog in the region and updated each year, and (b) cumulated GED based on inverted events only (yellow line), and based on expectations for all GEs (blue line) with their uncertainty (gray area).

Fig. 13.

Same as in Fig. 12 for Jakobshavn Isbrae, Rinks glacier and Upernavik Isström (West Greenland).

Fig. 14.

Same as in Fig. 12 for Alison, Sverdup and Kong Oscar glaciers (North West Greenland).

GEs and calving patterns depend on the terminus front position which determines the flow speed and flotation level (Nick and others, Reference Nick, Vieli, Howat and Joughin2009; Joughin and others, Reference Joughin, Smith, Shean and Floricioiu2014). The tendency of a glacier to produce full-glacier thickness capsizing icebergs is related to the exploitation of existing basal crevasses which tend to develop in the vicinity of the grounding line (Van der Veen, Reference Van der Veen1998). Murray and others (Reference Murray2015c) suggest that the propagation of crevasses to the glacier surface is strongly enhanced when the glacier base is deeper than buoyant equilibrium. This condition was met at Jakobshavn Isbrae in 2009–13 and at Helheim since 2004/05, as both glaciers were retreating and advancing over a basal depression. Since 2005, most northwestern glacier retreat has been accelerating (McFadden and others, Reference McFadden, Howat, Joughin, Smith and Ahn2011; Carr and others, Reference Carr, Vieli and Stokes2013; Moon and others, Reference Moon, Joughin and Smith2015) followed by increasing occurrences of GEs and calving volumes since 2010. Such changes in calving behavior may then be related to fresh new boundary conditions as the termini have to reach near-grounded positions and could approach near-floating and deep water conditions.