# Monitoring Greenland ice sheet buoyancy-driven calving discharge using glacial earthquakes

Amandine Sergeant, Anne Mangeney, Vladislav A. Yastrebov, Fabian S.A.

Walter, Jean Paul Montagner, Olivier Castelnau, Eléonore Stutzmann,

Pauline Bonnet, Velotioana Jean-Luc Ralaiarisoa, Suzanne Bevan, et al.

► To cite this version:

Amandine Sergeant, Anne Mangeney, Vladislav A. Yastrebov, Fabian S.A. Walter, Jean Paul Montagner, et al.. Monitoring Greenland ice sheet buoyancy-driven calving discharge using glacial earthquakes. Annals of Glaciology, 2019, pp.1-21. 10.1017/aog.2019.7. hal-02139311v2

HAL Id: hal-02139311

https://hal.science/hal-02139311v2

Submitted on 25 Jun 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Monitoring Greenland ice sheet buoyancy-driven calving discharge using glacial earthquakes

Amandine SERGEANT,<sup>1,2</sup> Anne MANGENEY,<sup>2,3</sup> Vladislav A. YASTREBOV,<sup>4</sup>

Fabian WALTER,<sup>1</sup> Jean-Paul MONTAGNER,<sup>2</sup> Olivier CASTELNAU,<sup>5</sup>

Eléonore STUTZMANN,<sup>2</sup> Pauline BONNET,<sup>2,4,5</sup> Velotioana Jean-Luc RALAIARISOA,<sup>2*</sup>

Suzanne BEVAN,<sup>6</sup> ◐ Adrian LUCKMAN<sup>6</sup>

<sup>1</sup>Laboratory of Hydraulics, Hydrology and Glaciology, ETH Zürich, Zürich, Switzerland.

E-mail: sergeant@vaw.baug.ethz.ch

<sup>2</sup>Institut de Physique du Globe de Paris, CNRS UMR 7154, Université Paris Diderot-Paris 7, Paris, France

<sup>3</sup>ANGE team, INRIA, Laboratoire Jacques-Louis Lions, Paris, France

<sup>4</sup>MINES ParisTech, PSL Research University, Centre des Matériaux, CNRS UMR 7633, Evry, France

<sup>5</sup>Laboratoire PIMM, Arts et Métiers, CNRS, CNAM, HESAM Université, Paris, France

<sup>6</sup>Geography Department, College of Science, Swansea University, SA2 8PP, UK

**ABSTRACT.** Since the 2000s, Greenland ice sheet mass loss has been accelerating, followed by increasing numbers of glacial earthquakes (GEs) at near-grounded glaciers. GEs are caused by calving of km-scale icebergs which capsize against the terminus. Seismic record inversion allows a reconstruction of the history of GE sources which captures capsize dynamics through iceberg-to-terminus contact. When compared with a catalog of contact forces from an iceberg capsize model, seismic force history accurately computes calving volumes while the earthquake magnitude fails to uniquely characterize iceberg size, giving errors up to 1 km³. Calving determined from GEs recorded at eight glaciers in 1993–2013 accounts for up to 21% of the associated discharge and 6% of the Greenland mass loss. The proportion of discharge attributed to capsizing calving may be underestimated by at least 10% as numerous events could not be identified by standard seismic detections (Olsen and Nettles, 2018). While calving production tends to stabilize in East Greenland, Western glaciers have released more and larger icebergs since 2010 and have become major contributors to Greenland dynamic discharge. Production of GEs and calving behavior are controlled by glacier geometry with bigger icebergs being produced when the terminus advances in deepening water. We illustrate how GEs can help in partitioning and monitoring Greenland mass loss and characterizing capsize dynamics.

**KEYWORDS:** climate change, glacier calving, glacier discharge, iceberg calving, seismology

## INTRODUCTION AND BACKGROUND

The Greenland ice sheet (GrIS) mass loss and its contribution to sea-level rise has more than quadrupled from 1991–2001 to 2002–11 (Shepherd and others, 2012) due to both increased discharge of ice to the ocean and decreased surface mass balance (Van den Broeke and others, 2009, 2016; Enderlin and others, 2014; Velicogna and others, 2014; Khan and others, 2014). The partitioning of mass loss between dynamic (i.e. changes in ice flow, thinning and calving rates) and surface processes is important because these losses indicate different forcing (Howat and others, 2011). While surface mass balance is driven by atmospheric processes, tidewater glacier dynamics are driven by changes in resistive stress at the terminus (e.g. Nick and others, 2009), likely triggered by changing ocean conditions (Holland and others, 2008; Straneo and others, 2010), and further modulated by glacier geometry (Howat and others, 2005, 2007; Joughin and others, 2008b, 2014; Enderlin and others,

2013; Moon and others, 2014; Felikson and others, 2017; Kehrl and others, 2017).

Between one third and one half of the annual GrIS (see Appendix A for definition of acronyms) mass is lost at the termini of tidewater glaciers (Van den Broeke and others, 2009; Enderlin and others, 2014). Iceberg calving and sub-marine melting are collectively known as dynamic dis-charge. They are notoriously difficult to quantify separately (Benn and others, 2017b). One goal of this study is to provide independent measurements of the GrIS discharge related to the calving of large rolling icebergs from emitted seismic signals, known as glacial earthquakes (GEs). The characterization of icebergs is important as calving rates are projected to increase during the coming decades and play a significant role in the stability of the entire GrIS (Rignot and Kanagaratnam, 2006; Rignot and others, 2011). Icebergs carry and release freshwater far from the glacier as they drift offshore (Enderlin and others, 2016; Wagner and others, 2017) and can potentially affect the large-scale ocean overturning circulation and ocean temperature (Fichefet and others, 2003; Holland and others, 2008; Wilton and others, 2015; Yang and others, 2016; Stern and others, 2016).

*Present address: Institut de Physique de Rennes, Université Rennes, CNRS UMR 6251, Rennes, France.
## Calving processes

Calving and related dynamic processes are poorly represented in the current generation of ice-sheet models (Benn and others, 2017a). Owing in large part to a lack of observational data (e.g. James and others, 2014), 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, 1998; Benn and others, 2007; Murray and others, 2015c; Benn and Åström, 2018). Although smaller calving events may happen and still contribute to the discharge (e.g. Rignot and others, 2015), 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, 1987; Benn and others, 2007). 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, 2010; James and others, 2014). 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, 1968).

Buoyancy-driven or capsizing calving appears to be a major component of current glacier dynamics in Greenland (e.g. Murray and others, 2015c). 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, 2014; Wagner and others, 2016; Benn and others, 2017a). 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, 2008a; Nettles and others, 2008; Rosenau and others, 2013; Holland and others, 2016).

## Glacial earthquakes

GEs are magnitude $M_{SW} \sim 5$ events that are recorded globally (Ekström and others, 2003). 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, 2008; Murray and others, 2015a). 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

<figure><img src="image_2.png" /><figcaption>**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.</figcaption></figure>

(Joughin and others, 2008a; Veitch and Nettles, 2012) and/or in the absence of a stiff ice-mélange layer in the proglacial fjord (Amundson and others, 2010; Olsen and Nettles, 2017). 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 (2012) and Olsen and Nettles (2017). 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, 2007; Howat and Eddy, 2011; Moon and others, 2012; Carr and others, 2013; Velicogna and others, 2014; Murray and others, 2015b).
GEs constitute a promising tool to track changes in individual glacier dynamics, terminus retreat, shape and buoyant conditions (Veitch and Nettles, 2017; Olsen and Nettles, 2017). 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, 2017) 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, 2008; Pratt and others, 2014) or iceberg capsizing processes (Tsai and others, 2008). 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, 2008; Sergeant and others, 2016; Veitch and Nettles, 2017). 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, 2015a). The iceberg capsize perturbates deep water and lowers the water pressure beneath the floating tongue according to Murray and others (2015a). 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 (2003), GE sources have been commonly represented with a CSF model (cendroid 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, 2012) but their relation to iceberg sizes is not straightforward. Observations (Walter and others, 2012; Sergeant and others, 2016) and modeling results (Tsai and others, 2008; Amundson and others, 2012a; Sergeant and others, 2018) 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.

