NUMERICAL SIMULATION OF TSUNAMIS GENERATED BY ICEBERG CALVING

Calving icebergs falling into water can generate large tsunamis, so-called iceberg-tsunamis. This phenomenon poses a threat for the fishing and shipping industries and coastal communities. Examples in Greenland include a 50 m amplitude wave recorded during an iceberg calving event at the Eqip Sermia glacier in 2014 and a tsunami generated by a capsizing iceberg in 1995 damaging a harbour. This work aims to numerically model the generation and propagation of such iceberg-tsunamis with the open source code Foam-extend 4.0. The original multiphase flow solver relying on the Immersed Boundary Method in this code has been modified to handle moving immersed boundaries, and it was then coupled with a motion solver to determine the iceberg motion. This method is validated with an own large-scale iceberg-tsunami test conducted in a 50 m × 50 m basin. The results show that the numerical iceberg motion and tsunamis generally agree with the laboratory observations. Further, the presented model can capture the laboratory wave heights and decay well. Future work will focus on simulating additional iceberg calving scenarios as well as modelling turbulence.


INTRODUCTION
Iceberg calving is the detachment of an ice mass from a larger ice volume such as a glacier or ice sheet. This phenomenon is a major reason for ice mass loss in Greenland and the Antarctica (Depoorter et al., 2013;Benn et al., 2017). Iceberg calving into water ( Figure 1) may generate large waves (Lüthi and Vieli, 2016). Such waves are called iceberg-tsunamis herein, short for iceberg-generated tsunamis, and can reach tens of meters in height (Heller et al., 2019c). Iceberg-tsunamis are generated by different iceberg calving mechanisms (Benn et al., 2007;Heller et al., 2019c), e.g. fall (icebergs fall into water), overturning (icebergs rotate around their base) and capsizing (floating icebergs capsize). Examples of iceberg-tsunamis in Greenland include a 50 m large wave observed at Eqip Sermia in 2014 (Lüthi and Vieli, 2016), a capsizing iceberg causing severe damage in a local harbour in 1995 (Mendsonboaz, 2009) and an event where the inhabitants of Innaarsuit had to be evacuated due to iceberg-tsunami hazards from a large iceberg floating offshore of the village in July 2018 (The Guardian, 2018). A number of recent studies confirm these potential hazards of iceberg-tsunamis (Levermann, 2011;MacAyeal et al., 2011;Burton et al., 2012;Lüthi and Vieli, 2016;Heller et al., 2019c). However, only a few field measurements and experimental studies have been conducted to quantify the generation and propagation of iceberg-tsunamis thus far. Experimental investigations include the small-scale laboratory wave flume tests of Burton et al. (2012) investigating iceberg calving and large-scale tests conducted by Heller et al. (2019a,b,c) in a 50 m × 50 m wave basin. The latter study addressed iceberg-tsunamis involving five iceberg calving mechanisms: A: capsizing, B: gravity-dominated fall, C: buoyancy-dominated fall, D: gravitydominated overturning and E: buoyancy-dominated overturning. Gravity-dominated icebergs essentially fall into the water body whereas buoyancy-dominated icebergs essentially rise to the water surface (Heller et al., 2019c). Lüthi and Vieli (2016) reported probably the best documented iceberg-tsunami field event thus far. They observed a 50 m large iceberg-tsunami near the Eqip Sermia glacier terminus which run-up of 10-15 m on the opposite shore at a distance of 4 km. Another field study was conducted by Vaňková and Holland (2016) to investigate iceberg-tsunami propagation through the Sermilik Fjord, Greenland, and they found that the height of one of the iceberg-tsunamis was still 24 cm at a distance of 30 km from the glacier terminus. They further used a numerical model to simulate the power spectrum density of the iceberg-tsunamis by using a prescribed method to simulate the generation process. This appears to be the only numerical simulation of icebergtsunamis thus far. Therefore, it is necessary to develop a numerical model to advance the investigation of such iceberg-tsunamis, which is the main aim of this work.
Since iceberg-tsunamis generated by the gravity-dominated fall mechanism are related to landslide-tsunamis (Heller and Hager, 2010;Heller and Spinneken, 2015), numerical models reproducing subaerial landslidetsunamis accurately may also be capable to simulate iceberg-tsunamis. Options are models based on the meshbased Euler method such as OpenFOAM, Thetis, REEF3D and the Lagrangian method Smoothed Particle Hydrodynamics (SPH). A weakness of SPH is the numerical energy dissipation, which often results in a strong overestimation of the wave decay for wave propagation (Violeau and Rogers, 2016) such that it is commonly coupled with a wave propagation model such as SWASH or XBeach to overcome this issue (Abadie et al., 2012;Tan et al., 2018;Ruffini et al., 2019).
The mesh-based Euler method is well capable of modelling fluid-structure interactions once the challenges of mesh adaptivity and free surface tracking are overcome. OpenFOAM is a common open source mesh-based computational fluid dynamics code containing numerous solvers and utilities to implement complex physical problems efficiently and flexibly (Jasak 2009). Therefore, OpenFOAM has the potential to reproduce the full range of iceberg calving mechanisms and iceberg-tsunamis. Additionally, OpenFOAM takes advantage of message passing interface parallelism, allowing single models to be run on multiple cores, which is convenient for large domains.
The handling of large displacements of bodies such as icebergs and the associated remeshing is a challenging key requirement in the context of iceberg-tsunami generation, which has not yet been fully resolved in OpenFOAM. Mesh-manipulation methods such as mesh deformation cause strong mesh distortion and topological changes may lead to an unusable mesh (Jasak et al., 2014). The Immersed Boundary Method (IBM) shows some key advantages over mesh-manipulation methods and avoids these problems: the boundary of the moving body is thereby represented by cells in the mesh at each time step, thus, the mesh itself does not need to be modified with the movement of the body.
An IBM toolbox has been implemented in Foam-extend (a branch of OpenFOAM) by Jasak et al. (2014), where the boundary of a solid body is represented by immersed boundary cells in the mesh as shown in Figure  2. This is achieved by recognising cells overlapped by the body boundary and then, based on their positions, to generate a mask to mark the fluid, solid and immersed boundary cells. The values of the variables of the immersed boundary cells, such as the velocity and pressure, are then evaluated by special interpolation schemes. The aim of this work is to set up and validate a numerical model to simulate both the generation and propagation of iceberg-tsunamis by coupling IBM with an (iceberg) motion solver in Foam-extend. To validate and calibrate this new approach, novel results of the large-scale tests of Heller et al. (2019a,b,c) are used. Section 2 provides the details about the numerical model, followed by a description of the experimental setup, numerical results and discussion in Section 3. The conclusions and future work are presented in Section 4.

