A glacier–ocean interaction model for tsunami genesis due to iceberg calving

Glaciers calving icebergs into the ocean significantly contribute to sea-level rise and can trigger tsunamis, posing severe hazards for coastal regions. Computational modeling of such multiphase processes is a great challenge involving complex solid–fluid interactions. Here, a new continuum damage Material Point Method has been developed to model dynamic glacier fracture under the combined effects of gravity and buoyancy, as well as the subsequent propagation of tsunami-like waves induced by released icebergs. We reproduce the main features of tsunamis obtained in laboratory experiments as well as calving characteristics, the iceberg size, tsunami amplitude and wave speed measured at Eqip Sermia, an ocean-terminating outlet glacier of the Greenland ice sheet. Our hybrid approach constitutes important progress towards the modeling of solid–fluid interactions, and has the potential to contribute to refining empirical calving laws used in large-scale earth-system models as well as to improve hazard assessments and mitigation measures in coastal regions, which is essential in the context of climate change. Dynamic glacier fracture and the subsequent generation and propagation of iceberg-induced tsunamis are reproduced using a unified numerical glacier–ocean model, in line with observations at the Eqip Sermia glacier in Greenland, as well as laboratory experiments.

G lacier calving into the ocean (Fig. 1) is predicted to be one of the largest contributions to sea-level rise in the future [1][2][3] . This process corresponds to~50% of the mass loss from ice sheets in Greenland and Antarctica 4,5 . Glacier calving is caused by first, second, and third-order processes. Firstorder processes correspond to the formation of surface crevasses in the ice owing to spatial variation in flow velocity. Second and third-order processes include crack propagation owing to local stress concentrations, ice stretching in the vicinity of the ice front, and oceanic erosion and torque induced by buoyant forces 6,7 . Depending on the shape of the glacier outlet, this may lead to different calving scenarios 8 . Glacier calving can have dramatic consequences, as falling or capsizing icebergs can generate large tsunamis, which are a threat to coastal infrastructure, ecology, and people [9][10][11][12][13][14][15] . In high mountain proglacial lakes, calvinginduced waves may further pose major hazards through triggering lake outburst floods with high destructive potential in downstream valleys 16 .
Most existing numerical approaches for marine-terminating glaciers were developed to study the slow creep of ice using continuum Eulerian methods (e.g., Elmer/Ice 17 ), and the calving rate is generally evaluated using simplified and empirically based calving laws or simple analytic models 18 . In general, the validation of these models against observations remains relatively limited 19 and mostly excludes temporal scales of single calving events. Åström et al. 5,20 and Bassis and Jacobs 21 developed purely Lagrangian particle-based models, based on a discrete version of Newton's equations of motion, to study the dynamics of sea ice and glacier calving. Despite several approximations, including a simplified water-ice interaction law, their simulations were able to reproduce the fractal nature of the debris size distribution 20 and diverse calving features based on glacier geometry 21 . However, the discrete nature of these models makes them computationally very expensive and therefore limited to single events. Furthermore, the water in their simulation is not explicitly modeled, which prevents the study of the tsunamigenic potential of glacier calving. More recently, Mercenier et al. 22,23 developed a transient multiphysics finite element model to simulate the effect of oceanic melt on ice break-off at the terminus of a marine glacier. They showed that a von Mises stress criterion led to realistic calving front geometries. Although this model successfully couples slow glacier flow with a damage-based calving criterion, it cannot simulate tsunamis induced by calving.
Tsunamis generated by landslides were extensively studied 24,25 . Yet, only a few studies focused on tsunamis induced by calving glaciers. Lüthi and Vieli 10 reported a 45-50 m tsunami generated by the calving of a 200 m high ice cliff of the calving front of the Eqip Sermia glacier in Greenland. The wave was still 3.3 m large 4.5 km from the calving outlet and led to a 20 m run-up on the opposite shore. Recently, Heller et al. 13 performed large-scale laboratory experiments to study the characteristics of waves generated by different calving mechanisms. These authors showed that empirical equations established for landslidesinduced tsunamis 25 overestimated wave amplitudes and generally fail to reproduce the physics of calving-induced tsunamis. Recently, Chen et al. 26 were able to reproduce the characteristics of waves reported in Heller et al. 13 using foam-extend and the Immersed Boundary Method. Yet, an approach to simulate fracture processes during dynamic glacier calving and tsunamis in a unified manner is still missing.
Here, we report coupled glacier calving and tsunami experiments and develop a continuum damage Material Point Method (MPM) to explicitly simulate ice fracture and hydraulic interactions. This new model accurately reproduces dynamic ice fractures and generated tsunami characteristics for different calving mechanisms for laboratory and real-world scales.