## Proposed method

Broadband seismic inversion accurately retrieves not only the force direction (Olsen and Nettles, 2017) but also its time-evolution (force history) which captures the capsize dynamics (Sergeant and others, 2016, 2018). 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, 2018). 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, 2010). 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.

## COMPILATION OF FORCE CATALOGUES USED FOR CALVING VOLUME INVERSION

### 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 (2007); Nettles and Ekström (2010); Veitch and Nettles (2012) and Olsen and Nettles (2017) 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

<table><caption>**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, 2007; Nettles and Ekström, 2010; Veitch and Nettles, 2012; Olsen and Nettles, 2017). The values of A are from mechanical modeling of capsizing iceberg contact forces (Sergeant and others, 2018) 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.</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>Value</th><th>Unit</th><th>Description</th><th>Filter</th></tr></thead><tbody><tr><td>Surface-wave magnitude</td><td>M<sub>SW</sub></td><td>4.6–5.1</td><td>–</td><td>Conventional earthquake surface-wave magnitude M<sub>S</sub> but computed on long-period bandpass filtered seismograms (Ekström, 2006a)</td><td>35–150 s</td></tr><tr><td>CSF magnitude</td><td>A<sub>CSF</sub></td><td>0.1 × 10<sup>14</sup><br>2 × 10<sup>14</sup></td><td>kg m</td><td>Result of the twice-time-integration of the CSF model which best reproduces the seismograms (e.g. Ekström and others, 2003; Veitch and Nettles, 2012)</td><td>40–100 s</td></tr><tr><td>Force history</td><td>F(t)</td><td>1.6 × 10<sup>10</sup><br>1 × 10<sup>12</sup></td><td>N</td><td>Time-evolution of the source force function inverted from seismograms (Sergeant and others, 2016) characterized by varying amplitudes over time and one maximum amplitude F<sub>max</sub></td><td>20–100 s</td></tr><tr><td>Contact force magnitude</td><td>A</td><td>0.5 × 10<sup>14</sup><br>1 × 10<sup>15</sup></td><td>kg m</td><td>Result of the twice-time-integration of the contact force generated by capsizing iceberg on the terminus</td><td>–</td></tr></tbody></table>
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, 2006a). 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 (2008) and Walter and others (2012) 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 (2016) 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, 2010, 2012b; Walter and others, 2013; Sergeant and others, 2016), 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 $\psi$ (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, 2017; Olsen and Nettles, 2017) and we force the azimuth angles $\psi$, when aligned with the glacier flow, to be consistent with up-glacier directions (i.e. in general, $0° \le \psi \le 180°$ for GEs located at Western Greenland glaciers, and $-180° \le \psi \le 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, 2015a). 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.

## Iceberg capsize model and contact forces

The iceberg capsize model of Sergeant and others (2018) 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

<figure><img src="image_3.png" /><figcaption>**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 (2015a) and James and others (2014). 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 $\epsilon H$) which tips over in bottom-out style in water with surface elevation $z_B + \Delta 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 superbuoyant as it experiences a positive $\Delta 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$.</figcaption></figure>

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 $\epsilon = W/H$ has to be lower than a critical value which depends on the ice-to-water density ratio (MacAyeal and others, 2003). 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
<figure><img src="image_4.png" /><figcaption>**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).</figcaption></figure>

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, 2014; Murray and others, 2015c), 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 $\Delta z$. Icebergs which experience a negative $\Delta z$ are then uplifted by buoyancy compared with their buoyant position at the initiation of the capsize. Conversely, for positive $\Delta 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 (1977). A catalog of contact force histories<fnref n="1" /> is then computed with varying sets of ($\epsilon$, $H$, $\Delta 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 \le \epsilon \le 0.6$, 600 m $\le H \le 1100$ m, $500 \text{ m} \le L \le 6000 \text{ m}$ and $-20 \text{ m} \le \Delta z \le 20 \text{ m}$.

## ICEBERG VOLUME ESTIMATION FROM FORCE HISTORY

Horizontal GE force components are compared with the catalog of modeled contact forces, with the latter filtered in

the same frequency band that was used for seismic inversion (20–100 s). The similarity of each pair of force time series is quantified via the variance reduction (VR) of their misfit. The calving volume inversion procedure consists of two steps. We first perform a ($\epsilon$, $H$, $\Delta z$) parameter space exploration with a contact force catalog corresponding to a given iceberg length $L_i$ and calculate an iceberg volume $V_0$ given by a best-fitting aspect ratio $\epsilon_0$ and height $H_0$ as $\epsilon_0 H_0^2 L_i$. The uncertainty on the retrieved volume, noted $\delta V_0$, is computed from the variation of the misfit function on $\epsilon$ and $H$. This step is then repeated for every $L$ and we compute the $V_0(L)$ function which is then used to calculate and average a final iceberg volume $V$ with statistical uncertainties, when consistent $V_0(L_i)$ values are obtained with $L$. All variables used in this study are reported in Appendix A.

## 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 $\epsilon$ and $H$ but one given length $L$ and at hydrostatic balance (i.e. $\Delta 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_0$ (Fig. 3a). When using noisy forces, maximum volume errors (15% for SNR=0.5 and 4% for

<footnote n="1">¹More precisely, since we carry out 2-D simulations, we recover linear force density which is computed in Newton per units of iceberg lateral length.</footnote>
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_0$ plotted in Fig. 3b is constructed by taking all $(ϵ, H)$ combinations that yield VR values >98% of maximum VR. While $δV_0$, as defined here, is always <15% for icebergs with $ϵ ≤ 0.2$, $δV_0$ 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, 2018, 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_0$ error or uncertainty $δV_0$ 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.

## 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 \sim 800$ m, an across-glacier length $L_d \sim 2500$ m, and an aspect ratio $ϵ_d \sim 0.23$ (Murray and others, 2015a). 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 \sim 800$ m, $ϵ_d \sim 0.37$ and $4000$ m $< L_d < 5000$ m (James and others, 2014). The GE horizontal force components shown in Fig. 4a (blue and green lines) have maximum amplitudes of $3 \times 10^{10}$ N and $4.6 \times 10^{10}$ N for the 2013 and 2010 events, respectively.

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 $(ϵ_0, H_0)$ and true iceberg dimensions $(ϵ_d, H_d)$, respectively. For the 2013 event, the best VR value (0.95) is achieved for $ϵ_0 = 0.2$, $H_0 = 850$ m and $Δz_0 = -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 (2015c). For the 2010 event, maximum VR (0.98) corresponds to a 11 m-uplifted iceberg with $ϵ_0 = 0.33$ and $H_0 = 825$ m. For both events, minimization-determined parameters are in good agreement with the true iceberg geometry and yield iceberg volumes $V_0$ of $\sim 0.36$ km³ and $\sim 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_0 < 7\%$) for the thinner iceberg, it increases up to 40% of $V_0$ for the larger calving event. Retrieved $δV_0$ 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, 2015a). 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, 2016).

The results of all inversion run with other $L$ values are reported in Fig. 4d. Best retrieved parameters $(ϵ_0, H_0)$ and volumes $V_0$ 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^3$ and $L$ (Tsai and others, 2008; Amundson and others, 2012a; Sergeant and others, 2018). 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 = ϵ_0H_0^2L_i$ that are estimated for a given iceberg length $L_i$, when accepted. This means that $V_0$, to be included in the final calculation of $V$, must (1) be stable around a certain value (i.e. $V_0(L)$ function must have reached a plateau), while (2) $ϵ_0$ is also stabilized and (3) $H_0$ does not exceed the actual glacier thickness $± 15\%$. In this study, we consider glacier thicknesses estimated in previous studies (e.g. Joughin and others, 2008a,b, 2014; Stearns and Hamilton, 2007; Howat and others, 2011; Murray and others, 2015c) 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_0$ 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, 2012) or recent ice thickness measurements from IceBridge surveys (Rignot and others, 2013). 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

