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

Glacier calving is a major contributor to sea-level rise, accounting for roughly half of the mass loss from the Greenland and Antarctic ice sheets. Calving arises from surface crevasse formation and subsequent fracture processes driven by stress concentrations, ice stretching near the front, and oceanic erosion and buoyancy-induced torques, with scenarios dependent on glacier outlet geometry. Calving can generate large, hazardous tsunamis impacting coastal infrastructure, ecosystems, and populations, and in mountain lakes can trigger outburst floods. Despite its importance, modeling calving and its tsunamigenic potential at event timescales remains challenging due to complex, rapid solid–fluid interactions. This study addresses the need for a unified, physics-based framework capable of simulating dynamic ice fracture and the generation and propagation of calving-induced tsunamis, and validating against laboratory and field observations.

Most marine-terminating glacier models focus on slow viscous creep using Eulerian continuum approaches (e.g., Elmer/Ice), estimating calving rates from simplified empirical laws with limited validation and little focus on single-event dynamics. Lagrangian particle-based models (Åström et al.; Bassis and Jacobs) simulated sea ice and calving dynamics and captured debris size distributions and geometry-dependent calving features but omitted explicit water modeling and are computationally expensive, limiting tsunami studies. Transient multiphysics finite element models (Mercenier et al.) coupled ice flow with damage-based criteria and realistic front geometries but did not simulate calving-induced tsunamis. Landslide-generated tsunami research is extensive, yet calving-induced tsunami studies are fewer: Lüthi and Vieli documented a 45–50 m tsunami at Eqip Sermia; large-scale experiments by Heller et al. showed that landslide-based empirical equations overestimate calving-induced wave amplitudes; Chen et al. reproduced laboratory wave characteristics using immersed boundary CFD. A unified approach simulating both dynamic glacier fracture and tsunamis has been lacking.

The authors develop a continuum damage Material Point Method (MPM) framework to couple dynamic ice fracture with water dynamics and tsunami generation. Ice mechanics: A non-associative elastoplastic constitutive model based on the Cohesive Cam Clay (CCC) yield surface is adopted, suitable for brittle fracture of low-porosity ice. The yield surface depends on pressure p and Mises equivalent stress q, with cohesion parameter β, critical state slope M, and pre-consolidation pressure p0. A non-associative plastic flow rule ensures near volume preservation in plasticity. Softening/hardening is modeled by evolving p0 using a deviatoric plastic strain-based hardening variable; a q-based non-associative return mapping projects states to the yield surface, with softening for pp0. Elastic response uses a split Neo-Hookean hyperelasticity (deviatoric–dilational split) to accommodate large deformations. Water mechanics: Water is modeled as a nearly incompressible fluid via a weakly compressible equation of state penalizing density deviations, commonly used in SPH, to approximate incompressibility efficiently within MPM. Numerical solution: MPM discretizes continua into Lagrangian material points and solves mass and momentum conservation on a background Eulerian grid using an explicit time-stepping scheme (symplectic Euler) with APIC transfers to better conserve linear and angular momentum. Stability is enforced via CFL and elastic wave speed constraints. Validation and simulation setups: (1) Laboratory-scale tsunami generation: 2D tanks (16 m × 1 m) simulate three calving mechanisms from Heller et al.: gravity-dominated fall (GF), buoyancy-dominated fall (BF), and capsizing (CS), using rigid blocks with ice-like density. Simulated wave amplitudes are compared to experiments using a 2D-to-3D amplitude transformation (Heller & Spinneken) accounting for Froude number, relative slide thickness, and mass, with reported uncertainty. Wave celerity is compared to linear shallow-water theory v=√(gD). (2) Iceberg size vs. buoyancy: A sliding ice block of height H interacts with water at submergence depth D to vary buoyancy; dynamic fractures and iceberg lengths are measured and compared to an analytical bending model derived from beam theory, yielding L/H = 2σ_t / [3 ρ_i g H (1 − D/H)]. (3) Real-world case: Eqip Sermia (Greenland) 2D simulation uses observed fjord bathymetry (35–135 m depths) and a 200 m ice cliff. The calving front geometry is initialized from observations; surface crevasses are mimicked using simplex Perlin noise to emulate crevasse depth/opening statistics. Model parameters (within laboratory-measured ice property ranges) are calibrated to match observed iceberg geometry and failure angle. Wave amplitudes are measured at specified distances and converted to 3D ranges with the 2D/3D transformation for comparison to field measurements. Wave speeds are compared with shallow-water theory for the fjord’s depth. (4) 3D demonstrations: Large-scale 3D simulations (20–36 million particles) explore 3D fracture patterns and tsunami propagation with realistic geometries and crevasse fields, using frictionless basal sliding to induce calving. Computational performance is reported (<1 day on a workstation). Data and code: Laboratory and field datasets are publicly available; the CD-MPM code is open-source on GitHub.