Results
Ice and water mechanics. To model the dynamic fracture of the glacier ice, we developed a non-associative elastoplastic model based on the Cohesive Cam Clay (CCC) yield surface used by Gaume et al. 27 to simulate snow and avalanche mechanics. A mixed-mode yield surface such as this was shown to adequately model brittle ice fracture based on experimental data 28 . However, the previously chosen associative flow rule was only adequate owing to the porous nature of snow, allowing for volume change (compaction hardening). Conversely, in the case of a significantly less-porous material such as ice, choosing a non-associative flow rule 29 is key owing to its natural volume-preserving qualities 30 . As such, we adopt a non-associative flow rule 31  softening law to model the dynamic ice fracture. The mechanical behavior of water is simulated using a nearly incompressible equation of state 32 . The mass and momentum balance equations are solved numerically using the MPM 27,33 (see Methods for a complete model description).
Calving mechanism and tsunami wave characteristics. In a first set of experiments, we model the tsunami wave response to three main calving mechanisms according to the laboratory experiments of Heller et al. 13 : (i) gravity-dominated fall (GF), (ii) buoyancy-dominated fall (BF), and (iii) capsizing (CS). The numerical setup consists of the interaction between a rigid block of density similar to that of ice (920 kg/m 3 ) and a D w = 1 m deep water tank (see Methods) in two dimensions. For all simulated cases ( Fig. 2 and Supplementary Movie 1), a rapid increase of the wave amplitude is reported until it reaches a maximum value. Then, the wave amplitude significantly decreases. A 2D/3D amplitude transformation factor based on an extensive amount of experimental data 34 (see Methods) was used to compare our results to those of Heller et al. 13 . A good agreement between our numerical solution and experimental data is observed. In particular, the maximum wave amplitude is correctly reproduced in both capsizing and gravity fall cases and slightly overestimated in the buoyant fall case (Fig. 2b). Concerning the wave velocity, our simulation results perfectly follow the theoretical limit prediction for linear shallow-water waves v w ¼ ffiffiffiffiffiffiffiffi gD w p ; furthermore, a fair agreement with experimental data is reported. Consistent with laboratory experiments, our model results also suggest that gravity fall cases are the most hazardous as they induce the largest tsunami amplitudes.
Iceberg sizes: interplay between buoyancy and ice weight. To evaluate the capabilities of our new constitutive model for ice, we perform glacier calving MPM simulations using a simple geometry for which analytical solutions exist for validation: an ice block of height H i is sliding at a constant velocity and enters a water tank at varying submergence depth D to modify the buoyancy contribution and thus, ice fracture features. An example of simulation result for D/H i = 0.2 is shown in Fig. 3a and Supplementary Movie 2. In this case, the main contribution to instability is gravity. Hence, we observe ice fractures starting to form at a certain length behind the front, and that propagate from the top of the glacier owing to gravity-induced bending. Once released, the iceberg generates a tsunami that propagates to the other side of the water tank (see Supplementary Movie 2). If the simulation is run for a long time, we observe subsequent detachments of icebergs of similar sizes, hence this length measure seems a robust output quantity of the model. The first iceberg length L i is computed and plotted in Fig. 3b as a function of the submergence depth-ice height ratio D/H i . The iceberg length first increases as submergence depth increases owing to a progressive stabilizing effect of buoyancy. For D/H i~0 .92 = ρ i /ρ w , no ice failure is reported as the bending force due to the ice weight is compensated by the effect of buoyancy. For D/H i > 0.92, we observe ice fractures starting from the bottom of the ice slab related to the effect of buoyant forces which overcome body weight forces. In turn, the iceberg length decreases with increasing D/H i . Additional simulation results are shown in the Supplement for D/H i = 0.7 and for lower ice strength values. Numerical results were then compared to a simple analytical model which compares the bending stress induced by the interplay between buoyancy and gravity and the ice tensile strength. This model (see Methods) leads to the following analytical expression for the iceberg length: where σ t is the tensile strength of ice. Our MPM simulation results are in line with the predictions of this theoretical bending model (Fig. 3b) and are consistent with the position of the modeled stress maxima by Mercenier et al. 22 . The theoretical prediction for linear shallow-water waves is given by Calving-induced tsunami at the Eqip Sermia glacier. Given the ability of our model to quantitatively reproduce wave and ice fracture characteristics, we simulate one past calving event of the Eqip Sermia, a tidewater outlet glacier in Greenland, for which numerous observations and field measurements are available for both the glacier and tsunamis 10 . The geometry of the glacier front has been initialized according to the real shape of the glacier 10 . The undulating "crevassed" surface structure has been produced according to a Simplex Improved Perlin noise 35 to mimic the crevasse distribution of the glacier ( Fig. 4 and Supplementary Movies 3 and 4). The detailed description of this event, the numerical setup, and the parameter calibration are given in the Methods. The failed portion of the glacier has a width of 75 m in our simulations which is slightly larger than the reported value measured between 50 and 60 m. The angle of the failure plane in our simulation is 54°, in good agreement with the measurements (between 45 and 60°). Note that the size of the initial iceberg is not only dependent on the relative water level, but also on the crevasse distribution, as shown in Fig. 4 where the failure is initiated at the location of the first deeper crevasse. The wave amplitude measured in our MPM simulation ( Fig. 4) 250 m from the calving front is 50 m, which is in good agreement with the wave amplitude observed at this location (between 45 and 50 m) 10 . A wave amplitude of 3.3 m was measured at the tide gauge on the shore opposite of the glacier. In our 2D simulation, the wave amplitude decreases from~40 to 30 m at a distance of 2000 m from the calving front and then slowly decreases further. We measure a wave amplitude of 27 m at the distance of the tide gauge, which corresponds to a wave amplitude between 2.5 and 6.6 m in 3D (see Methods) also in good agreement with the reported field measurement. Finally, we measure an average wave speed of 35 m/s, which also agrees well with the measured value of 32 m/s; furthermore, this agrees with theoretical predictions for a fjord water depth estimated between 100 and 150 m 36 , leading to a theoretical shallow-water velocity range between 31 and 38 m/s.