$$δV = \frac{1}{N} \sqrt{\sum_{i}^{N} δV_0(L_i)^2} \quad (1)$$

where $δV_0(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
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 4.** Iceberg volume inversion for two well documented calving events of dimensions ($\epsilon_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 ($\epsilon_0$, $H_0$, $\Delta z_0$, $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 $\epsilon$ and $H$ for the catalog of contact forces scaled by the actual iceberg length ($L = L_d$). For this representation, the third model parameter $\Delta z$ has been fixed to the value reached for maximum VR, i.e. $\Delta z_0$ indicated above each panel. The real and calculated iceberg dimensions are ($\epsilon_d$, $H_d$) and ($\epsilon_0$, $H_0$), 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 ($\epsilon$, $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 ($\epsilon_0$, $H_0$) 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 $\epsilon_0$ is stable and $H_0$ does not exceed $\pm 15\%$ of the glacier thickness $H_d$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 5.** Size-frequency distributions of glacial earthquakes and inverted icebergs at individual Greenland glaciers for 1993–2013. Glacier locations are indicated in Fig. 8. Blue bars indicate CSF magnitudes $A_{CSF}$ from the Nettles catalog. Orange bars indicate the maximum amplitudes $F_{max}$ of the force time series that constitute our GE sub-catalog. Uncertainties on calving volume distributions come from the bin width (0.2 km³) used to generate histograms which is set to twice the median volume uncertainty obtained for all events.</figcaption></figure>

of $\delta V$ is 0.07 km³ across all inverted GEs. This represents an average volume uncertainty of 6%.

## RESULTS: ICEBERG PRODUCTION ASSOCIATED WITH GLACIAL EARTHQUAKES

Out of 15 GE producing glaciers in Greenland, seven are not included in this study as we fail to confidently determine GE sources and because not enough calving events could be characterized to make statistically relevant statements. We then focus our analysis on two Eastern major glaciers, Kangerdluqssuaq and Helheim, and six Western glaciers including Jakobshavn Isbrae, Rinks glacier and most active northwest glaciers, Upernavik Isström, Kong Oscar, Sverdrup, and Alison glaciers. Altogether, these glaciers cover 90% of the Greenland seismicity. Out of 203 GEs, the calving volumes of 28 events could not be estimated due to inconsistently calculated $V_0(L)$ and large VR

discrepancy. Finally, our sub-catalog of computed icebergs covers 40% of the total Greenland GE production and almost 45% of the production at considered glaciers.

As we aim to compute the evolution of capsizing calving volumes in Greenland and their contribution to the discharge, we need to assess the representativity of our sub-catalog of events. To do so, we look at the differences in magnitude patterns between our GE force database and the solutions for the complete event catalog. Fig. 5 shows the size-frequency distribution of GEs and computed calving volumes at each glacier. The proportion of inverted GEs with respect to the Nettles catalog varies from one glacier to another between 18 and 64% (Table 2, columns 1–3). Even if numerous events are missed (between 8 and 69 depending on the glacier), the distributions of the inverted seismic force amplitudes used in this study ($F_{max}$, orange bars) cover the entire range of CSF magnitudes and resemble $A_{CSF}$ histograms of the Nettles catalog (blue bars). Our GE
<figure><img src="image_12.png" /><figcaption>**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, 2007; Veitch and Nettles, 2012; Olsen and Nettles, 2017). $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 $\epsilon$. 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 $\epsilon$ are presented in (i), as well as the evolution of the magnitude of the contact force $A_c \propto H (1 - \epsilon) V$ (Eqn (3), color lines) when hydrodynamic effects are neglected and $H$ is fixed to 900 m.</figcaption></figure>

<table><caption>**Table 2.** Characteristics of glacial earthquake (GE) production and associated discharge (GED) computed at Greenland glaciers. Columns give number of events, proportion of inverted events, most common CSF magnitude $A_{CSF}$, most common inverted volume, GED in 1993–2013, and proportion of glacier dynamic discharge (DD) attributed to GEs for 2000–12. Columns 6–7 give results from GED median projections. Uncertainties are not reported here but are illustrated in Figs 7 and 12–14. GE contribution to the dynamic discharge DD averaged over East and West Greenland only accounts for DD measurements at the eight investigated glaciers.</caption><thead><tr><th>Glacier/<br>Region</th><th># of GEs<br>1993–2013</th><th>% of<br>inverted GEs</th><th>A<sub>CSF</sub> peak<br>(× 10<sup>14</sup> kg m)</th><th>V peak<br>(km<sup>3</sup>)</th><th>Associated discharge<br>GED (Gt)</th><th>GED contribution<br>to discharge DD (%)</th></tr></thead><tbody><tr><td>Helheim</td><td>104</td><td>38</td><td>0.4</td><td>0.7</td><td>79</td><td>12.4</td></tr><tr><td>Kangerdluqssuaq</td><td>84</td><td>18</td><td>0.7</td><td>1.5</td><td>107</td><td>15.5</td></tr><tr><td>Eastern Greenland:</td><td>188</td><td>27%</td><td>–</td><td>–</td><td>204 Gt</td><td>17%</td></tr><tr><td>Jakobshavn Isbrae</td><td>71</td><td>45</td><td>0.3</td><td>0.9</td><td>66</td><td>5.3</td></tr><tr><td>Rinks Isbrae</td><td>22</td><td>64</td><td>0.3</td><td>0.3</td><td>15</td><td>5.5</td></tr><tr><td>Upernavik Isström</td><td>39</td><td>34</td><td>0.3</td><td>0.3</td><td>16</td><td>16</td></tr><tr><td>Alison</td><td>14</td><td>65</td><td>0.4</td><td>0.7</td><td>36</td><td>13</td></tr><tr><td>Sverdrup</td><td>13</td><td>62</td><td>0.4</td><td>0.5</td><td>7</td><td>12.2</td></tr><tr><td>Kong Oscar</td><td>49</td><td>37</td><td>0.5</td><td>0.7</td><td>10</td><td>29.5</td></tr><tr><td>Western Greenland:</td><td>208</td><td>44%</td><td>–</td><td>–</td><td>163 Gt</td><td>10%</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>**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.</figcaption></figure>

database covers almost every seismically active year at each glacier and no specific temporal trend in GE magnitudes is observed at individual glaciers (Veitch and Nettles, 2012). This indicates that our sub-dataset would be representative of each glacier's overall seismicity in terms of temporal trend and GE sizes. It is worth noting that the two quantities that are compared, $F_{max}$ (with unit N) and $A_{CSF}$ (in kg m) have different natures and are obtained from different processing. We thus believe that discrepancies between the Nettles $A_{CSF}$ distributions and our $F_{max}$ is not necessarily related to incompleteness of our event sub-catalog.

## Earthquake sizes and calving volumes

As noticed by Veitch and Nettles (2012), 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, 2010; Veitch and Nettles, 2012). 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.

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 $\epsilon$ 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
<figure><img src="image_14.png" /><figcaption>**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.</figcaption></figure>

$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.

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, 2005, 2008; Joughin and others, 2008a; Bevan and others, 2012). 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
<table><caption>**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 (2014).</caption><thead><tr><th>GED contribution to DD</th><th>2000</th><th>2005</th><th>2012</th><th>2000–12</th></tr></thead><tbody><tr><td>Eight investigated glaciers</td><td>3–11% (6.5%)</td><td>12.5–21% (14%)</td><td>16.5–37% (20%)</td><td>8.5–21% (12%)</td></tr><tr><td>All Greenland tidewater glaciers</td><td>0.2–0.4% (0.3%)</td><td>1.4–2% (1.5%)</td><td>3.8–9% (4.7%)</td><td>2.5–6.5% (3.6%)</td></tr></tbody></table>

melting below the surface (Benn and others, 2017b). 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 (2014) 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.

## 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

$$C = \sum_{k}^{N} (V_{k} - \langle V \rangle) \quad (2)$$

where $V_k$ denotes iceberg volumes associated with N GEs within a given time period (e.g., 1 year) at one glacier, and $\langle V \rangle$ is the mean iceberg volume, averaged over all events produced at the glacier between 1993 and 2013. For constant discharge, $C \approx 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.

### Western glaciers

As observed by Veitch and Nettles (2012) and Olsen and Nettles (2017), GE production has increased since the early 2000s while Greenland glaciers accelerated during a period of above-average oceanic and atmospheric temperatures (Rignot and others, 2004; Holland and others, 2008; Murray and others, 2010; Straneo and others, 2010; Bevan and others, 2012). North West Greenland was also affected with an increased number of glaciers which began to retreat and accelerate. The fresh occurrence of GEs in North West Greenland is then likely related to changes in sea-ice and ice-mélange buttressing strength, glacier speeds, thinning rates and any successions of advancing and retreating phases (Moon and Joughin, 2008; Moon and others, 2015; Joughin and others, 2010; Howat and Eddy, 2011; McFadden and others, 2011; Carr and others, 2013) which may have brought the glacier calving fronts to near-grounded positions, facilitating seismogenic buoyancy-driven calving (Amundson and others, 2010). The growing numbers of GEs in West Greenland suggest that these glaciers still approach near-grounded and floating positions each year, mostly in summer (Olsen and Nettles, 2017). While before 2010, most discrete calving volumes were still small (i.e. $C<0$), it appears that iceberg sizes in this region tend to become larger also contributing to accelerated mass loss on the Western coast.

### 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, 2012). 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
<figure><img src="image_15.png" /><figcaption>**Fig. 9.** 1993–2013 front positions (orange) and speed (black) at (a) Helheim, (b) Kangerdluqssuaq and (c) Jakobshavn. Vertical bars indicate timing of glacial earthquakes (GEs). (d, e, f) Evolution of the cumulative GE-associated discharge anomaly C at corresponding glacier. 1993–2011 Eastern glacier positions are from Bevan and others (2012) and updated to 2013. 2009–13 Jakobshavn positions are from Joughin and others (2014).</figcaption></figure>

Eddy, 2011), 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, 2017). 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, 2014; Kehrl and others, 2017).