Governing equations
In this numerical model, the motion of the iceberg and the surrounding water can be resolved respectively and coupled. The multiphase fluids are assumed to be incompressible, viscous and Newtonian. Under these conditions the governing equations are (adapted from Jasak et al., 2014) Here the two unknown variables and are the fluid velocity vector and pressure, respectively, and denotes the density, the dynamic viscosity and the gravitational acceleration vector. The Volume of Fluid method is applied to track the interface between the two fluids. The phase fraction (0 ≤ ≤ 1) is introduced with = 1 denoting one fluid (water), = 0 the other one (air), and 0 < < 1 the interface. The physical parameters, such as and of the two fluids, are then evaluated in function of as where the subscript and denote water and air, respectively. Once the velocity field is obtained, can be updated over time by solving the transport equation [5] In this work, the iceberg motions can be translation, rotation or a combination of the two and are flow-induced rather than prescribed. The equations of motion for the blocks are given as where and are the acceleration and angular acceleration vectors, respectively, is the force vector acting on the block, which is defined as where denotes the pressure force, the viscosity force caused by the two fluids, the drag force, −maa is the virtual force caused by the added mass and is the gravity force. is the total torque in relation to the centre of rotation, is the moment of inertia and denotes the mass of the iceberg. Details about the force calculation are presented in Section 2.2. and are directly calculated using field data in the current time step (pressure and velocity gradient). and are given as (Enet and Grilli, 2007) where is the drag force coefficient, denotes the iceberg main cross section perpendicular to the direction of velocity, denotes the submerged iceberg volume and is the added mass coefficient. The values of and are discussed in Section 3.2. Therefore, and in Eqs.
[6] and [7] are calculated as The subscript in Eqs.
[11] and [12] stands for the immersed boundary cell, denotes the vector of the immersed boundary cell area and is the shear stress along the immersed boundary cell which can be calculated by multiplying the immersed boundary cell's dynamic viscosity and its gradient of velocity. In Eq. [12], is the vector pointing from the immersed boundary cell to the centre of rotation and and are the position vector of the centre of mass and rotation of the block, respectively.