Discussion
We proposed a new model based on the MPM, a hybrid Eulerian-Lagrangian numerical technique, to simulate glacier calving and tsunami generation. Our nearly incompressible water model was able to reproduce well the effects of different iceberg calving mechanisms on the amplitude and speed of the generated tsunamis measured in laboratory experiments 13 . Deviations in maximum wave amplitude can be attributed to the empirical nature of the 2D/ 3D transformation (±50% deviation reported in Heller and Spinneken 34 ). In addition, slight wave velocity deviations from the linear shallow-water limit theory can be attributed to the intermediate character and non-linearity of the waves as well as frequency dispersion for which variations around the linear wave theory expression can reach around ±10% for a w /D w > 0.03 37 . This wave velocity difference may contribute to the faster amplitude decrease in the gravity-dominated fall and capsizing cases. Our model is intended to simulate the complex multiphase interplay between ice and water mechanics during glacier calving. Although our water model has been extensively used in Smoothed Particle Hydrodynamics (SPH) free surface flow simulations of water and has been validated for several applications 32,38-40 , solving incompressible Navier-Stokes equations using computational fluid dynamics would lead to more accurate results for the water dynamics 41 . Yet, it would also computationally be significantly more expensive and challenging to implement within the current hybrid solid-fluid numerical framework. Our new ice mechanical model has been validated by simulating a simple glacier calving geometry for which the competing effects of the glacier weight and buoyant forces lead to different ice fracture mechanisms and iceberg lengths that agreed well with theoretical predictions. When applying the model to a real-world outlet glacier, the model quantitatively reproduced the iceberg geometry and dimensions as well as the main characteristics of the generated tsunamis. The simulated iceberg length largely depends on the crevasse distribution and the chosen value of the ice's tensile strength. We back-calculated the latter such that the model results are in good agreement with the observations. Yet, we verified that the obtained value was within a realistic range based on extensive laboratory research on ice mechanics 28,42 . Nevertheless, glacial ice has high variability in its mechanical properties as well as deep crevasses leading to strong anisotropy, which both have not been accounted for here. In turn, one can expect smaller iceberg size distributions in reality 43 (see Supplementary Figs. 2 and 3). It is important to note that our model only applies to short time scales for which the ice behaves as a brittle material. We do not simulate viscous creep, nor the complex basal processes leading to crevasse formation. We could implement an ice creep law in our model, but the time-steps required to capture fast ice fracture are incompatible with such slow, viscous deformations. However, one could easily couple the present model with a standard continuum glacier flow model 17,22 . In addition, note that an assumption in the formulation of our ice plastic flow rule gives rise to an analytical solution that speeds up simulation (see Supplementary Material and Simo 44 ). This simplified formulation is based on the hypothesis of a relatively stiff material (see Eq. (5) in the Supplement), which is fully appropriate for ice. An iterative return mapping algorithm would be required for very soft matter (see Gaume et al. 27 ).
The major advantage of our approach compared with existing methods 5,20,22,23 is the ability to simulate, in a unified manner, all of the processes related to fast glacier calving including dynamic ice fracture, tsunami formation, and propagation. In addition, we can apply our method to any type of geometry, including complex surface crevasses and boundary conditions. This allows us to analyze one or multiple subsequent calving events that can be measured, which is crucial for model validation. This unified and multiphysics approach prevents us from using model chains 5,8 , which suffer from error propagation. Furthermore, the hybrid Eulerian-Lagrangian continuum framework leads to a tremendous reduction in the computational time compared with the . These advantages and the good agreements with laboratory and field data make our approach extremely well suited for glacier tsunami hazard assessment in coastal regions. Finally, our model can readily be applied to study lake outburst of glacial or mountain lakes induced by glacier calving 16 , landslides, snow, ice, or rock avalanches [45][46][47] .