## 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, 2004). 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, 2014). In 2009–13, Jakobshavn remained near-grounded as the calving front was resting over a bed depression (Joughin and others, 2014). 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, 2016). 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, 2009; Joughin and others, 2012, 2014) and the development of a super-buoyant tongue whose length increases as ice advances (Benn and others, 2017a). 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, 2009) or glacier advance in deepening water (Wagner and others, 2016). 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.

GEs and calving patterns depend on the terminus front position which determines the flow speed and flotation level (Nick and others, 2009; Joughin and others, 2014). 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, 1998). Murray and others (2015c) 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
<figure><img src="image_16.png" /><figcaption>**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 (2014). 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 (2014).</figcaption></figure>

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, 2011; Carr and others, 2013; Moon and others, 2015) 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.

# DISCUSSION

## Calving volume inversion from glacial earthquakes

We have demonstrated that estimating robust calving volumes from GE force history is possible when averaged over consistent iceberg inversion results (i.e. $V_0(L)$). The fact that we fail to constrain iceberg shape (i.e. all three individual dimensions) does not result from our processing. It would rather be due to the difficulty to distinguish very similar contact forces that can be generated by capsizing icebergs of different geometry, and especially different height and length. Due to the symmetrical nature of the contact force amplitudes around $\epsilon \sim 0.4$ (Sergeant and others, 2018), the

comparison between seismic and modeled force histories for one given $L$ may give two peaks of the misfit function in the case of icebergs with a small aspect ratio or conversely, a large aspect ratio. For the 2013 Helheim event presented in Fig. 4c, we obtain two regions of a similar force fit, around $\epsilon = 0.2$ and $\epsilon = 0.47$ corresponding to VR values of 0.98 and 0.9, respectively. Volume error associated with the secondary VR peak would be 66% of the true value. Although the difference between the extreme VR peak values may be small, the best force model is always achieved in the right $\epsilon$ region as shown by synthetic inversions of calving volumes, even when noisy forces are considered (Fig. 3). The ambiguity in the iceberg volume is therefore possible to be resolved thanks to the aspect ratio's control on the force history.

In the following we discuss the errors we could obtain when trying to fit a calving volume to the earthquake magnitude $A_{CSF}$. Considering non-deforming icebergs and glacier termini, the contact force magnitude $A$ can be expressed as the sum of a geometric term $A_c$ which comes from the (twice-)integrated contact reaction, and a hydrodynamic term $A_w$. $A_w$ comes from the integration of drag forces and forces that arise from water pressure gradients that may be created by water deflection during iceberg tipping. The
<figure><img src="image_17.png" /><figcaption>**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 (2012) and updated to 2013. Glacier geometry in (c) is from Kehrl and others (2017).</figcaption></figure>

analytical form of $A_w$ is difficult to assess, though numerical modeling of drag forces (Amundson and others, 2012a; Sergeant and others, 2018) and hydrodynamic fields (Bonnet and others, 2018) show that accounting for water effects is critical for capturing capsize dynamics. However, following Tsai and others (2008) and Amundson and others (2012a, see equation 19), when neglecting hydrodynamics effects and considering thin enough icebergs ($\epsilon < 0.5$), the contact force magnitude $A$ would be equal to $A_c$ which scales as

$$A_c \approx \frac{1}{2} \rho_i H (1 - \epsilon) V, \quad \text{when } \epsilon < 0.5 \quad (3)$$

where $\rho_i$ is the density of the ice. The central question is, if we assume the iceberg height $H$ to be equal to the full glacier thickness $H_d$ we could deduce from additional glaciological measurements, is it possible to derive an appropriate scaling factor between GE magnitudes and calving volumes? If thin icebergs are considered (i.e. $1 - \epsilon \sim 1$), $A_c/V$ would be of the order of $10^6 \text{ kg.m}^{-2}$. However, any variation in $\epsilon$ alone may cause variations in the $A_c/V$ scaling factor of one order of magnitude leading then to errors in estimated volumes which can reach up to $1 \text{ km}^3$ when trying linear regressions between $A_{CSF}$ and $V$ (see also volume deviations with $\epsilon$ for same $A_c$ quantities in color lines, Fig. 6i). Then the extrapolation of GE magnitudes to calving volumes is difficult without high uncertainties. The comparison of GE force histories with adapted mechanical modeling of contact forces is then needed to extract accurate calving volumes.

One control on the time evolution of the contact force is the parameter couple ($\epsilon$, $\Delta z$). Quantifying the influence of $\Delta z$ relative to the one of $\epsilon$ or other iceberg dimension on the inversion results is not straightforward. Buoyant flexure of the glacier plays a significant role in calving (Wagner and others, 2016). Glacier elevation with respect to flotation level varies with the distance to the grounding line, and is sensitive to glacier thickness, glacier flow and water depth. The physical meaning of parametrization $\Delta z$ used here is not straightforward to link to glacier flotation level. Indeed the elevation of the nascent iceberg (captured in $\Delta z$) is not necessarily the same as for the glacier beam as it depends also on the shape of the iceberg submerged part. To better describe iceberg and glacier buoyancy states, more complex models which deal with viscous ice flow and subglacial topography including the complex crevasse network need to be used. Therefore, exact values of calculated $\Delta z$ are difficult to integrate into our analysis. Nevertheless, the most inverted $\Delta z$ describe uplifted icebergs (Fig. S3c). This result is reasonable and corresponds to observations of calving dynamics at Helheim (Joughin and others, 2008a; James and others, 2014; Murray and others, 2015c).