Newly implemented solver in Foam-extend
A new solver has been implemented based on interIbFoam in Foam-extend 4.0, which is an already available incompressible multiphase fluid solver with IBM support. In contrast to interIbFoam based on a static mesh, the newly implemented solver interDyMIbFoam within this work, can handle dynamic immersed boundaries in order to describe various kinds of motion of moving bodies, such as icebergs.
In order to determine the iceberg motion, a motion solver is required and should be coupled with the fluid solver. This motion solver is a modification of the already provided solver sixDoFMotion. Note that not all 6 Degrees of Freedom (DoF) motions are required in the iceberg-tsunami mechanism modelled within this study. Some restrictions are applied by resetting some moment and force components to zero at each time step for the translational and rotational motion, respectively, to replicate the laboratory conditions.
The two-way coupling method between the fluid solver and motion solver is achieved by a new dynamic mesh handling class, via which the velocity and pressure field data are read and used to calculate the new force. In turn, the new position of the immersed boundary may change the velocity and pressure field according to the no-slip condition. Figure 3 shows the steps applied in the new interDyMIbFoam solver. When the solver is executed, a small initial time step t 1 is set which is then controlled by the Courant-Friedrichs-Lewy number. The forces in Eqs.
[6], [7] and [8] on the iceberg are calculated for each time step before the motion solver is called to determine the new position of the block. The immersed boundary is then updated by regenerating the immersed boundary mask. This step is analogous to moving meshes in other dynamic mesh solvers, where the face fluxes are adjusted. Thereafter, governed by the PIMPLE loop algorithm (Issa, 1986), interDyMIbFoam solves the velocity and pressure equations to obtain the velocity and pressure fields successively, and then Eq.
[5] to reconstruct the current interface between the two fluids. Finally, a turbulence correction function could be called for each time step. However, turbulence is not considered herein. This may not be fully appropriate for all iceberg-tsunami experiments conducted by Heller et al. (2019a,b,c), as further discussed in Section 3.2.

Experimental setup
Large-scale experiments have been conducted in a 50 m × 50 m large basin with an effective area of 40.3 m × 33.9 m. While the iceberg-tsunamis were generated by five different iceberg-calving mechanisms (Heller et al., 2019a,b,c), only the gravity-dominated fall mechanism (Figure 4) is used herein to calibrate and validate the numerical model.
A purpose-built steel frame to support the iceberg block was fixed to the basin side wall. A block made of polypropylene homopolymer with a density ≈ 920 kg/m 3 was used to mimic an iceberg. The block was hold in position with an electromagnet prior to release. This electromagnet was attached to a steel plate integrated into the block. Figure 4(a) shows the side view of the gravity-dominated fall experiment. The water depth h was 1.00 m and the basin bottom was horizontal. The block length l, width b and thickness s were 0.500, 0.800 and 0.500 m, respectively. The front release position in Figure 4(a), corresponding to the distance of the bottom face of the block from the still water surface, was 0.0 m. The iceberg-tsunami features were measured with 35 resistance-type wave probes with a sampling frequency of 100 Hz. The probes were only placed in a quarter circle as shown in Figure 4(b), because the wave field was symmetric relative to the block axis. A cylindrical coordinate system (r, z, ) has been adopted with the origin located at the steel frame centre on the water surface (Figure 4). The wave propagation angle (angular angle) is defined positive in clockwise direction. The locations of the wave probes are shown in Table  1, together with the location of the 5 MP camera used for general observations. A 9-Degree of Freedom (DoF) motion sensor was fixed on the top face of the block to record the block motion. The 9-DoF motion sensor measured accelerations along three local axes, three global angles and three components from the Earth's geomagnetic field. The time array was shifted to a common starting point for the motion sensor, wave probes and the camera via a synchronisation system. Low-pass filters with cut-off frequencies between 3 and 11 Hz were applied to remove noise in the wave probe data (Heller, 2019). Details about the trajectory inference method to extract the block velocity and position based on the 9-DoF motion sensor can be found in Chen et al. (2019). Table 1. Locations of the wave probes and camera in the laboratory experiment. Wave probes marked with * were also used in the numerical wave basin (Heller, 2019).  The convergence tests were conducted in a smaller wave basin than shown in Figure 5, namely in a basin of length and width of 6.0 m. Note that the resolution plays an important role in the force calculation, affecting both the iceberg motion and iceberg-tsunamis. To keep the motion constant with different resolutions in the convergence tests, the motion was prescribed with the motion measured in the laboratory experiments. Figure  6 shows the waves obtained in the convergence tests at wave probe B1 together with the laboratory results. The convergence tests show that the wave profiles reasonably well collapse on a curve for resolutions higher than 5.0 cm × 2.5 cm × 2.5 cm. Therefore, a cell dimension of 2.5 cm × 2.5 cm × 2.5 cm was chosen for the main tests for both the iceberg-tsunami generation and propagation zones.