Methods
A non-associative elastoplastic model for dynamic ice fracture. To model the dynamic hardening and softening of the glacial ice, we adapt and expand a nonassociative elastoplastic model based on the aforementioned CCC yield surface used by Gaume et al. 27 . This yield surface showed great success for modeling the plasticity of snow by incorporating its natural cohesive properties; however, the previously chosen associative flow rule was only adequate owing to the porous nature of snow, allowing for some volume change. Conversely, in the case of a significantly less-porous material such as ice, choosing a non-associative flow rule is key owing to its natural volume-preserving qualities 30 . As such, we adopt the non-associative flow rule proposed by Wolper et al. 31 to model the dynamic fracture of a variety of solids while ensuring volume preservation. Unfortunately, this approach brings with it a unique problem in the discrete treatment of plasticity: when stresses are projected to the surface orthogonally to the hydrostatic axis, there is no change in pressure and, as such, no analytic way to compute the hardening/softening. Wolper et al. found success in treating this issue using a geometric intersection approach; however, we propose a more physically based approach that uses quantities we have analytic expressions for.
First, we outline the elastoplastic theory behind our return mapping approach. We use the deformation gradient, F = ∂ϕ/∂X, where ϕ(X, t) denotes a deformation map from the undeformed coordinate X to the current configuration x. Furthermore, we decompose F into an elastic part, F E , and a plastic part, F P , as F = F E F P .
We adopt a split Neo-Hookean hyperelasticity model to predict nonlinear stress-strain relations suitable to simulating materials undergoing large deformations 48 ; this deviatoric-dilational split energy allows for the separate penalization of shearing and volume change and is written as: ð tr ðF T FÞ À dÞ and Ψ κ ðJÞ ¼ κ 2 with d = 2 or 3 being the problem dimension, J ¼ detðFÞ, and μ and κ being the shear and bulk modulus, respectively. Using this model, we can also write out the Kirchoff stress, τ, and its deviatoric part, s, as: Our yield surface is defined based on critical state soil mechanics 49 and depends on the first (pressure p) and second (Mises equivalent stress q) invariants of the Kirchoff stress tensor τ. We use the convention that compression corresponds to p > 0. Recent triaxial loading laboratory experiments showed that the yield surface of ice has an elliptical shape in the p−q space 50 . As a consequence, we chose the following CCC yield surface similar to that used for snow by Gaume et al. 27 and Meschke et al. 51 : where β is the cohesion coefficient, M is the slope of the cohesionless critical state line and p 0 is the pre-consolidation pressure. Softening or hardening is performed by shrinking or expanding the yield surface through variations in p 0 (Fig. 6). The major difference compared with Gaume et al. 27 is in the plastic flow rule and the hardening variable used for the hardening/softening relationship. Snow is a highly porous material, which justified the use of an associative plastic flow rule and a hardening law based on the volumetric plastic strain. In contrast, the failure of ice is not associated with large changes in volume 29 . Hence, a nonassociative plastic flow rule has been developed, and a relevant formulation for the deviatoric plastic strain, α (see definition in Supplementary Note 1 and illustration in Supplementary Figure 1), is used in the hardening rule as follows: where κ ¼ 2 3 μ þ λ is the bulk modulus and ξ is the hardening factor. A q-based non-associative projection to the yield surface is performed if the trial elastic pressure p tr is between −βp 0 and p 0 . If this is not the case, projection is made to the tips of the yield surface (Fig. 6). The rate of deviatoric plastic deformation is positive ( _ α > 0) if p < p c , leading to material softening; conversely, the rate is negative ( _ α < 0) if p > p c , inducing material hardening (see details in Supplementary Note 1). The proposed deviatoric hardening variable is also different from that proposed by Wolper et al. 31 who used a non-physically based variable formed from simple geometric considerations.
A nearly incompressible fluid model for water. Water is modeled as a nearly incompressible fluid 32 with stress tensor σ w and partial stress p w that depends on the current water density, ρ, and initial water density, ρ w , according to: where κ w is the water bulk modulus and θ is a parameter that penalizes deviations from incompressibility. This relationship, which is classically used in water simulations using SPH, another particle-based method 32 , is designed to stiffly prevent volume change of the water phase.
Material Point Method. MPM discretizes a continuum into Lagrangian material points to track mass, momentum, and deformation gradient, and uses a background Eulerian grid to solve mass and momentum conservation as follows: where ρ is density, t is time, v is velocity vector, σ is the Cauchy stress tensor, and g is the gravitational acceleration vector. In MPM, mass is automatically conserved as the mass of each particle is constant. The momentum balance is solved on the background grid through discretization of its weak form. The explicit MPM algorithm by Stomakhin et al. 52 is applied with a symplectic Euler time integrator. Details of the adopted MPM time-stepping algorithm can be found in 27,52,53 . In contrast with Gaume et al. 27 , we use the Affine Particle-In-Cell (APIC) method 54 to perform MPM transfer operations, which allows for better conservation of both linear and angular momentum. The Cauchy stress tensor σ in Eq. (8) is related to the strain as follows where Ψ is the elastoplastic potential energy density, F E is the elastic part of the deformation gradient, F, and J ¼ detðFÞ.
Iceberg length theoretical model. Our choice to simulate glacier calving on a simple geometry is motivated by the existence of a simplified analytical solution for model validation. This solution has been derived based on the bending stresses induced by gravitational forces (from the interplay between the ice weight and buoyancy) using beam theory 55 . Ice fracture requires that the total bending stress induced by this interplay exceeds the ice's tensile strength, σ t . Bending stresses can be computed according to σ = Mz/I where M is the bending moment, z is the vertical coordinate (above the neutral surface), and I is the moment of inertia. This leads to the following instability criterion: which leads to Eq. (1).
Experimental data 2D/3D transformation. The general expression for the transformation of the 2D to 3D primary wave amplitudes is given by Heller and Spinneken 34 : where r is the radial distance from the impact, D w is the water depth, γ is the wave propagation angle, and the parameters F, S, and M are the slide Froude number, the relative slide thickness, and the relative slide mass, respectively, and defined as follows: where v s is the slide velocity, s is the slide thickness, m s is the slide mass, and b is the slide width. The function f γ is defined as Iceberg-tsunami laboratory experiments. We conducted the experiments in the Delta Basin of Deltares, Delft, with effective size, excluding wavemakers and absorbing beaches, of 40.3 m × 33.9 m. The total number of experiments was 66.
The following provides an overview of the three herein selected experiments with further details for all experiments given by Heller et al. 13,37 . Table 1 shows the relevant parameters. We modeled the icebergs with polypropylene homopolymer blocks with a density similar to that of ice. The 0.500 m × 0.800 m × 0.500 m block in the gravity-and buoyancy-dominated fall cases was released at the vertical boundary of the basin and the 0.500 m × 0.800 m × 0.250 m block in the capsizing case was released offshore. The blocks weighed up to 187.1 kg and the still water depth was 1.00 m.
Three types of experiments were conducted: (i) GF, (ii) BF, and (iii) CS. For GF cases, we held the block in position with an electromagnet prior to release, which was connected to a rope. We fixed the supporting frame for this electromagnet and the block to a steel plate at the basin wall. The block was moved in the vertical direction towards the release position using a winch system fixed to a support structure outside the basin. For BF cases, the block was pulled underwater to the release position with a rope attached to the center of the block base. We further stabilized the block with a steel beam from above ( Fig. 7) and released both the beam and the rope simultaneously to initiate waves. For CS cases, we held the block in position with a wooden rod guided through the center of the block. Steel profiles on both sides accommodated this rod such that the block was able to heave and pitch, but not to sway and surge. We initiated capsizing by removing a fitting that stabilized the block and then simply releasing the block by hand (Fig. 7).
The velocity v s corresponds to the fastest moving section of the block derived from the recordings of a nine Degree of Freedom motion sensor (Adafruit BNO055) with a procedure described in Chen et al. 26 . The sensor was located in a watertight enclosure and attached to the block surface. We extracted the still images shown in Fig. 7 from camera footage recorded with a 5 MP PointGrey ZBR2-PGEHD-50S5C-CS camera at 15 Hz. We used up to 35 resistance-type wave probes to record the wave features. Relative to the coordinate origins at the front of the steel plate in the center of the block in the cross-shore direction (GF and BF cases) and at the block center (CS case), we arranged these probes at different points in the range of relative radial distances r/D w = 2-35 and wave propagation angles γ = 0 ∘ (block axis) to −75 ∘ (fall cases) and r/D w = 2-15 and γ = 0 ∘ to −180 ∘ (capsizing case). The angle γ is thereby defined positively in the clockwise direction. We shortened the water surface time series individually to remove data affected by reflection from the basin boundaries and filtered most wave probe data with a lowpass filter with a cutoff frequency at 9-11 Hz 37 .
Eqip Sermia glacier calving-tsunami field experiments. Glacier calving activity at Eqip Sermia was investigated with a terrestrial radar interferometer (TRI), surveys with unmanned aerial vehicles (drones), pressure sensors, and time-lapse cameras 56 . The ice wall collapse and the ensuing tsunami were filmed by tourists on a tour boat. The exact thickness of removed ice was determined at high spatial resolution (10 m pixel size) from TRI interferograms, as was the terminus geometry before, during and after the collapse (Lüthi and Vieli 10 ). The position of the tour boat was registered with the TRI, such that the tsunami wave height from the videos could be geometrically determined. The wave height in the impact zone on the shore opposite of the glacier was registered with a pressure sensor at a 4 s interval.
Concerning the 2D/3D transformation, the following ranges of F, S, and M have been used 10 in the article: Setup of the simulations. In all our simulations, the resolution of the MPM background mesh was chosen 40-50 times smaller than the minimum characteristic geometrical length of the system (rigid block length, glacier height, or water depth). In 2D and 3D simulations, we used approximately four and eight material points per element, respectively. In addition, the time step was evaluated based on both the CFL and elastic wave speeds to ensure the stability of the simulation.
Tsunami (Fig. 2). A 2D water tank of 16 m × 1 m is simulated. For GF and BF simulations, a solid block of 0.5 m × 0.5 m is used. In the GF experiment, the block is released on the left side of the tank with the bottom at the level of the water's free surface. For the BF case, the block is released with the bottom at a depth of −0.83 m corresponding to the experimental setup (Table 1). For the CS simulation, a solid block of 0.8 m × 0.25 m is used and released from a neutrally floating position with an initial angle of 35°to activate capsizing in a similar way as in the experiments (small initial movement by hand). For all cases, we use a block density of 920 kg/m 3 (similar to that of ice). The numerical and mechanical parameters are presented in Table 2.
Effect of buoyancy on iceberg lengths (Fig. 3). A 2D water tank of length 250 × N m and depth 50 × N + D w m is simulated with N ∈ {1, 3, 5, 10}. The ice slab has a length of 150 × N m and a cliff height of H i = N × 40 m. The ice slab slides on a frictionless boundary condition with a speed of 1 ffiffiffiffi N p m/s. Simulations are performed with a constant water level, i.e., the water level increase due to iceberg generation is compensated by water particles exiting the water tent at the water  Table 2 Numerical and mechanical parameters used in the simulations. tank outlet. The water tank has frictionless walls. The numerical and mechanical parameters are presented in Table 2. In this case, generic mechanical properties falling within a realistic range of reported values for ice 28 were chosen for the sake of the comparison with the analytical model. An elastic modulus of 1 GPa was chosen following Astrom et al. 5,20 .
Calving-induced tsunami at the Eqip Sermia glacier (Fig. 4). We simulated a 2D water tank of length 4500 m with variable depth. The first 500 m of the fjord have an average depth of 35 m and smoothly increases to a depth of~135 m in the main fjord, in agreement with field observations 10 . The ice cliff is 200 m high and the real-world shape of the calving ice mass measured in Lüthi and Vieli 36 has been modeled. The crevasse distribution has been simulated according to the Perlin simplex terrain noise model 35 leading to surface crevasses typically between 20 and 50 m deep and crack openings between 5 and 10 m. Parameters are presented in Table 2. In this case, model parameters were back-calculated within a range of realistic values for ice 28,42 to achieve good agreement (iceberg size and failure angle) with observations. An elastic modulus of 1 GPa was chosen following Astrom et al. 5,20 .
3D glacier calving (Fig. 5a). We simulated a glacier that is 1 km long, 2 km wide and has a 150 m irregular vertical cliff. The water tank is 200 m deep. The submergence depth is 40 m corresponding to H i /D~0.26. The glacier is supported by a frictionless base that slides backward at a speed of 3 m/s, in order to progressively remove basal support at the calving front and thus induce calving. Initial cracks have been placed within the glacier to mimic crevasses. These initial cracks follow a Voronoi distribution with cell size r = (30, 80, 160) m. Parameters are presented in Table 2. Generic yield surface parameters were chosen in agreement with experimental data 28 . Concerning the elastic modulus, most studies related to glacier calving numerical simulations use elastic moduli for ice lower than experimentally measured values to speed up simulation 5,20 (as the time step depends on the elastic wave speed). In this simulation, to test the capabilities of our model in simulating large-scale calving events with realistic properties, we used an elastic modulus of E = 10 GPa based on experimentally measured values 42 .
3D glacier calving and tsunami generation (Fig. 5b). We simulated a water tank of variable depth: the first 500 m are 300 m deep, and then the next 1250 m are 135 m deep. The ice slab is 600 m long and has a 200 m irregular vertical cliff. The crevasse distribution is the same as in the Eqip Sermia 2D simulation. The submergence depth is 50 m corresponding to H i /D~0.2. The glacier slides over a frictionless boundary condition at a speed of 3 m/s. Parameters are presented in Table 2. In this case, a low elastic modulus value was chosen to speed up the simulation (which suffers from the large size of the water tank). Generic values of the yield surface parameters were chosen within a realistic range 28,42 . An elastic modulus of 0.1 GPa was chosen following Astrom et al. 20 to speed up this simulation and focus on tsunami characteristics.

Data availability
The data corresponding to the iceberg-tsunami laboratory experiments can be found at https://hydralab.eu/research--results/ta-projects/project/hydralab-plus/11/. The data corresponding to the Eqip Sermia glacier calving-tsunami field experiments can be found in a previous publication at https://tc.copernicus.org/articles/10/995/2016/.