We were not able to compute the calving volumes of 14% of the 203 inverted GEs as no consistent volumes were reached when varying $L$. In this study we used a catalog of contact forces for bottom-out capsizes only. As hydrodynamics affect top-out and bottom-out capsizes differently, contact forces differ with calving style (Amundson and others, 2012a; Sergeant and others, 2018). Such deviations may justify our viability to invert coherent iceberg parameters for a portion
<figure><img src="image_18.png" /><figcaption>**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 (2014), (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).</figcaption></figure>

of the events. Indeed, top-out events are observed to be much less frequent in nature than bottom-out icebergs (Amundson and others, 2010). The discrimination between top-out and bottom-out events is difficult when only based on the horizontal seismic force. Further modeling efforts including real glacier terminus conditions should help to better distinguish between calving styles, when using also the capsize-generated vertical forces in the inversion procedure.

## Representativity of calving discharge

The estimation of Greenland discharge attributed to buoyancy-driven calving is reliable only if (1) our subset of seismic events is statistically representative for Greenland glacier seismicity as discussed in the previous sections, and

if (2) the Nettles catalog of GEs itself is also representative of capsizing calving in Greenland. In the following we discuss this further.

### Seismic detection of nontabular iceberg calving

Since they have been studied routinely, GE datasets have been complemented with different versions of the detection algorithm. Including regional Greenland stations in the detection process improved earthquake detections (Veitch and Nettles, 2012; Olsen and Nettles, 2017). However, many studies which use concurrent observations such as time-lapse photography, satellite imagery, local seismometers, or even water pressure sensors near the glacier report numerous capsizing calving events which remained undetected by standard processing (e.g. Amundson and
<figure><img src="image_19.png" /><figcaption>**Fig. 13.** Same as in Fig. 12 for Jakobshavn Isbrae, Rinks glacier and Upernavik Isström (West Greenland).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 14.** Same as in Fig. 12 for Alison, Sverdup and Kong Oscar glaciers (North West Greenland).</figcaption></figure>
others, 2010, 2012b; Walter and others, 2013). For example, Kehrl and others (2017) observed in the fjord of Helheim glacier ~70 nontabular icebergs between 2008 and 2013 while only 37 GEs were detected. Olsen and Nettles (2018) use cross-correlation templates at regional stations to identify smaller calving events (0.04–0.3 km³ measured from satellite imagery) which did not pass the detection threshold on teleseismic data. Olsen and Nettles (2018) estimate that undetected events could represent ice loss totaling an additional 10% of that accounted for by the glacial-earthquake standard catalog.

## Calving discharge

We estimated the volume of ice discharged in Greenland through GEs in 1993–2013 to be 398 km³ with lower and upper bounds being 270 km³ and 700 km³, respectively. Disparities between lowest, highest and median GED expectations come from the shape of the iceberg size-frequency distributions at individual glaciers and are difficult to reduce due to our analysis of a limited number of GEs. As outlined with Figs 9–11, the interannual variability of calving volumes is related to glacier geometry and buoyant state, flow velocity and ice-front position. When a glacier experiences large fluctuations in its fjord position from one year to another and speed acceleration, the annual average and predominant calving volume should shift. Estimates of GE-associated discharge could then be refined when averaged over different timescales.

Calving discharge was not calculated at seven glaciers which produced altogether 46 GEs (i.e. 10% of Greenland GE production). As predominant calving volumes vary from one glacier to another, it is difficult to assess the amount of GE-associated discharge that is missed in our study.

Finally, as discussed above, our calculated discharge (Tables 2 and 3) should give a lower bound to the discharge that is attributed to capsize events in Greenland. Besides, tidewater glaciers also produce numerous km-scale tabular icebergs and smaller icebergs which fall from ice cliffs or detach from underwater ice protrusions (Rignot and others, 2015; Wagner and others, 2016). It is then expected that calving is responsible for a much larger part of the GrIS dynamic mass loss.

## SUMMARY AND CONCLUSIONS

We developed a seismo-mechanical procedure for calculating calving volumes of seismogenic buoyancy-driven calving events which generate GEs in Greenland. The use of passive seismology is on its way to become a standard tool for investigating glacier processes and monitoring glacier changes (Podolskiy and Walter, 2016; Aster and Winberry, 2017). This study shows that for calving-related seismicity, classical seismic quantities such as the earthquake magnitude are difficult to interpret in terms of glaciological processes and mass loss as the seismic source is influenced by sea/ice interactions, iceberg and glacier geometries and hydrodynamic effects (Amundson and others, 2012a; Sergeant and others, 2018; Bartholomaus and others, 2012, for small non-capsizing icebergs). In the case of GEs and capsizing icebergs, the magnitude of the seismic source force cannot be used to uniquely define iceberg volumes, without errors of the order of 1 km³.

To refine the quantification of calving volumes, we propose to use the force history which captures the capsize dynamics that are primarily controlled by the iceberg shape. Whereas all three dimensions of icebergs are difficult to resolve well without any a priori postulation on at least one of them, iceberg volume can still be inverted. Indeed a robust calving volume can be estimated by averaging consistent inversion results run for several sets of iceberg dimensions. The confidence interval on the iceberg size varies with the iceberg geometry and can reach up to 40% for icebergs with intermediate width-to-height ratios. Synthetic tests give a maximum error on estimated volumes of 15% when run on poorly inverted force time series, which represents 0.3 km³ for the largest Greenland capsizing icebergs. We confidently inverted the iceberg sizes of ~200 capsizing events with an average uncertainty of 6% of the calving volumes.

We calculated the calving volumes of 40% of all GEs that were referenced in Greenland between 1993 and 2013. Thanks to preferred iceberg size-frequency distributions at individual glaciers, we could extrapolate the results to complete time series of GE-associated discharge at eight tidewater glaciers. Despite large uncertainties that are difficult to reduce due to our limited GE analysis, the total mass loss related to such events was estimated to be at least ~250 Gt (most probable value 367 Gt) over 1993–2013. While the iceberg discharge cumulated over the two decades has been so far equally partitioned between East and West Greenland, the contribution of Western glaciers has been accelerating since 2010 with an increasing number of calving events (Olsen and Nettles, 2017), but also increasingly larger produced icebergs. The comparison between ice-front position, GE occurrence and iceberg size indicate that glacier produced larger icebergs when they were advancing into deepening water. While the increasing numbers of GEs and related discharge that are observed since the early 2000s follows the general trend of widespread glacier acceleration, retreat, increased calving rates and discharge across Greenland, the variability in GE and calving production over time and space is to be related to the glacier geometry control on ice flow dynamics.

One third to one half of GrIS mass loss has been attributed to dynamic discharge at outlet tidewater glaciers, including calving and submarine melting. We estimate that buoyancy-driven calving events that produced GEs could be responsible for 8 to 21% of the total discharge recorded at the same glaciers (most probable value 12%) and between 3 and 6% of the overall ice-sheet loss. These estimations do not include other types of calving such as tabular icebergs and smaller icebergs. A large number of capsizing events which have produced earthquakes but remained undetected by the standard algorithms has recently been reported (e.g. Kehrl and others, 2017). This implies that our calculations on the contribution of such calving events to the Greenland mass loss could be far underestimated (at least 10% according to Olsen and Nettles, 2018). At Western glaciers, ice discharge through GEs has accelerated more rapidly than dynamic mass loss at the same glaciers. If trends continue and tidewater glaciers keep their near-grounded positions, buoyancy-driven calving and calving in general could become a major component of the GrIS mass loss.