Device
The main tests were based on resolved iceberg motion. The computational domain shown in Figure 5 contained 20,768,000 cells and required 54 cores, and the main test took approximately 50 h to complete. For the gravity-dominated fall case, two parameters by given resolution affect the numerical results: the drag force coefficient and the add mass coefficient . Some indications for the values of and are available from Enet and Grilli (2007) who simulated rigid submarine landslide-tsunamis with = 0.36 and = 0.61. Figure 7 shows the vertical iceberg displacements for different and values, and the corresponding iceberg-tsunamis recorded at wave probe B1 are shown in Figure 8. As expected, increasing and reduces the iceberg motion and tsunami heights, and the waves occur later in time. The best agreement between the numerical and experimental iceberg movements and tsunamis are obtained for = 0.6 and = 0.4.   The first wave amplitude and wave height are well captured, and the difference between the simulated and measured wave amplitude and wave height are 15.5% and 12.6% (Table 2), respectively, for Cd = 0.6 and Cm = 0.4 (Figure 8). The reason for this underestimation may be that the numerical acceleration of the iceberg in Figure 7 is larger than in the experiment when the block is fully submerged, making the iceberg moving back to the water surface faster and inhibiting the growth of the first wave. This is related to the IBM method where the boundary of the iceberg is represented by cells, resulting in a slightly larger iceberg than the real geometry increasing the buoyancy force.
Wave decay during iceberg-tsunami propagation is important for hazard assessment for off-and onshore structures. Figure 9 shows the time history of the water surface elevation η at wave probes B1, B7 and B13 from both the numerical and laboratory models. In order to calculate the decay of the wave amplitude a and height H, the relative values of a, H and r are used, i.e. a/h, H/h and r/h. Based on the assumption that the waves decay with a power function of the form ( /ℎ) − , the wave decay exponent c comparison of the numerical and laboratory model is shown in Table 2, together with a/h and H/h of the first wave. Table 2 shows that both the numerical wave amplitude and wave height decay agree well with the laboratory tests, especially the wave height decay showing only a 2.2% deviation in the exponent of the wave decay power function.

CONCLUSIONS
This article presented a numerical model to simulate the generation and propagation of tsunamis generated by calving icebergs, so-called iceberg-tsunamis. This numerical model is based on the Immersed Boundary Method (IBM) in Foam-extend. The multiphase solver interIbFoam has been coupled with a newly implemented motion solver. This enables to handle dynamic immersed boundaries to resolve iceberg motion under a wide range of iceberg calving mechanisms. Coupling between the motion and flow solvers is achieved by simulating the fluid-solid interaction including the calculations of pressure force, viscose force, drag force and virtual force due to the added mass. The wave basin dimensions of 15.0 m × 18.0 m × 1.7 m allowed the first three waves measured up to 5 m from the block impact location to be recorded without reflections. In the laboratory model, the movement of a 0.500 m × 0.800 m × 0.500 m large "iceberg" was recorded with a 9-DoF motion sensor, and the iceberg-tsunamis were measured with 35 resistance-type wave probes up to a distance of 35 m from the iceberg impact location.