• Laboratory-scale tsunami generation: The model reproduces the temporal evolution and maxima of primary wave amplitudes for three calving mechanisms (GF, BF, CS), with good agreement to Heller et al.’s experiments after applying the 2D/3D transformation. Maximum amplitudes are matched for GF and CS and slightly overestimated for BF. Wave fronts travel at velocities consistent with linear shallow-water theory v≈√(gD) and in fair agreement with experiments. Consistent with experiments, gravity-dominated falls generate the largest tsunami amplitudes, indicating highest hazard potential.
• Iceberg size dependence on buoyancy: In sliding-block calving simulations, the first iceberg length L increases with submergence depth D/H due to buoyant stabilization, reaching no-failure near D/H≈0.92 (≈ρ_i/ρ_w), then decreases for deeper submergence where buoyancy overcomes body forces and fracture initiates from the base. Measured L/H agrees with the analytical bending model L/H = 2σ_t / [3 ρ_i g H (1 − D/H)].
• Eqip Sermia case study: Simulated calving geometry matches observations: failure width 75 m (observed 50–60 m) and failure plane angle 54° (observed 45–60°). Simulated wave amplitude at 250 m from the front is 50 m (observed 45–50 m). At the tide gauge distance on the opposite shore, the 2D simulation yields 27 m amplitude, corresponding (via 2D/3D conversion) to 2.5–6.6 m, consistent with the measured 3.3 m. Average simulated wave speed is 35 m/s, matching the measured 32 m/s and theoretical shallow-water predictions for fjord depths 100–150 m (31–38 m/s).
• Computational feasibility: 3D calving and tsunami simulations with 20–36 million particles run in under a day on a workstation, enabling practical hazard assessment studies.

The study demonstrates that a unified continuum damage MPM framework can capture the coupled dynamics of glacier calving and tsunami generation. Validations against controlled laboratory experiments confirm that the model reproduces key wave amplitude and celerity characteristics for distinct calving mechanisms, with discrepancies largely attributable to the empirical 2D/3D amplitude transformation uncertainty and expected deviations from linear shallow-water theory due to nonlinearity and dispersion. The iceberg size analysis reveals the controlling balance between buoyancy and ice weight, identifying a submergence threshold near D/H≈0.92 where bending stresses are neutralized, and validates against a simple bending stress analytical model. Application to Eqip Sermia shows quantitative agreement with observed failure geometry, tsunami amplitude decay over kilometers, and wave speeds, indicating the model’s capacity to simulate real-world events. Collectively, these results address the need for event-scale, multiphysics simulations of calving and tsunamis, providing a tool to interpret observations, refine empirical calving laws for Earth system models, and inform coastal hazard assessment and mitigation.

This work introduces a continuum damage MPM that unifies dynamic glacier fracture with tsunami generation and propagation, validated against large-scale laboratory experiments and a well-observed calving event at Eqip Sermia. The model reproduces mechanism-dependent tsunami amplitudes and speeds, predicts iceberg sizes consistent with an analytical bending model, and matches field observations of wave heights and velocities. Its computational efficiency and geometric flexibility make it suitable for hazard assessment of calving-induced tsunamis and related outburst events. Future research should (i) couple the present brittle-fracture MPM with continuum glacier flow models to account for viscous creep and basal processes over longer timescales; (ii) incorporate more detailed water dynamics (e.g., incompressible Navier–Stokes) where feasible to improve wave predictions; (iii) represent anisotropy and spatial variability in ice mechanical properties and explicit crevasse fields to better capture iceberg size distributions; and (iv) extend to fully 3D, site-specific hazard studies with data assimilation.

• Water modeling uses a weakly compressible approximation rather than fully incompressible Navier–Stokes, potentially limiting accuracy for complex wave dynamics.
• The framework targets short timescales with brittle ice behavior; it does not simulate viscous creep or detailed basal processes that form crevasses.
• Ice mechanical properties are assumed within realistic ranges but do not capture the full variability and anisotropy of glacial ice, which can affect fracture and iceberg size distributions.
• Comparisons of 2D simulations to 3D observations rely on an empirical 2D/3D amplitude transformation with reported ±50% uncertainty.
• Wave speed comparisons may deviate from linear shallow-water theory due to nonlinearity and dispersion, especially for intermediate-depth, finite-amplitude waves.
• Certain numerical simplifications (e.g., stiff-material assumption in the plastic flow rule enabling analytical return mapping) may not generalize to very soft materials.