We have illustrated in this study that GEs can help to quantify the partitioning of Greenland mass loss and monitor tidewater calving behavior which is crucial for a
better representation of calving processes (Benn and others, 2017a). Improvements in seismic detections and characterization of GEs are needed to expand the results of this study to a complete catalog of seismic and calving events. For a decade, increasingly high rate measurements of glacier characteristics and dynamics (e.g. James and others, 2014) have been collected. Reverse motion and acceleration of the glacier trunk were observed during and after capsize events (Murray and others, 2015a; Rosenau and others, 2013; Holland and others, 2016). The dynamics have been so far attributed to ice elastic deformation under capsizing iceberg forces (Murray and others, 2015a) and viscous flow adjustment to the loss of resistive stress at the terminus front (Howat and others, 2005; Joughin and others, 2008a). The processes controlling glacier adjustment to calving are still not resolved. They could involve frictional sliding or stick-slip events as observed in Antarctica (Bindschadler and others, 2003; Winberry and others, 2009). Coupling sophisticated models of glaciers with glaciological and seismic measurements offer a high potential for characterizing calving processes at tidewater glaciers and basal processes which are key unknowns for future sea-level rise estimation (Ritz and others, 2015).

## SUPPLEMENTARY MATERIAL

The supplementary material for this article can be found at https://doi.org/10.1017/aog.2019.7

## ACKNOWLEDGMENTS

The numerical modeling was done using Z-set (www.zset-software.com), a non-linear finite-element suite developed by the Centre des matériaux, Mines ParisTech and Onera (France). Seismic stations used in this study are part of the Greenland Ice-sheet Monitoring Network (Anderson and others, 2009) and seismic records can be accessed through IRIS data center (www.iris.edu). Seismic data were collected via ObsPy routines. Figures in this paper were created with Matlab. We acknowledge Douglas MacAyeal as the scientific editor and two anonymous reviewers who helped to improve the quality of the manuscript. We thank Meredith Nettles and Kira Olsen for providing the catalogs of seismic events (www.globalcmt.org/Events), Ellyn Enderlin who kindly shared ice discharge data, Clément Hibert and Guillaume Jouvet for fruitful discussions. This work and authors have been partially funded by the ANR contract ANR-11-BS01-0016 LANDQUAKES; CNCS-UEFISCDI project PN-II-ID-PCE-2011-3-0045, Institut Universitaire de France; the ERC contract ERC-CG-2013-PE10-617472 SLIDEQUAKES; the Swiss National Science Foundation for project GlaHMSeis PP00P2-157551; and a DGA-MRIS scholarship. This is IPGP contribution 4024.

## REFERENCES

* Amundson J and 5 others (2008) Glacier, fjord, and seismic response to recent large calving events, Jakobshavn Isbræ, Greenland. *Geophys. Res. Lett.*, **35**, L22501

* Amundson JM and 5 others (2010) Ice mélange dynamics and implications for terminus stability, Jakobshavn Isbræ, Greenland. *J. Geophys. Res.: Earth Surf. (2003–2012)*, **115**, F01005

* Amundson JM, Burton JC and Correa-Legisos S (2012a) Impact of hydrodynamics on seismic signals generated by iceberg collisions. *Ann. Glaciol.*, **53**(60), 106–112

* Amundson JM and 5 others (2012b) Observing calving-generated ocean waves with coastal broadband seismometers, Jakobshavn Isbræ, Greenland. *Ann. Glaciol.*, **53**(60), 79–84

* Anderson K and 10 others (2009) The Greenland ice sheet monitoring network (glisn). In *AGU Fall Meeting Abstracts*, volume 1, 0032

* Aster R and Winberry J (2017) Glacial seismology. *Rep. Prog. Phys.*, **80**(12), 126801

* Bartholomaus T, Larsen C, O'Neel S and West M (2012) Calving seismicity from iceberg-sea surface interactions. *J. Geophys. Res.: Earth Surf. (2003–2012)*, **117**(F4), F04029

* Benn DI and Åström JA (2018) Calving glaciers and ice shelves. *Adv. Phys.* **X**, **3**(1), 1513819

* Benn DI, Warren CR and Mottram RH (2007) Calving processes and the dynamics of calving glaciers. *Earth-Sci. Rev.*, **82**(3), 143–179

* Benn DI and 7 others (2017a) Melt-under-cutting and buoyancy-driven calving from tidewater glaciers: new insights from discrete element and continuum model simulations. *J. Glaciol.*, **63**(240), 691–702

* Benn DI, Cowton T, Todd J and Luckman A (2017b) Glacier calving in Greenland. *Curr. Clim. Change. Rep.*, **3**(4), 282–290

* Bevan S, Luckman A and Murray T (2012) Glacier dynamics over the last quarter of a century at helheim, kangerdlugssuaq and 14 other major Greenland outlet glaciers. *Cryosphere*, **6**(5), 923–937

* Bindschadler RA, King MA, Alley RB, Anandakrishnan S and Padman L (2003) Tidally controlled stick-slip discharge of a west Antarctic ice. *Science*, **301**(5636), 1087–1089

* Bonnet P and 7 others (2018) Mechanical modelling of iceberg capsize constrained by seismic inversion. In *Polar 2018, Open Science Conference*

* Carr JR, Vieli A and Stokes C (2013) Influence of sea ice decline, atmospheric warming, and glacier width on marine-terminating outlet glacier behavior in Northwest Greenland at seasonal to interannual timescales. *J. Geophys. Res.: Earth Surf.*, **118**(3), 1210–1226

* Ekström G (2006) Global detection and location of seismic sources by using surface waves. *Bull. Seismol. Soc. Am.*, **96**(4A), 1201–1212

* Ekström G, Nettles M and Abers GA (2003) Glacial earthquakes. *Science*, **302**(5645), 622–624

* Enderlin E, Howat I and Vieli A (2013) High sensitivity of tidewater outlet glacier dynamics to shape. *Cryosphere*, **7**(3), 1007–1015

* Enderlin EM and 5 others (2014) An improved mass budget for the Greenland ice sheet. *Geophys. Res. Lett.*, **41**(3), 866–872

* Enderlin EM, Hamilton GS, Straneo F and Sutherland DA (2016) Iceberg meltwater fluxes dominate the freshwater budget in Greenland's iceberg-congested glacial fjords. *Geophys. Res. Lett.*, **43**(21), 11–287

* Felikson D and 10 others (2017) Inland thinning on the Greenland ice sheet controlled by outlet glacier geometry. *Nat. Geosci.*, **10**(5), 366

* Fichefet T and 5 others (2003) Implications of changes in freshwater flux from the Greenland ice sheet for the climate of the 21st century. *Geophys. Res. Lett.*, **30**(17), 2003

* Griggs J and 10 others (2012) A new bed elevation dataset for Greenland. *Cryos. Discuss*, **6**, 4829–4860

* Hilber HM, Hughes TJ and Taylor RL (1977) Improved numerical dissipation for time integration algorithms in structural dynamics. *Earthq. Eng. Struct. Dyn.*, **5**(3), 283–292

* Holland DM, Thomas RH, De Young B, Ribergaard MH and Lyberth B (2008) Acceleration of Jakobshavn Isbrae triggered by warm sub-surface ocean waters. *Nat. Geosci.*, **1**(10), 659–664

* Holland DM and 10 others (2016) An intensive observation of calving at helheim glacier, East Greenland. *Oceanography*, **29**(4), 46–61

* Howat IM and Eddy A (2011) Multi-decadal retreat of Greenland's marine-terminating glaciers. *J. Glaciol.*, **57**(203), 389

* Howat IM, Joughin I, Tulaczyk S and Gogineni S (2005) Rapid retreat and acceleration of helheim glacier, East Greenland. *Geophys. Res. Lett.*, **32**(22), L22502
Howat IM, Joughin I and Scambos TA (2007) Rapid changes in ice discharge from Greenland outlet glaciers. *Science*, **315**(5818), 1559-1561

Howat IM, Joughin I, Fahnestock M, Smith BE and Scambos TA (2008) Synchronous retreat and acceleration of Southeast Greenland outlet glaciers 2000–06: ice dynamics and coupling to climate. *J. Glaciol.*, **54**(187), 646–660

Howat IM and 5 others (2011) Mass balance of Greenland's three largest outlet glaciers, 2000–2010. *Geophys. Res. Lett.*, **38**(12)

James TD, Murray T, Selmes N, Scharrier K and OLeary M (2014) Buoyant flexure and basal crevassing in dynamic mass loss at helheim glacier. *Nat. Geosci.*, **7**(8), 593–596

Joughin I, Abdalati W and Fahnestock M (2004) Large fluctuations in speed on Greenland's Jakobshavn Isbrae glacier. *Nature*, **432** (7017), 608–610

Joughin I and 8 others (2008a) Ice-front variation and tidewater behavior on helheim and kangerdlugssuaq glaciers, Greenland. *J. Geophys. Res.: Earth Surf.* (2003–2012), **113**(F1), F01004

Joughin I and 7 others (2008b) Continued evolution of Jakobshavn Isbrae following its rapid speedup. *J. Geophys. Res.: Earth Surf.*, **113**(F4), F04006

Joughin I, Smith BE, Howat IM, Scambos T and Moon T (2010) Greenland flow variability from ice-sheet-wide velocity mapping. *J. Glaciol.*, **56**(197), 415–430

Joughin I and 6 others (2012) Seasonal to decadal scale variations in the surface velocity of Jakobshavn Isbrae, Greenland: Observation and model-based analysis. *J. Geophys. Res.: Earth Surf.*, **117**(F2), F02030

Joughin I, Smith BE, Shean DE and Floricioiu D (2014) Brief communication: Further summer speedup of Jakobshavn Isbræ. *Cryosphere*, **8**(1), 209–214

Kehrl L, Joughin I, Shean D, Floricioiu D and Krieger L (2017) Seasonal and interannual variabilities in terminus position, glacier velocity, and surface elevation at helheim and kangerlus-suaq glaciers from 2008 to 2016. *J. Geophys. Res.: Earth Surf.*, **122**(9), 1635–1652

Khan SA and 10 others (2014) Sustained mass loss of the Northeast Greenland ice sheet triggered by regional warming. *Nat. Clim. Chang.*, **4**(4), 292

MacAyeal DR, Scambos TA, Hulbe CL and Fahnestock MA (2003) Catastrophic ice-shelf break-up by an ice-shelf-fragment-capsize mechanism. *J. Glaciol.*, **49**(164), 22–36

McFadden EM, Howat IM, Joughin I, Smith BE and Ahn Y (2011) Changes in the dynamics of marine terminating outlet glaciers in West Greenland (2000–2009). *J. Geophys. Res.: Earth Surf.*, **116**(F2), F02022

Meier M and Post A (1987) Fast tidewater glaciers. *J. Geophys. Res.: Solid Earth*, **92**(B9), 9051–9058

Moon T and Joughin I (2008) Changes in ice front position on Greenland's outlet glaciers from 1992 to 2007. *J. Geophys. Res.: Earth Surf.* (2003–2012), **113**(F2), F02022

Moon T, Joughin I, Smith B and Howat I (2012) 21st-century evolution of Greenland outlet glacier velocities. *Science*, **336**(6081), 576–578

Moon T and 6 others (2014) Distinct patterns of seasonal Greenland glacier velocity. *Geophys. Res. Lett.*, **41**(20), 7209–7216

Moon T, Joughin I and Smith B (2015) Seasonal to multiyear variability of glacier surface velocity, terminus position, and sea ice/ice mélange in Northwest Greenland. *J. Geophys. Res.: Earth Surf.*, **120**(5), 818–833

Murray T and 10 others (2010) Ocean regulation hypothesis for glacier dynamics in Southeast Greenland and implications for ice sheet mass changes. *J. Geophys. Res.: Earth Surf.* (2003–2012), **115**(F3), F03026

Murray T and 10 others (2015a) Reverse glacier motion during iceberg calving and the cause of glacial earthquakes. *Science*, **349**, 305–308

Murray T and 10 others (2015b) Extensive retreat of Greenland tidewater glaciers, 2000–2010. *Arct. Antarct. Alp. Res.*, **47**(3), 427–447

Murray T and 9 others (2015c) Dynamics of glacier calving at the ungrounded margin of helheim glacier, Southeast Greenland. *J. Geophys. Res.: Earth Surf.*, **120**(6), 964–982

Nettles M and Ekström G (2010) Glacial earthquakes in Greenland and Antarctica. *Annu. Rev. Earth. Planet. Sci.*, **38**, 467–491

Nettles M and 10 others (2008) Step-wise changes in glacier flow speed coincide with calving and glacial earthquakes at helheim glacier, Greenland. *Geophys. Res. Lett.*, **35**(24), L24503

Nick FM, Vieli A, Howat IM and Joughin I (2009) Large-scale changes in Greenland outlet glacier dynamics triggered at the terminus. *Nat. Geosci.*, **2**(2), 110–114

Olsen KG and Nettles M (2017) Patterns in glacial-earthquake activity around Greenland, 2011–13. *J. Glaciol.*, **63**(242), 1077–1089

Olsen K and Nettles M (2018) Analysis of regional seismic data reveals dominance of buoyancy-driven calving at Greenland glaciers. In *AGU Fall Meeting Abstracts*

Podolskiy EA and Walter F (2016) Cryoseismology. *Rev. Geophys.*, **54**(4), 708–758. (doi: 10.1002/2016RG000526)

Pratt MJ, Winberry JP, Wiens DA, Anandakrishnan S and Alley RB (2014) Seismic and geodetic evidence for grounding-line control of Whillans ice stream stick-slip events. *J. Geophys. Res.: Earth Surf.*, **119**(2), 333–348

Reeh N (1968) On the calving of ice from floating glaciers and ice shelves. *J. Glaciol.*, **7**, 215–232

Rignot E and Kanagaratnam P (2006) Changes in the velocity structure of the Greenland ice sheet. *Science*, **311**(5763), 986–990

Rignot E and 5 others (2004) Accelerated ice discharge from the Antarctic peninsula following the collapse of larsen b ice shelf. *Geophys. Res. Lett.*, **31**(18), L18401

Rignot E, Velicogna I, van den Broeke MR, Monaghan A and Lenaerts JT (2011) Acceleration of the contribution of the Greenland and Antarctic ice sheets to sea level rise. *Geophys. Res. Lett.*, **38**(5), L05503

Rignot E, Mouginot J, Larsen CF, Gim Y and Kirchner D (2013) *Icebridge wise I2 geolocated ice thickness and surface elevation*. National Snow and Ice Data Center, Boulder Colorado University

Rignot E, Fenty I, Xu Y, Cai C and Kemp C (2015) Undercutting of marine-terminating glaciers in West Greenland. *Geophys. Res. Lett.*, **42**(14), 5909–5917

Ritz C and 5 others (2015) Potential sea-level rise from antarctic ice-sheet instability constrained by observations. *Nature*, **528** (7580), 115

Rosenau R, Schwalbe E, Maas HG, Baessler M and Dietrich R (2013) Grounding line migration and high-resolution calving dynamics of Jakobshavn Isbræ, West Greenland. *J. Geophys. Res.: Earth Surf.*, **118**(2), 382–395

Sergeant A and 6 others (2016) Complex force history of a calving-generated glacial earthquake derived from broadband seismic inversion. *Geophys. Res. Lett.*, **43**(3), 1055–1065

Sergeant A and 5 others (2018) Numerical modeling of iceberg capsizing responsible for glacial earthquakes. *J. Geophys. Res.: Earth Surf.*, **123**, 3013–3033

Shepherd A and 10 others (2012) A reconciled estimate of ice-sheet mass balance. *Science*, **338**(6111), 1183–1189

Stearns LA and Hamilton GS (2007) Rapid volume loss from two East Greenland outlet glaciers quantified using repeat stereo satellite imagery. *Geophys. Res. Lett.*, **34**(5), L05503

Stern A, Adcroft A and Sergienko O (2016) The effects of Antarctic iceberg calving-size distribution in a global climate model. *J. Geophys. Res.: Oceans*, **121**(8), 5773–5788

Straneo F and 7 others (2010) Rapid circulation of warm subtropical waters in a major glacial fjord in East Greenland. *Nat. Geosci.*, **3** (3), 182–186

Sulak DJ, Sutherland DA, Enderlin EM, Stearns LA and Hamilton GS (2017) Iceberg properties and distributions in three Greenlandic fjords using satellite imagery. *Ann. Glaciol.*, **58**(74), 92–106

Tsai VC and Ekström G (2007) Analysis of glacial earthquakes. *J. Geophys. Res.: Earth Surf.* (2003–2012), **112**, F03S22
Tsai VC, Rice JR and Fahnestock M (2008) Possible mechanisms for glacial earthquakes. *J. Geophys. Res.: Earth Surf.* (2003–2012), **113**(F03014)

Van den Broeke M and 8 others (2009) Partitioning recent Greenland mass loss. *science*, **326**(5955), 984–986

Van den Broeke MR and 7 others (2016) On the recent contribution of the Greenland ice sheet to sea level change. *Cryosphere*, **10**(5), 1933–1946

Van der Veen C (1998) Fracture mechanics approach to penetration of bottom crevasses on glaciers. *Cold. Reg. Sci. Technol.*, **27**(3), 213–223

Veitch SA and Nettles M (2012) Spatial and temporal variations in Greenland glacial-earthquake activity, 1993–2010. *J. Geophys. Res.: Earth Surf.* (2003–2012), **117**(F002412), 1–20

Veitch SA and Nettles M (2017) Assessment of glacial-earthquake source parameters. *J. Glaciol.*, **63**(241), 867–876

Velicogna I, Sutterley T and Van Den Broeke M (2014) Regional acceleration in ice mass loss from Greenland and Antarctica using grace time-variable gravity data. *Geophys. Res. Lett.*, **41**(22), 8130–8137

Wagner TJ, James TD, Murray T and Vella D (2016) On the role of buoyant flexure in glacier calving. *Geophys. Res. Lett.*, **43**(1), 232–240A

Wagner TJ, Stern AA, Dell RW and Eisenman I (2017) On the representation of capsizing in iceberg models. *Ocean Modelling*, **117**, 88–96

Walter F and 5 others (2012) Analysis of low-frequency seismic signals generated during a multiple-iceberg calving event at Jakobshavn Isbræ, Greenland. *J. Geophys. Res.: Earth Surf.* (2003–2012), **117**(F1), F01036

Walter F, Olivieri M and Clinton JF (2013) Calving event detection by observation of seiche effects on the Greenland fjords. *J. Glaciol.*, **59**(213), 162–178

Wiens DA, Anandakrishnan S, Winberry JP and King MA (2008) Simultaneous teleseismic and geodetic observations of the stick-slip motion of an antarctic ice stream. *Nature*, **453**(7196), 770–774

Wilton DJ, Bigg GR and Hanna E (2015) Modelling twentieth century global ocean circulation and iceberg flux at 48 n: implications for West Greenland iceberg discharge. *Prog. Oceanogr.*, **138**, 194–210

Winberry JP, Anandakrishnan S and Alley RB (2009) Seismic observations of transient subglacial water-flow beneath macayeal ice stream, west Antarctica. *Geophys. Res. Lett.*, **36**(11), L11502

Yang Q and 7 others (2016) Recent increases in Arctic freshwater flux affects labrador sea convection and atlantic overturning circulation. *Nat. Commun.*, **7**, 10525

# APPENDIX

## NOTATIONS

<table><thead><tr><th>Symbol</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>A</td><td>Magnitude of the contact force. Calculated from iceberg capsize modeling</td><td>[kg m]</td></tr><tr><td>A<sub>CSF</sub></td><td>Magnitude of the centroid single force. Calculated from seismic waveform inversion</td><td>[kg m]</td></tr><tr><td>F</td><td>Force history produced by the iceberg at the terminus</td><td>[N]</td></tr><tr><td>H</td><td>Iceberg height</td><td>[m]</td></tr><tr><td>H<sub>0</sub></td><td>Calculated iceberg height given one iceberg length L<sub>i</sub></td><td>[m]</td></tr><tr><td>H<sub>d</sub></td><td>True iceberg height from field measurements or input data</td><td>[m]</td></tr><tr><td>L</td><td>Iceberg across-glacier length</td><td>[m]</td></tr><tr><td>L<sub>d</sub></td><td>True iceberg length from field measurements or input data</td><td>[m]</td></tr><tr><td>M<sub>SW</sub></td><td>Magnitude of glacial earthquake surface waves</td><td>—</td></tr><tr><td>V</td><td>Final calculated iceberg volume computed as V = 1/N &sum;<sub>i=1</sub><sup>N</sup> V<sub>0</sub>(L<sub>i</sub>) when volumes V<sub>0</sub>(L<sub>i</sub>) are consistent between each other</td><td>[m<sup>3</sup>]</td></tr><tr><td>V<sub>0</sub></td><td>Calculated iceberg volume given one iceberg length L<sub>i</sub>, V<sub>0</sub> = &epsilon;<sub>0</sub>H<sub>0</sub><sup>2</sup>L<sub>i</sub></td><td>[m<sup>3</sup>]</td></tr><tr><td>V<sub>d</sub></td><td>True iceberg volume from field measurements or input data, V<sub>d</sub> = &epsilon;<sub>d</sub>H<sub>d</sub><sup>2</sup>L<sub>d</sub></td><td>[m<sup>3</sup>]</td></tr><tr><td>W</td><td>Iceberg along-glacier width, W = &epsilon;H</td><td>[m]</td></tr><tr><td>z<sub>B</sub></td><td>Water level for the iceberg to be at flotation at the initiation of its capsize</td><td>[m]</td></tr><tr><td>z<sub>w</sub></td><td>Water level which defines the hydrostatic balance of the iceberg at the initiation of its capsize, z<sub>w</sub> = z<sub>B</sub> + &Delta;z</td><td>[m]</td></tr><tr><td>&epsilon;</td><td>Iceberg aspect ratio, &epsilon; = W/H</td><td>—</td></tr><tr><td>&epsilon;<sub>0</sub></td><td>Calculated iceberg aspect ratio given one iceberg length L<sub>i</sub></td><td>—</td></tr><tr><td>&epsilon;<sub>d</sub></td><td>True iceberg aspect ratio from field measurements or input data</td><td>—</td></tr><tr><td>&Delta;V</td><td>Uncertainty on the final calculated iceberg volume V and computed as &Delta;V = 1/N &sum;<sub>i=1</sub><sup>N</sup> &Delta;V<sub>0</sub>(L<sub>i</sub>)<sup>2</sup> with &Delta;V<sub>0</sub>(L<sub>i</sub>) corresponding to the volume uncertainty associated with V<sub>0</sub>(L<sub>i</sub>)</td><td>[m<sup>3</sup>]</td></tr><tr><td>&Delta;V<sub>0</sub></td><td>Uncertainty on the iceberg volume V<sub>0</sub> calculated for one iceberg length L = L<sub>i</sub></td><td>[m<sup>3</sup>]</td></tr><tr><td>&Delta;z</td><td>Perturbation of water level z<sub>B</sub> that is needed for the iceberg to be out of flotation at the initiation of its capsize, &Delta;z = z<sub>w</sub> - z<sub>B</sub></td><td>[m]</td></tr><tr><td>&Delta;z<sub>0</sub></td><td>Calculated &Delta;z given one iceberg length L<sub>i</sub>. Captures the buoyancy state or the iceberg prior to its capsize, i.e. &Delta;z<sub>0</sub> ~ 0 when the iceberg is at flotation, &Delta;z<sub>0</sub> > 0 when the iceberg is super-buoyant (below its flotation level), &Delta;z<sub>0</sub> > 0 when the iceberg is aerial (above its flotation level)</td><td>[m]</td></tr><tr><td>&psi;</td><td>Azimuth of glacial earthquake inverted force</td><td>[degree]</td></tr></tbody></table>