# Rheology and bathymetry effects captured by a multilayer landslide-tsunami model

A Marbœuf, Anne Mangeney, Anne Le Friant, Enrique Domingo Fernández-Nieto, Antoine Lucas, P Poulain, Annabelle Moatty, M Silver, Rodrigo Pedreros, Anne Lemoine, et al.

► To cite this version:

A Marbœuf, Anne Mangeney, Anne Le Friant, Enrique Domingo Fernández-Nieto, Antoine Lucas, et al.. Rheology and bathymetry effects captured by a multilayer landslide-tsunami model. Geophysical Journal International, 2025, 10.1093/gji/ggaf232. hal-05134424

HAL Id: hal-05134424

https://hal.science/hal-05134424v1

Submitted on 28 Jun 2025

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
# Rheology and bathymetry effects captured by a multilayer landslide-tsunami model

A. Marbœuf¹, A. Mangeney¹, A. Le Friant¹, M. J. Castro², E. Fernández-Nieto³,
A. Lucas¹, P. Poulain¹, A. Moatty⁴, M. Silver⁵, R. Pedreros⁶, A. Lemoine⁶, F. Løvholt⁷

¹ Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, F-75005 Paris, France

² Departamento Análisis Matemático, Estadística e Investigación Operativa y Matemática Aplicada.
Universidad de Málaga, Campus Teatinos S/N, Spain.

³ Departamento de Matemática Aplicada I, Universidad de Sevilla.E.T.S. Arquitectura.
Avda, Reina Mercedes, Sevilla s/n. 41012, Spain

⁴ UMR LGP (Laboratoire de Géographie Physique) - CNRS, Paris 1 Panthéon Sorbonne, UPEC

⁵ US Department of Agriculture (USDA) Natural Resources Conservation Service (NRCS)

⁶ BRGM, French Geological Survey, Direction Risks, Orléans, France

⁷ Norwegian Geotechnical Institute,
Sandakerveien 140, 0484 Oslo (Street Address), P.O. Box 3930 Ullevl Stadion, NO-0806 Oslo, Norway

Received tbd; in original form tbd

## SUMMARY

Numerical models are a powerful tool for understanding and predicting the impact of landslide-generated tsunamis. We use here the HySEA code which incorporates a multilayer structure and non-hydrostatic pressure to simulate the tsunami generated by a potential submarine landslide located offshore of the Mayotte island. The island is surrounded by a lagoon and steep slopes weakened by the 2018 seismo-volcanic crisis. The influence of the input parameters and of the model assumptions is shown to change by a factor 2 the predicted maximum water
free surface elevation, velocity and maximum inundation depth. This demonstrates the need of using numerical models for building local scaling laws to relate tsunami and landslide properties. Our results highlight the necessity of incorporating high-resolution bathymetry, in-depth variations through multilayer modeling and relevant landslide rheology to accurately predict tsunami impact. In case of strong topography variations as in Mayotte, using 4 layers seems to be a good compromise between accuracy and computational cost. Accounting for these effects would enable to refine hazard maps by identifying safe and high-risk coastal zones and to improve wave arrival time estimates, thus reducing tsunami-related risks in regions like Mayotte.

**Key words:** Numerical modeling, Submarine landslides, Tsunamis

# 1 INTRODUCTION

Neumann et al. (2015) estimate that 625 million people live in low-elevation coastal zones, i.e. altitudes below 10 m. That number is expected to increase in the near future and reach 1 billion by 2030, mostly concerning population in least developed countries. Tsunami risk is growing due to climate change and the denser concentrations of people living in these area. While earthquakes are responsible of about 80 % of all historical tsunamis, landslides come second with 5 % as reported by Harbitz et al. (2014). In addition, when landslides are not the sole source of tsunamis, they are also reported as a contributing factor in a further 13 % jointly with other sources, sometimes as the main tsunamigenic source. Hence, the landslide contribution is often strong when both earthquakes and landslides are involved. Identifying and assessing landslide-induced tsunami hazards is essential for reducing the risk of catastrophic events on communities. Such assessments often involves hazard or probability maps, which inform early-warning systems, evacuation plans, and the development of natural or engineered coastal defenses.

## 1.1 Landslide tsunami hazard assessment

Landslide-induced tsunami hazard identification and assessment is a challenging task due to its multidisciplinary aspect (Roger et al., 2024). This paper focuses on one of the major component of tsunami hazard identification and assessment: numerical modeling of landslide-generated
tsunamis. Calibrated towards laboratory experiments (e.g. (Grilli et al., 2002; Ma et al., 2015; Poulain et al., 2023)) and past landslides and tsunamis in the field (e.g. (Crosta et al., 2016; Gylfadóttir et al., 2017; Brunet et al., 2017; Paris et al., 2019; Svennevig et al., 2024)), numerical models are a powerful tool capable of quantifying tsunami height and velocity over a realistic bathymetry for a given landslide scenario. Yavari-Ramshe and Ataie-Ashtiani review the available approaches in the literature and show the large variety of existing models (Yavari-Ramshe & Ataie-Ashtiani, 2016, Table 3). Each model implements different degrees of complexity when resolving the water part: (i) 2D and 3D hydrostatic models with no dispersion effect that include non-linear Shallow Water equations (Harbitz, 1992; Allgeyer et al., 2019; Poulain et al., 2022); (ii) Boussinesq or one-layer weakly dispersive models (Lynett & Liu, 2002; Aïssiouene et al., 2020); (iii) Multilayer non-hydrostatic dispersive models (Macías et al., 2021a,b); (iv) Reynolds-Averaged Navier-Stokes (RANS) (Liu et al., 2005) and Large Eddy Simulations (LES) (Mokhtarzadeh et al., 2021) models; and (v) full Navier-Stokes models (Mangeney et al., 2000; Abadie et al., 2010; Allgeyer et al., 2019). Note that the performances of these models against observations may be evaluated thanks to techniques such that the Nash–Sutcliffe model efficiency coefficient (NSE) (Mathevet et al., 2023). These models are implemented in numerical codes that are currently in use and undergo continuous development: see e.g. FUNWAVE (Kirby et al., 1998), Basilisk (Popinet, 2003), HySEA (Macías et al., 2015, 2021a,b), OpenFOAM (Rauter et al., 2022) or SWASH (Ruffini et al., 2021). The reader is referred to Kirby et al. (2022); Yavari-Ramshe & Ataie-Ashtiani (2016); Marras & Mandli (2021) for more details. Concerning landslide modeling, the increasing complexity in the models aim at reproducing the behavior of the mixture between the granular material and the fluid. More accurate models include advanced rheological formulations concerning particle interaction and pore fluid stresses (e.g. (Rauter et al., 2022)). One finds in the literature: (i) depth-averaged models (e.g. (Fernández-Nieto et al., 2008; Grilli et al., 2017; Delgado-Sánchez et al., 2020; Gueugneau et al., 2021)); (ii) Navier-Stokes models (e.g. (Heinrich et al., 1999; Rauter et al., 2021)); and (iii) particle-sized models like Smoothed Particle Hydrodynamics (SPH) or Discrete Element Method (DEM) (e.g. (Guan & Shi, 2023; Radjai & Dubois, 2011)). A comparison between these different methods is performed by Martin et al. (2023).
## 4 *A. Marbœuf*

Regardless of their complexity, numerical models have well-known limitations, especially when simulating real-world events. The two primary challenges stem from limited computational power (see e.g. (Windows-Yule et al., 2015)) and gaps in data or understanding of the underlying physical processes. Masson et al. (2006) emphasize that the reliability of numerical outputs are sensitive to the quality of available field data. In particular, field data related to a geophysical description of the underburden, including different geophysical horizons, geotechnical properties of soils or composition and stability of joints and faults (in the case of rock slopes), are all of first order importance. The generation and release probability mechanism also involves a high degree of uncertainty that should be embedded in prognostic forecasting (Løvholt et al., 2020). Bathy-metric data is another important source of uncertainties (Hare et al., 2011; Svennevig et al., 2024) since high resolution data (lower than 1 m) are rarely available (Cannat et al., 2013). Probabilistic and statistical methods are well suited for dealing with these uncertainties. Several groups have attempted sensitivity analyses (Gonzalez-Vida et al., 2019; Sabeti & Heidarzadeh, 2022; Tinti et al., 2000; Esposti Ongaro et al., 2025) sometimes coupled with advanced probabilistic / statistical techniques (Sánchez-Linares et al., 2015; Zhao et al., 2021; Behrens & Dias, 2015; Mulia et al., 2022; Tozato et al., 2023) to investigate the uncertainties related to input data. For earthquake tsunamis, Probabilistic Tsunami Hazard Assessments (PTHA) approaches have become a standard (Grezio et al., 2017, and reference therein). For landslide tsunamis, probabilistic methods are less developed, but have been employed for local settings where there are rich availability of past landslide data (Lane et al., 2016; Esposti Ongaro et al., 2025) or detailed assessment of local slope stability (Løvholt et al., 2020, 2025). To this end, improved modeling capabilities for the tsunami generation will reduce modeling uncertainty.

## 1.2 Context and background on Mayotte

The volcanic island of Mayotte is located in the Indian ocean and is composed of Petite Terre and Grande Terre. The island catches the attention of researchers and the French government since the beginning of an intense seismo-volcanic crisis in 2018 (Mercury et al., 2022). The location of earthquakes near steep slopes and the construction of a new volcanic structure (Feuillet et al.,
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Mayotte bathymetric data and computational domain considered in this study. The two islands Petite Terre and Grande Terre are surrounded by a well-developed shallow submarine shelf defining a lagoon. The zoom on the right indicates the sensible areas in terms of strategic infrastructures.</figcaption></figure>

2021) may trigger submarine instabilities, in particular east of Mayotte (Roger, 2019). The two islands Petite Terre and Grande Terre are surrounded by a well-developed shallow submarine shelf defining a lagoon (Figure 1). The shelf corresponds to a significant topography change, from shallow slopes on the shelf (< 9°) to flanks with maximum local slopes of 25° to 60° (Lemoine et al., 2020b; Poulain et al., 2022). The flanks are prone to tsunamigenic gravitational instabilities given the repetitive seismic activity located 5 to 15 km from the coast east of Petite Terre. In deeper water, gullies and canyons form large valleys may control the circulation of sediments. Studies on Mayotte concern mainly the volcanic and seismic activities (Roger, 2019; Thinon et al., 2022; Mercury et al., 2022, 2020; Lemoine et al., 2020a; Saurel et al., 2021; Sultan et al., 2023). After a preliminary overview considering 62 tsunami scenarios performed by Lemoine et al. (2020b), both of seismic and non-seismic sources, Poulain et al. (2022) focus on the modeling of the tsunami impact related to submarine landslides. The authors couple two numerical codes: hydrostatic HYSEA (Fernández-Nieto et al., 2008; Macías et al., 2020; Escalante et al., 2019) and FUNWAVE-TVD (Wei et al., 1995; Chen, 2006; Kennedy et al., 2001; Shi et al., 2012). They test hydrostatic/Boussinesq (i.e. non dispersive and weakly dispersive) numerical models and the influence
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Bathymetry and landslide sources from Lemoine et al. (2020b). Circled submarine landslide scenarios are those considered by Poulain et al. (2022).</figcaption></figure>

of the initial landslide volume through 8 different submarine scenarios. However, this coupling introduces additional errors and reduces the computational efficiency. We will focus here on one of the most threatening submarine scenario, identified by Lemoine et al. (2020b) and considered by Poulain et al. (2022), using only the multilayer code HySEA describing the landslide dynamics as well as tsunami generation and propagation.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** (Left) Initial released mass, represented in red, for the considered scenario called Piton 200. (Right) Cross section (AA') of the bathymetry and of the initial released mass.</figcaption></figure>
**1.3 Motivation and objectives**

While a few landslide-tsunami studies have included in-depth variation in the models using mul-
tiple layers in the vertical direction (Macías et al., 2017; Pedrosa-González et al., 2022; Bonilauri
et al., 2024), few have numerically tested the influence of the number of layers, with the notable
exception of Ongaro et al. (2021). The present paper builds on Poulain et al. (2022), delving
deeper into the impact of input parameters. Leveraging recent advancements in HySEA (Macías
et al., 2021a,b), we conduct a robust sensitivity analysis focusing on newly introduced features
controlling the dispersive properties of the model: multilayer structure and non-hydrostatic pres-
sures. Since these new features are computationally demanding, evaluating their impact on the
final results is crucial to determine whether their inclusion is necessary. In particular, are these
new features of great influence on the computed waves hitting the Mayotte's coasts? We are also
interested in the influence of crucial parameters such as the grid resolution or the landslide related
parameters. This sensitivity analysis will greatly help in the results interpretation and the choice of
many input parameters, given the large uncertainties presented in the literature. As a first attempt
to investigate the sensitivity to input parameters on Mayotte simulations, we restrict ourselves to a
single volume configuration, such as in Gonzalez-Vida et al. (2019). We also keep statistical and
probabilistic techniques for next studies.

Although landslide-tsunami models are currently employed for hazard assessment and risk
analysis (Leone et al., 2011), significant uncertainties remain unquantified. These uncertainties
are highly dependent on factors specific to the simulated event, such as topography, material prop-
erties, landslide volume, and water depth. Our objective here is to quantify the variability of the
sea elevation, runup and tsunami impact velocity induced by: (i) the hydro versus non-hydrostatic
approximation; (ii) the precision of the flow description in the vertical direction (from depth-
averaged to multilayer equations); (iii) the bathymetry resolution; and (iv) the landslide rheology.
The Mayotte configuration is a generic case study involving a coral reef, which often happens in
oceanic islands. The coral reef induces strong variations of the bathymetry and related effects on
the tsunamigenesis.

The rest of the paper is organized as follows: Mayotte's field data and the selected scenario are
8 *A. Marbœuf*

presented in Section 2; the multilayer HySEA code is introduced and its Shallow Water models for
both the landslide and the water waves are described in Section 3; the reference case is presented
in Section 4; the sensitivity analysis is performed in Sections 5, 6 and 7; this work is concluded in
Section 8.

## 2 SCENARIO AND FIELD DATA

### 2.1 Mayotte's bathymetric and sea floor data

The accurate description of the bathymetry is known to have a major influence on tsunami propa-
gation. Bathymetric data are shown in Figures 1 and 2 and are taken from Lemoine et al. (2020b).
The data are based on Gebco 2014 (https://www.gebco.net), HOMONIM SHOM DTM (100 m
resolution, https://data.shom.fr), MAYOBS 1 (Feuillet et al., 2021) (30 m resolution), bathymetric
surveys of SHOM (25 m resolution, https://data.shom.fr), and litto3D (lidar data at 1 m resolu-
tion, https://data.shom.fr). All these data are gathered to build a 10 m resolution Digital Elevation
Model (DEM) on the computational domain shown in Figure 1. The coordinate system used and
displayed for all the figures is EPSG:4471-RGM04/UTM zone 38S. The Histolitt coastline from
SHOM is represented on all the figures by a bold black full line. The main specificity of this
bathymetry is the steep slope at the coral reef with typical water depths of 1000 m east of the reef
and 30 m in the lagoon (Figure A1).

Figure A2 shows the spatial distribution of the Manning coefficient involved in the bed-water
friction law (Section 3.1). The Manning coefficient is based on sea floor morphology as character-
ized in the 'Corine Land Cover' database and local coral reef (Lemoine et al., 2020b). Given the
condition of the sea floor, a corresponding Manning coefficient has been taken from the literature
(Bunya et al., 2010). Note that the Manning coefficient is higher along the topography change at
the border of the coral reef compared to its value elsewhere below the ocean.

### 2.2 Selected scenario

A key challenge in landslide-tsunami hazard assessment is constructing realistic landslide scenar-
ios. This difficulty arises primarily from limited data on tsunamigenic landslide volumes, locations,
occurrence rates, digital elevation models, and gaps in understanding landslide dynamics and material behavior (Roger et al., 2024; Løvholt et al., 2020). In their exploratory study, Lemoine et al. (2020b) simulated tsunamis generated by 32 submarine landslide scenarios around Mayotte. The authors identify the most threatening events to the local population and infrastructures from: (i) initial sea surface deformation modeled using TOPICS (Tsunami Open and Progressive Initial Conditions System) tool; (ii) the empirical GEOWAVE model (Watts et al., 2003), based on well-tested parameter values. Using more advanced landslide-tsunami simulations, Poulain et al. (2022) build and analyze 8 submarine landslide scenarios from the most impactful events with diverse volumes and depths (circled in Figure 2), providing hazard maps to the local French authorities. However, at that time, they use the hydrostatic depth-averaged version of HySEA coupled with the wave propagation code FUNWAVE-TVD.

Using a single and more precise numerical model, we focus our study on one submarine scenario called Piton 200 selected for the potential threat it represents to the local population and infrastructures (Poulain et al., 2022; Lemoine et al., 2020b). The collapse structure and extent of Piton 200, described in Figure 3, are defined by a geomorphological analysis of bathymetric surveys. The collapse structure is located at the shelf to slope transition close to Petite Terre (2.5 km) at depths between 50 and 600 m below sea level and it involves a volume of 200 Mm³. The Piton 200 scenario meets three of the four hazard characteristics of submarine landslide tsunami cited by Roger et al. (2024): (i) rapid arrival time or short tsunami travel time, (ii) relatively large wave amplitudes, and (iii) lack of forewarning.

## 2.3 Gauge locations

Figure 2 presents the locations of the six virtual wave gauges used in this study. Gauge n° 3 is positioned at the Piton 200 location (76.5 m depth) and gauge n° 4 is located 5 km further offshore to provide wave observation in front of the landslide and at greater water depth (1106 m depth). Three additional gauges are put closer to the coast. These gauges have been positioned around major assets for the island: the airport runway, the terminal for the ferry that links Petite Terre and Grande Terre, and the DEALM (Direction de l'Environnement, de l'Aménagement, du Logement et de la
Mer). The airport gauge n° 1 is located in the commune of Pamandzi at the end of the current runway (3.23 m depth). The gauge n° 6 is placed 2.3 km east of the airport gauge (same Y coordinate as the airport gauge) on the other side of the reef (466 m depth). This gauge corresponds to the location where the vertical velocity profile is animated: *animationEastAirportVelX.mp4* available in the supplementary material. The gauge n° 6 also defines the cross section at constant Y coordinate with the airport gauge: Figures A1, animations *animation1Layer.mp4* and *animation6Layer.mp4* available in the supplementary material. The pier's gauge n° 5 (Ponton) at 9.8 m depth is located in Dzaoudzi and also represents a critical infrastructure asset, as commuter traffic between Petite and Grande Terre is a daily occurrence. Finally, the gauge n° 2 on Grande Terre (DEALM) is located at 2.5 m depth near the DEALM building at M'tsapéré, south of Mamoudzou. These gauges allow a more local analysis and a closer look to the results at the selected critical locations. Only the 4 first gauges Airport, DEALM, Piton200 and Offshore are included in the main text. For convenience, the gauge data for gauge n° 5 (Ponton) are shown in Appendix D. Table A2 of Appendix D provides the exact coordinates (m) and the depth (m) of these 6 gauges.

## 3 NUMERICAL MODEL

Multilayer HySEA (Hyperbolic Systems and Efficient Algorithms) is part of a family of long wave landslide-tsunami models (Macías et al., 2021a,b). HySEA 1-layer models have been widely used and tested against laboratory experiments (e.g. (Poulain et al., 2023)) and field scale events (e.g. (Linares-Sánchez et al., 2011; Gonzalez-Vida et al., 2019)). HySEA models simulate both the granular mass and the water flow, as well as the coupling between the two *via* a friction term at the landslide/water interface (Macías et al., 2021a,b) (Figure 4). The landslide is described by the depth-averaged Savage-Hutter type equations. In 1-layer models, depth-averaged equations for the water motion are considered, including non-hydrostatic water pressures. In multilayer HySEA, on top of the landslide, several layers may be used to describe the water motion. This allows solving part of the vertical structure of the water flow, making it possible, for example, to describe velocity profiles (Fernández-Nieto et al., 2016; Garres-Díaz et al., 2020; Escalante et al., 2023). Adding vertical discretization (several layers) allows to better capture the dispersion effects as
shown in Macías et al. (2021a). In each layer, the shallow water equations are solved in the X-Y plane and the quantities are averaged in the vertical Z-direction over the layer thickness. Note that the shallow approximation and depth-averaging are made in the vertical direction for both the landslide and the water layers while it should be done in the direction perpendicular to the slope for the landslide and in the vertical direction for the water (Delgado-Sánchez et al., 2020). However, models accounting for both this effect and multilayer discretization do not presently exist.

## 3.1 Landslide-tsunami model with one water layer

Let us first recall the so-called 1-layer model which consists of 1 layer describing the landslide and 1 layer describing the water behavior. The equations include non-hydrostatic effects (i.e. dispersive effects) only in the water even though these terms also play a role in landslide dynamics (Garres-Díaz et al., 2021). The resulting equations read:

$$ 
\left\{
\begin{aligned}
\partial_t h + \partial_x (h u) + \partial_y (h v) &= 0 \\
\partial_t (h u) + \partial_x (h u^2 + \frac{1}{2} g h^2 + h p) + \partial_y (h u v) \\
&= (g h + p) \partial_x (H - h_s) - S_x \\
\partial_t (h v) + \partial_x (h u v) + \partial_y (h v^2 + \frac{1}{2} g h^2 + h p) \\
&= (g h + p) \partial_y (H - h_s) - S_y \\
\partial_t (h w) + \partial_x (h w u) + \partial_y (h w v) &= p \\
\nabla \cdot \mathbf{u} + \frac{w_\eta - w_0}{h} &= 0,
\end{aligned}
\right.
\quad (1) $$

where the unknowns of system (1) are the water surface height $h$, the horizontal water velocity vector $\mathbf{u} = (u \ v)$, the vertical water velocity $w$, the non-hydrostatic pressure $p$, the landslide thickness $h_s$ and its horizontal velocity $\mathbf{u}_s = (u_s \ v_s)$. In the sea, the total water height $h$ is defined as $h = H + \eta$ where $H(x, y)$ is the height of the water at rest above the bed and $\eta$ the sea surface elevation. Inland, $H(x, y)$ becomes negative and $-H > 0$ denotes the topography elevation above the sea level at rest. The water free surface elevation above the water level at rest becomes $h - H$ and $h$ is still the height of the water column above the topography. Figure 4 describes the notations, both inland and in the sea. Equation (1) simulates the water flow inland and offshore the same way.
No distinction is made by the model between inland and offshore zones and no additional model or treatment is needed to resolve the flow inland, outside the usual wet/dry transition treatment at small water depth (Macías et al., 2021a; Castro et al., 2005). The gravity acceleration is denoted by $g$. The last equation of system (1) is the incompressibility condition and involves the vertical velocity at the free surface $w_\eta$ and the vertical velocity at the bottom $w_0 = -\partial_t(H - h_s)$. The non-hydrostatic pressure is assumed to be zero at the free surface. The friction term

$$ \mathbf{S} = (S_x, S_y) = g h \frac{n^2}{h^{4/3}} \|\mathbf{u}\| \mathbf{u} - m_f (\mathbf{u}_s - \mathbf{u}) $$

involves the empirical Manning coefficient $n$, controlling the water/bottom friction, and the empirical granular mass/water friction coefficient $m_f$.

The landslide unknowns $h_s$, $u_s$, $v_s$ obey the reduced-gravity Savage-Hutter model (Savage & Hutter, 1989), which reads

$$ \left\{ \begin{aligned} & \partial_t h_s + \partial_x (h_s u_s) + \partial_y (h_s v_s) = 0 \\ & \partial_t (h_s u_s) \\ & \qquad + \partial_x (h_s u_s^2 + \frac{1}{2}(1-r)gh_s^2) + \partial_y (h_s u_s v_s) \\ & = (1-r)gh_s \partial_x (H-h_s) - S_{s,x} \end{aligned} \right. \qquad (2) \\ & \partial_t (h_s v_s) \\ & \qquad + \partial_x (h_s u_s v_s) + \partial_y (h_s v_s^2 + \frac{1}{2}(1-r)gh_s^2) \\ & = (1-r)gh_s \partial_y (H-h_s) - S_{s,y}, \end{aligned} \right. $$

where the friction term writes

$$ \mathbf{S}_s = (S_{s,x}, S_{s,y}) = g(1-r)\mu h_s \frac{\mathbf{u}_s}{\|\mathbf{u}_s\|} + r m_f (\mathbf{u}_s - \mathbf{u}). $$

In equations (2), buoyancy induces a factor $1-r$ that multiplies the gravity $g$ where $r$ is the density ratio between the constant fluid density $\rho_f$ and the constant grain phase density $\rho_s$

$$ r = \frac{\rho_f}{\rho_s}. \qquad (3) $$

The buoyancy effects are only taken into account when the landslide is underwater. In this paper, the fluid density is taken to be $\rho_f = 1000 \text{ kg/m}^3$ and the grain phase density is $\rho_s = \rho_g \Phi$ where
$\rho_g = 2500 \text{ kg/m}^3$ is the grain density and $\Phi$ the solid volume fraction. The solid volume fraction $\Phi$ is included in the sensitivity analysis and its values are discussed in Subsection 3.3. The friction coefficient $\mu$ plays a key role in the granular mass motion (Savage & Hutter, 1989; Greve & Hutter, 1993; Pouliquen & Forterre, 2002; Mangeney et al., 2007). Three different formula are tested in this paper, going from the constant Coulomb friction coefficient $\mu_C$ to the more complex Pouliquen and Forterre friction coefficient $\mu_P$ as well as the $\mu(I)$ rheology:

$$ \mu_C = \tan(\delta_1), \qquad (4) $$

$$ \mu_P = \begin{cases} \tan(\delta_1) + \frac{\tan(\delta_2) - \tan(\delta_1)}{1 + \beta \frac{h_s}{d Fr}} & \text{if } Fr \ge \beta \\ \mu_{start} + \left(\frac{Fr}{\beta}\right)^{\chi} (\mu_{stop} - \mu_{start}) & \text{if } Fr < \beta \end{cases}, \qquad (5) $$

$$ \mu(I) = \tan(\delta_1) + \frac{\tan(\delta_2) - \tan(\delta_1)}{1 + \frac{I_0}{I}}. \qquad (6) $$

The parameter $d$ is of the same order of magnitude as the grain diameter and is included in our sensitivity analysis (Table 1). The Froude number is defined as

$$ Fr = \begin{cases} \frac{\|\mathbf{u}_s\|}{\sqrt{g(1-r)h_s}} & \text{if } h_s > 0, \\ 0 & \text{otherwise.} \end{cases} \qquad (7) $$

We use $\beta = 0.136$ as in Pouliquen & Forterre (2002) even though this was calibrated for glass bead flows in the experiments. For $0 \le Fr < \beta$, the friction coefficient $\mu_P$ involves the following quantities

$$ \mu_{start} = \tan(\delta_3) + \frac{\tan(\delta_2) - \tan(\delta_1)}{1 + \frac{h_s}{d}} $$

$$ \mu_{stop} = \tan(\delta_1) + \frac{\tan(\delta_2) - \tan(\delta_1)}{1 + \frac{h_s}{d}} $$

and the empirical parameter $\chi = 0.001$ which is of small influence on the results as soon as $\chi < 0.01$ (Poulain et al., 2023). Equation (6) is a common expression resulting from steady flows on constant slopes (Jop et al., 2006) with the inertial number $I$ chosen here as in Martin et al. (2023)

$$ I = \frac{\frac{5}{2} \|\mathbf{u}_s\| d}{h_s \sqrt{g(1-r)h_s \Phi}}. \qquad (8) $$

The empirical value of the constant $I_0 = 0.279$ comes from Jop et al. (2006); Martin et al. (2023).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Multilayer structure and notations. The water depth $H$ is positive in the sea and is negative inland. Inland, $-H > 0$ is the topography.</figcaption></figure>

The solid volume fraction in the inertial number $I$ is empirically taken to be $\Phi = 0.6$ for all simulations. Changing $\Phi$ in the inertial number $I$ does not significantly affect the simulation results. Note that the $\mu(I)$ rheology make use of two angles $\delta_1$ and $\delta_2$, The Pouliquen and Forterre law $\mu_P$ involves three angles $\delta_1$, $\delta_2$ and $\delta_3$ while the Coulomb rheology only uses one friction angle $\delta_1$.

## 3.2 Multilayer discretization for the water motion

The landslide part is governed by the same equations since only the water column is divided into several layers. The water column of thickness $h$ is decomposed into $L$ sublayers of size $l_\alpha h$, where $\sum_{\alpha=1}^L l_\alpha = 1$ as illustrated in Figure 4. It is usually considered that $l_\alpha = 1/L$, although other definitions exist (Bonaventura et al., 2018). The unknowns of system (9) are now the total water thickness $h$, the horizontal water velocity vector $\mathbf{u}_\alpha = (u_\alpha v_\alpha)$, the vertical water velocity $w_\alpha$ for each layers $\alpha = 1, \dots, L$ and the non-hydrostatic pressure $p_{\alpha-1/2}$ at the lower interface of each layer

$$z = z_{\alpha-1/2} = \sum_{\beta=1}^{\alpha-1} l_{\beta}h + h_s - H.$$

The vertical averaged velocity is given by $(\bar{u}, \bar{v}) = \sum_{\alpha=1}^L l_\alpha \mathbf{u}_\alpha$.
The mass and momentum equations including non-hydrostatic (i.e. dispersive) effects read

$$ \left\{ \begin{array}{l} \partial_t U_f + \partial_x F_f + \partial_y G_f + B_f = H_f + K_{NH} - S_f, \\ \nabla_{NH} \cdot U_f = 0. \end{array} \right. \qquad (9) $$

The conservative terms $F_f$, $G_f$, and the pressure gradient $H_f$ are

$$ \begin{aligned} U_f &= \begin{bmatrix} h \\ \begin{bmatrix} h u_\alpha \\ h v_\alpha \\ h w_\alpha \end{bmatrix}_{\{\alpha=1,...,L\}} \end{bmatrix}, \\ F_f &= \begin{bmatrix} h \bar{u} \\ \begin{bmatrix} h u_\alpha^2 \\ h u_\alpha v_\alpha \\ h u_\alpha w_\alpha \end{bmatrix}_{\{\alpha=1,...,L\}} \end{bmatrix}, \\ G_f &= \begin{bmatrix} h \bar{v} \\ \begin{bmatrix} h u_\alpha v_\alpha \\ h v_\alpha^2 \\ h v_\alpha w_\alpha \end{bmatrix}_{\{\alpha=1,...,L\}} \end{bmatrix}, \\ H_f &= \begin{bmatrix} 0 \\ -g h \partial_x (h + h_s - H) \\ -g h \partial_y (h + h_s - H) \\ 0 \end{bmatrix}_{\{\alpha=1,...,L\}}. \end{aligned} $$

The breaking dissipation model and the friction terms are defined in $S_f$. The friction terms $S_f$ are only non-zero in the bottom layer $\alpha = 1$ and the breaking model is defined for the vertical
momentum equation at each layer

$$ S_f = \begin{bmatrix} 0 \\ g h \frac{n^2}{l_1 h^{4/3}} u_1 \sqrt{u_1^2 + v_1^2} - m_f(u_s - u_1) \\ g h \frac{n^2}{l_1 h^{4/3}} v_1 \sqrt{u_1^2 + v_1^2} - m_f(v_s - v_1) \\ m_1^w \\ 0 \\ 0 \\ m_\alpha^w \end{bmatrix}_{\{\alpha=2,...,L\}} $$

Following (Scala et al., 2024, Equation (8)), the coefficients for the breaking dissipation model are defined in the following form,

$$ m_{\alpha}^{w} = \frac{C}{l_{\alpha}^{2}} \sqrt{(\partial_{x}(hu_{\alpha}))^{2} + (\partial_{y}(hv_{\alpha}))^{2}} w_{\alpha}, \quad (10) $$

for $\alpha = 1, ..., N$, with $C = 35 \max(\frac{\sqrt{u^2+v^2}}{0.4\sqrt{gh}} - 1, 0)$ (Roeber et al., 2010).

The non-conservative products $B_f$ describe the momentum transfer across layer interfaces

$$ B_f = \begin{bmatrix} 0 \\ \left[ \frac{(u_{\alpha+1/2}\Gamma_{\alpha+1/2} - u_{\alpha-1/2}\Gamma_{\alpha-1/2})}{l_{\alpha}} \right] \\ \left[ \frac{(v_{\alpha+1/2}\Gamma_{\alpha+1/2} - v_{\alpha-1/2}\Gamma_{\alpha-1/2})}{l_{\alpha}} \right] \\ \left[ \frac{(w_{\alpha+1/2}\Gamma_{\alpha+1/2} - w_{\alpha-1/2}\Gamma_{\alpha-1/2})}{l_{\alpha}} \right] \end{bmatrix}_{\{\alpha=1,...,L\}} $$

with $\Gamma_{1/2} = \Gamma_{L+1/2} = 0$ and

$$ \Gamma_{\alpha+1/2} = \sum_{\beta=1}^{\alpha} l_{\beta} (\partial_x(h(\bar{u} - u_{\beta})) + \partial_y(h(\bar{v} - v_{\beta}))) \quad \alpha = 1, \dots, L-1. $$

Interface velocities are defined as $u_{\alpha+1/2} = \frac{1}{2}(u_{\alpha} + u_{\alpha+1})$ and analogous definitions for $v_{\alpha+1/2}$
and $w_{\alpha+1/2}$. The term $K_{NH}$ takes into account the non-hydrostatic effects

$$ K_{NH} = \left[ \begin{bmatrix} 0 \\ -\partial_x(h p_\alpha) + (p_{\alpha+1/2} \partial_x z_{\alpha+1/2} \\ -p_{\alpha-1/2} \partial_x z_{\alpha-1/2})/l_\alpha \\ -\partial_y(h p_\alpha) + (p_{\alpha+1/2} \partial_y z_{\alpha+1/2} \\ -p_{\alpha-1/2} \partial_y z_{\alpha-1/2})/l_\alpha \\ -(p_{\alpha+1/2} - p_{\alpha-1/2})/l_\alpha \end{bmatrix}_{\{\alpha\}} \right], $$

for $\alpha = 1, \dots, L$, where $p_\alpha = (p_{\alpha+1/2} + p_{\alpha-1/2})/2$ and

$$ z_{1/2} = h_s - H, \quad z_{L+1/2} = \eta = h + h_s - H $$

$$ p_{L+1/2} = 0. $$

Note that the pressure at the free surface is usually set to zero. Finally, the system is completed with the discrete incompressibility condition at each layer

$$ \nabla_{NH} \cdot \mathbf{U}_f = 0, $$

being

$$ h \nabla_{NH} \cdot \mathbf{U}_f = \left[ \begin{bmatrix} l_1 h \partial_x (h u_1) - l_1 h u_1 \partial_x h \\ + l_1 h \partial_y (h v_1) - l_1 h v_1 \partial_y h \\ + 2(h w_1 - h w_0) \\ (h u_{\alpha+1} - h u_{\alpha}) \partial_x z_{\alpha+1/2} \\ + (h v_{\alpha+1} - h v_{\alpha}) \partial_y z_{\alpha+1/2} \\ -(h w_{\alpha+1} - h w_{\alpha}) \end{bmatrix}_{\{\alpha\}} \right], $$

for $\alpha = 1, \dots, L-1$ and where $w_0$ is still defined by $w_0 = -\partial_t(H - h_s)$.

Multilayer HySEA implements a parallelized version of the described numerical method on GPU architectures (Castro et al., 2011; Escalante et al., 2018). The numerical discretization of the system is presented in Appendix F. The reader is referred to Fernández-Nieto et al. (2008); Macías et al. (2021a,b); Escalante et al. (2018) for more details about the numerical method.
3.3 Unknown parameters and model assumptions for the sensitivity tests

We perform a series of numerical simulations, listed in Table 1, by varying key parameters involved
in the model around a reference simulation (Section 4). In all simulations, the final time is set to
t<sub>f</sub> = 3600 s and open boundary conditions are prescribed thanks to sponge layers absorbing the
waves. The size of this sponge layer is defined to be 4 cells and thus depends on the grid resolution.
The reference simulation has the following features: (i) non-hydrostatic 3 layers simulation; (ii)
30 m bathymetry grid; (iii) μ(I) rheology; (iv) first friction angle δ<sub>1</sub> = 7°; (v) a characteristic
grain size diameter d = 1 m; (vi) a solid volume fraction Φ = 0.8; (vii) the spatialized Manning
coefficient of Figure 2 varying with the bed properties between 0.025 s.m<sup>-1/3</sup> and 0.2 s.m<sup>-1/3</sup>;
and (viii) a landslide/water interface friction m<sub>f</sub> = 2 × 10<sup>-5</sup> m.s<sup>-1</sup>. Each simulation differs from
the reference case by only one parameter for consistency and comparison purposes. The only
exception concerns the Hydro simulation which is hydrostatic with 1 layer in the water column.
For all the simulations, the same friction angles' differences are kept for all the sets: δ<sub>2</sub> − δ<sub>1</sub> = 10°
and δ<sub>3</sub> − δ<sub>1</sub> = 2° as done in Pouliquen & Forterre (2002); Poulain et al. (2023); Brunet et al.
(2017) (see also the discussion in Poulain et al. (2023) for the choice of the friction angle). In the
following, we refer to a set of friction angles only by the first one δ<sub>1</sub>. We varied the first friction
angle δ<sub>1</sub> by decreasing/increasing it by 4°, giving δ<sub>1</sub> = 3°, δ<sub>1</sub> = 7°, and δ<sub>1</sub> = 11°. Note that
δ = 11° approximately corresponds to the friction angle given by the law μ<sub>c</sub> = V<sup>-0.0774</sup>, with V
the landslide volume, fitted empirically by Lucas et al. (2014) on subaerial landslides. The effect
of the solid volume fraction Φ has been tested in the grain phase density. The order of magnitude
of the maximum packing fraction of mono-disperse granular media is 0.6 but this volume fraction
increases for poly-disperse materials (Andreotti et al., 2013). We have thus considered two typical
values Φ = 0.6 and Φ = 0.8 in this paper. We also considered different values for the parameter d
(0.1 m and 0.5 m) as well as for the Manning and interface friction parameters, which are generally
highly uncertain.

The 10 m resolution DEM from Lemoine et al. (2020b) is degraded for obtaining the 30 m
and 50 m grids. The resolution degradation is performed by removing rows and columns from
the raster of data to avoid additional treatments (like interpolations) and uncertainties. The 30 m
resolution for the grid is chosen for the reference case for its reasonable computation time given all the performed simulations and our computational capabilities: Appendix E Table A3 as well as Figure A10. To quantify bathymetry effect, we: (i) reduce the quality of the 30 m Digital Elevation Model (DEM) while keeping the same number of grid points in the numerical scheme; and (ii) decrease/increase the accuracy of the DEM with decreasing/increasing the number of grid points. We consider 10 m, 30 m and 50 m grids in this paper. Concerning the first point (i), we reduce the bathymetry quality by smoothing the DEM: 'Smoothed' case in Table 1. The smoothing process is implemented using a 3 × 3 kernel window Gaussian filter, which is commonly used for bathymetric data as it balances the need for noise reduction and the preservation of the spatial resolution. The kernel's coefficients were calculated based on the Gaussian function, which is defined as

$$G(x, y) = \frac{1}{2\pi\sigma^2}e^{-\frac{x^2+y^2}{2\sigma^2}}$$

where x and y are the distances from the kernel's center, and σ is the standard deviation of the Gaussian distribution. The value σ = 3 was chosen to ensure that the filter effectively smooths the data without overly blurring important features of the bathymetry.

Another source of uncertainty is the initial water level which varies with time due to the tide and which is expected to increase between 63 cm up to 1.32 m at the end of the century (Pörtner et al., 2019; Intergovernmental Panel on Climate Change (IPCC), 2023). We consider a case where the initial water level at rest is increased by 1 m: '+1 m water' case in Table 1.

Concerning the model assumptions, the hydrostatic approximation, leading to much simpler equations and comparatively very low computational cost, is not always valid. Non-hydrostatic pressures are sometimes required, as widely observed in the literature (Kirby et al., 2022; Poulain et al., 2023). Outside of the lagoon of Mayotte, it is likely to find wavelengths with the same order of magnitude as the water depth H (Poulain et al., 2022), thus challenging the hydrostatic approximation. We will investigate here the differences between hydrostatic and non-hydrostatic simulations.

Another key point here is to investigate the effect of adding layers in the vertical directions (i.e. increasing the vertical discretization) to describe the water motion in the non-hydrostatic case.
Three layers are recommended by Macías et al. (2021a) corresponding to the better compromise between accuracy and computational cost. Indeed, they show that a 3-layers model reproduce the phase velocity from the Airy's theory with a good accuracy for $0 \le k H \le 15$ where k is the wave number. Using less layers significantly increases the error in this range of k H. Here, the number of layers is taken between 1 and 10. Finally, the breaking wave mechanism described in Equation (10) has also been tested (Scala et al., 2024). We restrict ourselves to the parameters listed in Table 1 for the sake of clarity of this paper and these results are shown in the supplementary material. The maximum sea surface elevation is shown in Figure S10 (supplementary material) when the breaking model is active. For comparison purposes, Figure S10 also includes the maximum sea surface elevation for the reference case (Figure 9c). This feature affects the maximum sea surface elevation but the difference remains limited. In particular, the same amount of waves seems to pass the lagoon in both cases.

## 4 LANDSLIDE AND WATER WAVE DYNAMICS FOR THE REFERENCE CASE

Figure 5 shows the water wave propagation in the reference case. The landslide mass flows downwards the slope at a $54 \text{ m.s}^{-1}$ maximum velocity (Figure S3 in the supplementary material), almost stops at 250 s and reaches a runout distance of 8.14 km (Table A1). We clearly see the influence of the lagoon that dissipates and slows down the waves.

Water waves take around 60 s to reach the closest coast of Petite Terre, 170 s to reach the airport gauge n° 1, and 805 s to reach the DEALM gauge n° 2 on Grande Terre. Sea surface elevations of around 3 m height propagating at $2 \text{ m.s}^{-1}$ are recorded around 60 s after the event, when they first hit the eastern coast of Petite Terre (Figure 5a for the first impact location). This illustrates the threat represented by the scenario Piton 200. In that eastern part of Petite Terre where the depth $H \approx 2 \text{ m}$, the computed velocity of $2 \text{ m.s}^{-1}$ is of the same order of magnitude than the shallow velocity approximation $\sqrt{gH} \approx 4.5 \text{ m.s}^{-1}$. Maximum sea surface elevation and maximum surface velocities are shown later in the paper when discussing the other simulations (Figures 9 and Figure S1 in the supplementary material). The landslide dynamics and the maximum landslide velocity are also displayed in Figures S2 (top row) and S3 in the supplementary material. Figure
<figure><img src="image_6.png" /><figcaption>Figure 5. Propagation of water waves for the reference case: sea surface elevation at different times.</figcaption></figure>
<table><caption>**Table 1.** Performed simulations. Friction coefficients μ are respectively given by Equations (4), (5) and (6). The final time is set to $t_f = 3600$ s for all simulations. * the 30 m bathymetry is smoothed thanks to a Gaussian filter of kernel $3 \times 3$. Hydro type: Hydrostatic (H) or Non-hydrostatic (NH).</caption><thead><tr><th rowspan="2">Simulation</th><th rowspan="2">Hydro type</th><th rowspan="2">Layer n°</th><th rowspan="2">Grid res. (m)</th><th rowspan="2">Solid volum. fraction Φ</th><th rowspan="2">Grain size d (m)</th><th colspan="3">Friction angles</th><th rowspan="2">Rheology</th><th rowspan="2">Manning (s.m<sup>-1/3</sup>)</th><th rowspan="2">Mass water friction (m/s) × 10<sup>-5</sup></th><th rowspan="2">Initial water elevation (m)</th></tr><tr><th>δ (°)</th><th>δ<sub>1</sub></th><th>δ<sub>2</sub></th><th>δ<sub>3</sub></th></tr></thead><tbody><tr><td>Ref</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>Smoothed</td><td>NH</td><td>3</td><td>30*</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>Hydro</td><td>H</td><td>1</td><td>30</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>1 Layer</td><td>NH</td><td>1</td><td>30</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>6 Layers</td><td>NH</td><td>6</td><td>30</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>10 Layers</td><td>NH</td><td>10</td><td>30</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>10 m grid</td><td>NH</td><td>3</td><td>10</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>50 m grid</td><td>NH</td><td>3</td><td>50</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>Φ = 0.6</td><td>NH</td><td>3</td><td>10</td><td>0.6</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>0.1 m grain</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>0.1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>0.5 m grain</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>0.5</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>3° friction</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>1</td><td>3 13</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>11° friction</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>1</td><td>11 21</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>Coulomb</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>1</td><td>7</td><td>—</td><td>μ<sub>C</sub></td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>Pouliquen</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>1</td><td>7 17</td><td>9</td><td>μ<sub>P</sub></td><td>Fig. A2</td><td>2</td><td>0</td><td>0</td></tr><tr><td>Cst Manning</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>0.025</td><td>2</td><td>0</td><td>0</td></tr><tr><td>m<sub>f</sub> = 0</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>0</td><td>0</td><td>0</td></tr><tr><td>+1 m water</td><td>NH</td><td>3</td><td>30</td><td>0.8</td><td>1</td><td>7 17</td><td>—</td><td>μ(I)</td><td>Fig. A2</td><td>2</td><td>+1</td><td>0</td></tr></tbody></table>

6 represents the sea surface elevation at the gauge locations in the reference case (dark blue line). The sensitivity range of the sea surface elevation obtained from all the simulations listed in Table 1 displays the significant impact of the tested parameters. The sea surface elevation varies by a factor around 2-3 and up to 4 approximately compared to the reference case (see the first variation band peak on gauge n° 4).

Note that our results differ from what was obtained in Poulain et al. (2022) by coupling HySEA hydrostatic and FUNWAVE-TVD non-hydrostatic simulations. This may be due to the different approximation used (fully hydrostatic or fully non-hydrostatic here) and to the influence of the coupling time considered in Poulain et al. (2022). Poulain et al. (2022) simulate the landslide and generated near-field waves with hydrostatic HySEA and define a coupling time at which the com-
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Sea surface elevation η for the reference case and variation band given by all the simulations of Table 1.</figcaption></figure>

puted water state is transferred to FUNWAVE-TVD. At the coupling time, the computed waves in HySEA (hydrostatic simulation) are used to initialize FUNWAVE-TVD (non-hydrostatic simulation). These waves continue to propagate thanks to FUNWAVE-TVD.

## 5 KEY ROLE OF NON-HYDROSTATIC EFFECTS

Figures 7a and 7c show the maximum sea surface elevation $\zeta_{ij}$

$$ \zeta_{ij} = \max_{t} \eta_{ij}(t), \qquad (11) $$

for the hydrostatic (Hydro) and non-hydrostatic (1 Layer) cases, both dealing with 1 layer to describe the water motion. Inland on Petite Terre, the inundation pattern of Figures 7b and 7d correspond to the maximum inundation depth

$$ \max_{t} h_{ij}(t) \qquad (12) $$

Note that, inland, the sea surface elevation relative to the water at rest writes $\max_{t} h_{ij}(t) - H$ (Figure 4). Recall that multilayer HySEA resolves water flow inland the same way as offshore. No treatment is needed to create the inundation maps. The three resolutions considered in this paper do not allow for analyses of inundation patterns at building scale. The inundation maps are
<figure><img src="image_8.png" /><figcaption>**Figure 7.** a): maximum sea surface elevation $\ell_{ij}$ (Equation (11)) for the hydrostatic case. b): maximum inundation depth for the hydrostatic case. c): maximum sea surface elevation $\ell_{ij}$ (Equation (11)) for the non-hydrostatic case (1 layer). d): maximum inundation depth for the non-hydrostatic case (1 layer).</figcaption></figure>

intended to provide a general trend and a comparison point to determine the influence of input parameters.

Non-hydrostatic effects drastically change the result with major differences observed east of the lagoon. Overall, the hydrostatic simulation shows larger waves with smoother spatial variation of the maximum sea surface elevation. The waves overcome the reef and enter the lagoon more
a) Hydro (1 Layer)

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

b) Non-Hydro (1 Layer)

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

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

c) Spectral amplitude

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Spectrogram and spectral amplitude of the sea surface elevation signal at the gauge n° 4: Offshore in the hydrostatic (a-c) and non-hydrostatic 1-layer (b-c) cases.</figcaption></figure>

easily in the hydrostatic case. Around the Piton 200 collapse area, the wave pattern reflects the landslide trajectory very precisely only in the non-hydrostatic case. For both hydrostatic and non-hydrostatic cases, no more than 1.5 m waves are observed throughout the simulation inside the lagoon. The waves height passed the lagoon are mostly below 1 m; and even below 0.5 m in the larger part of the lagoon for the non-hydrostatic case. Waves are largely stopped by the coral reef which offers a great protection for the coasts of Mayotte. This artificial impact of the reef in the 1-layer models differs from what is obtained with the 3-layers reference case as discussed later (Section 6).

Waves however reach the east coast of Petite Terre but the high relief (+10 m) prevents the water from entering inland. We also see higher waves in the south east part of Petite Terre in the hydrostatic case. In Figures 7b and 7d, a larger part of the airport area is covered by water in the hydrostatic case (lower than 1.5 m). This is also illustrated at the airport gauge n° 1 in Figure S5 (supplementary material) where the maximum wave is almost 3 times higher in the hydrostatic case. Much higher maximum waves are obtained at the offshore gauge n° 4 in front
of the landslide, consistent with what is observed in laboratory experiments of granular flows and generated waves (Poulain et al., 2023). At this gauge n° 4, the waves are completely attenuated after 250 s in the hydrostatic case while the signal is still visible around 600 s in the non-hydrostatic case. Overall, the hydrostatic simulation overestimates the first waves (or the first water drop in the Piton 200 collapse area) up to a factor 3 compared to the non-hydrostatic 1-layer case. However, in the Mayotte configuration, the overestimation related to the hydrostatic approximation is not so much impactful on the coasts due to the artificial reef protection.

At early times, the wave signal contains much higher frequencies (>40 mHz) in the hydrostatic case as illustrated on the spectrograms at the offshore gauge n° 4 in Figure 8. This is due to the well-known sharp variations of the initial wave near the source in hydrostatic simulations (Poulain et al., 2023). The hydrostatic assumption strongly impacts the whole frequency content of the wave trains (Figure 8). This assumption is indeed no longer valid at the offshore gauge n° 4 for example. The shallow approximation assumes $H/λ << 1$, where λ is the tsunami wavelength. If we assume a velocity scaling with $\sqrt{gH}$, this leads to

$$f << \sqrt{\frac{g}{H}}$$

At the offshore gauge n° 4, the water depth is $H \approx 1100$ m, requiring $f << 95$ mHz while the dominant frequency $f \approx 20$ mHz is only 5 times smaller (Figure 8). On the opposite, the depth at the airport gauge n° 1 is $H \approx 3$ m meaning $f << 1800$ mHz. We observe a dominant frequency of around 1 mHz in Figure S9 (supplementary material) for all the presented simulations at the airport gauge, which means that the hydrostatic assumption is valid at this gauge. Note that for frequencies smaller than $\approx 10$ mHz (periods $T \ge 100$ s), the frequency content of the hydrostatic and non-hydrostatic cases well match (Figure 8c). Contrary to the early times, when $t \ge 100$ s, the frequency content in the non-hydrostatic case is higher than in the hydrostatic case, describing more precisely the wave train following the first waves (Figure 8a-b).
<figure><img src="image_13.png" /><figcaption>**Figure 9.** a): maximum sea surface elevation $t_{ij}$ (Equation (11)) for 1 layer. b): maximum inundation depth for 1 layer. c): maximum sea surface elevation $t_{ij}$ (Equation (11)) for 3 layers (reference case). d): maximum inundation depth for 3 layers (reference case). e): maximum sea surface elevation $t_{ij}$ (Equation (11)) for 6 layers. f): maximum inundation depth for 6 layers.</figcaption></figure>

# 6 CRUCIAL NEED OF MULTILAYER MODELS AND HIGH RESOLUTION BATHYMETRY

## 6.1 Effect of vertical discretization (number of layers)

One of the key result of this work is the crucial role of vertical discretization in capturing strong topography variations as those related to the presence of the coral reef.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Sea surface elevation $\eta$ for different number of layers superposed with the reference case.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Spectrograms and spectral amplitude of the sea surface elevation signal at the gauge n° 4: Offshore.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Vertical profiles at the airport gauge n° 1 of (a) the velocity in the X-direction and (b) the velocity in the Z-direction. The selected time for each simulation corresponds to $t_a + 75$ s, where $t_a = 165$ s is the arrival time of the first wave at the gauge. The vertical velocity profile is thus plotted at $t = 240$ s.</figcaption></figure>

### 6.1.1 Maximum sea surface elevation

Figure 9 shows the significant differences in maximum sea surface elevation $t_{ij}$ between simulations as the number of layers increases. The barrier power of the coral reef decreases with more layers. The main difference is observed between 1 layer and 3 layers. Indeed, with 3 layers, wave heights reach 0.5-1 m in the lagoon while they are mostly smaller than 0.5 m with 1 layer (Figures 9a and 9b). Using 6 and 10 layers still increases the wave capacity to overcome the coral reef barrier (comparing Figures 9b and 9c). As a result, the 1-layer case dramatically underestimates the waves past the lagoon: we observe a factor 2 comparing with 3 layers (Figures 9a and 9c) and 6 layers (Figures 9a and 9e). Signed differences lower than -100% are found in the Piton 200 area when comparing 1- and 3-layers simulations (Figure A3a). The 1-layer simulation also shows more than 100% of difference with the 3-layers reference case inside the lagoon. The differences become much smaller as we continue to increase the number of layers. For example, we still have differences up to 60% between the 3- and 6-layer cases. For a reasonable accuracy and computational cost, the optimal number of layers seems to be slightly higher than 3 layers, as recommended by Macías et al. (2021a), probably due to strong bathymetry variations around Mayotte. The 6- and 10-layers simulations show similar results and tend to converge on all gauges (Figure 10). Interestingly, the frequency content for a number of layers greater than 3 is very similar with frequencies up to 45 mHz. On the contrary, the 1-layer model only captures frequencies $f < 35$ mHz ($T > 28$ s) as shown in Figure 11.
Concerning the wave impact on the Mayotte coasts, Figures 9b, 9d and 9f presents the inundation pattern on Petite Terre depending on the number of layers. In the multilayer cases (above 3 layers), the water enters and propagates inside Petite Terre over approximately 100 m on a 2 km coastline south west at the airport location. In that area, recorded maximum water height are around 1 m but 50 cm of water covers the most part of the inundated zone. The path leading to the most western part is largely flooded by 40 cm and a small area north west is also covered by water. The water height reaches approximately 80 cm in that area. The 1-layer model leads to much less inundation in Figure 9b: the water penetrates the land over 50 m in average and observed water height are around 30 cm. The high relief (+10 m) protects the east coast of Petite Terre from inundation. The sea surface velocity patterns shown in Figure S1 (supplementary material) looks roughly similar between each other. As for the maximum sea surface elevation, the major differences are observed for the 1-layer case: in the landslide area, where the final deposit pattern is more visible, and inside the lagoon, where the reef offers a higher protection. The 3-layer model leads to higher velocities. For example, inland at the airport location, water penetrates with velocities around $0.5 \text{ m.s}^{-1}$ in the 1-layer case (Figure S1b in the supplementary material) and around $1.5 \text{ m.s}^{-1}$ in the 3-layer case (Figure S1d in the supplementary material). The overestimation of the barrier power of the coral reef with the 1-layer model is observed at the gauge n° 2 (Figure 10). More gauge results are given in Appendix D.

The increased accuracy with more water layers comes with the price of a higher computational cost as it is shown in Figure A10a. The real computational time is expected to grow linearly with the number of layers: a 2-layer simulation doubles the computational time compared to a 1-layer simulation. We observe indeed a linear trend in Figure A10a. The mean error of the simulated maximum sea surface elevation for 1 layer compared to the 10-layer simulation is non-negligible (Figure A10b). It stays below 20 cm as soon as we consider 2 layers or more, and it decreases rapidly as we increase the number of layers. We estimate a reasonable number of layers in term of accuracy and computational cost being 4 layers according to Figure A10b. We notice that the Mayotte configuration demands a number of layers higher than 3, while the gain from 6 layers or more seems low. In particular, the mean error decrease is non-negligible when going from 3 to
4 layers. Three layers is the recommended value from Macías et al. (2021a), based on analytical formula derived from the model to determine its dispersive properties.

### 6.1.2 Velocity profiles and power to overcome large reliefs

The capacity of the waves simulated by the multilayer model to overcome large reliefs depends on how the model captures vertical profiles of the velocity. Such information is intrinsically missing in depth-averaged 1-layer models. An example of velocity profiles is given in Figure 12 at the airport gauge n° 1. In multilayer models, even if the water flow is blocked by the bathymetry obstacle at depth, resulting in a strong drop of the horizontal velocity close to the sea floor, the surface velocity may stay large, allowing the water to overcome large reliefs. On the contrary, with only 1 layer, the bathymetric variations greatly impact the whole water column of constant velocity. The depth-averaged X-velocity in the 1-layer model is 1.5 times smaller than the averaged velocity obtained with 10 layers (if we remove the bottom layer) at the airport gauge n° 1 (Figure 12a). The bottom boundary layer sees a lower velocity due to friction while the rest of upper layers in the water column experiences a much higher velocity. While using 10 layers still increases the precision, the profiles of the horizontal and vertical velocities are close for 6 and 10 layers. It seems that we observe a convergence in the vertical profiles when increasing the number of layers. The velocity field is complex across the reef as shown in the movies *animation1Layer.mp4* and *animation6Layer.mp4* representing the velocity vectors along a 2D X-Z cross section between airport (left) and east airport (right) gauges (supplementary material). We may even observe higher horizontal velocities at the bottom layer than at the surface: see the animation *animationEastAirportVelX.mp4* between 120 s and 140 s, showing the time change of the velocity vertical profiles at the east airport gauge.

## 6.2 Bathymetry resolution and accuracy

The resolution of bathymetry data is expected to impact the simulation results even though there are not many in-depth studies of the resolution effects in literature. Figure 13c and 13e show similar patterns on the maximum sea surface elevation outside the lagoon for 30 m and 10 m bathymetry
<figure><img src="image_17.png" /><figcaption>**Figure 13.** a): maximum sea surface elevation $t_{ij}$ (Equation (11)) for a smoothed bathymetry on a 30 m grid. b): maximum inundation depth for a smoothed bathymetry on a 30 m grid. c): maximum sea surface elevation $t_{ij}$ (Equation (11)) for a 30 m grid (reference case). d): maximum inundation depth for a 30 m grid (reference case). e): maximum sea surface elevation $t_{ij}$ (Equation (11)) for a 10 m grid. f): maximum inundation depth for a 10 m grid. The dashed black line in a) and b) is the coastline calculated from the smoothed bathymetry while the full black line is the regular Histolitt coastline.</figcaption></figure>

grids. The capacity of waves to overcome the reef barrier is however higher on the 10 m grid. The east coast of Grande Terre mostly sees waves of around 1 to 1.5 m with the 10 m grid, as opposed to 0.5 to 1 m with the 30 m grid. In the west coast of Petite Terre, the 10 m grid noticeably lead to higher maximum inundation depth and higher velocities as illustrated in Figures 13d,f and Figure S1 (supplementary material). Inland at the airport location, water penetrates with high velocities
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Sea surface elevation $\eta$ for the different bathymetry superposed with the reference case.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Spectrograms and spectral amplitude of the sea surface elevation signal at the gauge n° 4: Offshore.</figcaption></figure>
around $1.5 \text{ m.s}^{-1}$ on the 30 m grid (Figure S1d in the supplementary material) and even $2.5 \text{ m.s}^{-1}$ on the 10 m grid (Figure S1f in the supplementary material). The high impact of the resolution of bathymetric data is also visible in Figure A4 showing differences in the sea surface elevation between the different grid simulations up to 60% – 80% in the lagoon.

We argue that these differences are related to the accuracy of the bathymetry and not to numerical dissipation. Indeed, the numerical dissipation is known to decrease with finer resolution meshes but is expected to be significant only after a few hours of simulations as recently shown in Svennevig et al. (2024). It would be interesting to run simulations on finer grids to more precisely quantify bathymetry impact on generated waves, but it is beyond the scope of this paper. Furthermore, it is quite rare to get bathymetry data in such a wide area with precision lower than 10 m. Finally, we check the effect of a smoothed, low resolution bathymetry described by a 30 m grid. Note that the dashed black coastline calculated from the smoothed bathymetry does not exactly fit with the Hystolitt coastline because of the smoothing process. While we observe roughly the same pattern on the maximum sea surface elevation between the smoothed and reference case (Figures 13a and 13c), the waves overcome the reef more easily with the smoothed grid. The differences on the simulated inundation areas is significant when looking locally around Petite Terre (Figures 13b and 13d). The inundation depth is much higher in the airport area due to lower topographic obstacles (Figure 13b): up to 3 m of water is observed. Taking into account the new coastline in dashed black defined by the smooth bathymetry, the east coast is almost totally under a few meters of water with penetration length going from 50 m to 100 m. As a result in our case, a low quality bathymetry may overestimate by up to 2 m the maximum water height inland and by 500 m the penetration length (airport location). The high relief located on the east coast of Petite Terre drops from 10 m to 3 m approximately due to the smoothing process.

Notably, the effect of the bathymetry strongly depends on the location and time. For example, the maximum positive wave observed at $t \sim 30 \text{ s}$ is much higher with the 10 m grid at the Piton 200 gauge $^\circ 3$ whereas the first simulated waves are the same for all the resolutions (Figure 14). On the contrary, the effect of bathymetry resolution is quite small at the offshore gauge $^\circ 4$,
but is significant at the DEALM gauge n° 2. The wave form is much smoother on the smoothed bathymetry on gauges n° 1 and 2.

A higher grid resolution allows to capture higher frequencies as shown by the spectrogram and spectral amplitude at the offshore gauge n° 4 (Figure 15). The typical wavelength associated with the wave train is $\lambda = \sqrt{gH}/f$, which gives $\lambda \approx 5$ km for $H \approx 1000$ m and $f \approx 20$ mHz. A numerical grid resolution of a few tens of meters is thus enough to handle such large wavelengths. We conclude that the higher frequency content with increasing resolution is related to the source (landslide) dynamics. The resolution of heterogeneities in the landslide dynamic thanks to higher grid resolution are indeed at a much smaller scale than 5 km (Figure S2 in the supplementary material). However, the grid resolution and bathymetry quality seem to weakly affect the contours of the flowing mass at 50, 100 and 3600 s. Differences appear on the flow thicknesses, especially at the early times 50 and 100 s.

Interestingly, a very long period oscillation $T \approx 600$ s is observed at the airport gauge n° 1 (spectral amplitude of Figure S9 in the supplementary material). This period corresponds to the typical period of a seiche that would go back and forth within the channel between north west of Petite Terre and Grande Terre, west of the airport gauge. Analytical formulas give a seiche period of $T_s = 2L/\sqrt{gH}$ in a closed rectangular basin of depth $H$ (Magdalena et al., 2020). In our case, the channel width is $L \approx 4$ km and the water depth $H \approx 15$ m, leading to $T_s \approx 650$ s. Only the multilayer simulations capture this frequency well.

# 7 LANDSLIDE RHEOLOGY, FRICTION PARAMETERS AND TIDE EFFECTS

## 7.1 Impact of landslide dynamics on the water waves

Despite lots of work, uncertainties still remain on the physical understanding of submarine landslide dynamics and associated parameters (Delannay et al., 2017; Roger et al., 2024). We will identify here which of these parameters mostly influence the simulated tsunami wave. The tsunami waveform is known to depend on the thickness distribution and velocity of the landslide, which in turn depend on the landslide rheology and associated friction coefficients. Figures 20 and S4 (supplementary material) illustrate that: the mass velocity and runout distance obviously increase when
<figure><img src="image_20.png" /><figcaption>**Figure 16.** a): maximum sea surface elevation $\ell_{ij}$ (Equation (11)) for a 3° friction angle. b): maximum inundation depth for a 3° friction angle. c): maximum sea surface elevation $\ell_{ij}$ (Equation (11)) for a 11° friction angle. d): maximum inundation depth for a 11° friction angle.</figcaption></figure>

lowering the friction coefficient $\mu(I)$, i.e. with decreasing friction angle $\delta_1$ and grain diameter $d$, as shown by Equation (6). The Pouliquen friction law (Equation (5)) gives a similar final deposit compared to the $\mu(I)$ rheology. The same conclusion holds for the solid volume fraction $\Phi = 0.6$ that gives a similar final deposit compared to the reference case with $\Phi = 0.8$. For all simulations, the front of the deposit contains the largest amount of mass with thicknesses ranging from 15 m up
<figure><img src="image_21.png" /><figcaption>**Figure 17.** a): maximum sea surface elevation $t_{ij}$ (Equation (11)) for a solid volume fraction $\Phi = 0.8$. b): maximum inundation depth for a solid volume fraction $\Phi = 0.8$. c): maximum sea surface elevation $t_{ij}$ (Equation (11)) for a solid volume fraction $\Phi = 0.6$. d): maximum inundation depth for a solid volume fraction $\Phi = 0.6$.</figcaption></figure>

to 40 m. The final runout, area, maximum thickness are collected in Table A1. The case $\delta_1 = 3^\circ$ leads to noticeably higher thicknesses on the front and therefore less mass left behind.

The maximum sea surface elevation (Figures A7 and A8) and the waveform at the gauge (Figure S6 and S7 in the supplementary material) only slightly depend on the grain diameter and
<figure><img src="image_22.png" /><figcaption>**Figure 18.** Sea surface elevation $\eta$ for the different sets of friction angles superposed with the reference case.</figcaption></figure>

on the type of friction law. However, the associated landslide deposit are quite different in Figures 20. On the contrary, the maximum sea surface elevation significantly depends on the values of the friction coefficients that have been tested: $\delta_1 = 3^\circ$ and $\delta_1 = 11^\circ$ (Figure 16). The maximum landslide velocities from Figure S3 (supplementary material) significantly differ between the $\delta_1 =$

<figure><img src="image_23.png" /><figcaption>**Figure 19.** Sea surface elevation $\eta$ for a solid volume fraction $\Phi = 0.6$ superposed with the reference case.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 20.** Final mass deposit for the different landslide parameters. White dashed lines indicate the final deposit for the reference case.</figcaption></figure>

3° and $\delta_1 = 11°$ cases while they appear to be closer when changing the grain diameter $d$ or the rheological law. Biggest differences in the maximum sea surface elevation are clearly found when comparing $\delta_1 = 3°$ and $\delta_1 = 11°$ simulations (Figure A5). The solid volume fraction value also has
<figure><img src="image_25.png" /><figcaption>**Figure 21.** Maximum inundation depth: a) +0 m water (reference case); b) +1 m water.</figcaption></figure>

a high impact on the maximum sea surface elevation (Figure 17), despite a similar final deposit for $Φ = 0.6$ compared to the reference case (Figure 20). We notice that increasing the solid volume fraction reduces the density ratio $r$ which causes a higher gravity force $(1-r)g$ for the landslide in the model. This is consistent with the observed maximum sea surface elevation in Figure 17 where the solid volume fraction $Φ = 0.6$ leads to higher waves, more concentrated around the landslide source, compared to the case $Φ = 0.8$. The high impact of the solid volume fraction $Φ$ suggests that using more complex and realistic rheology including dilatancy $(μ(I), Φ(I))$ would significantly affect the results (Iverson & George, 2014; Bouchut et al., 2016, 2025). The higher impact of both the friction angles and the solid volume fraction on the simulated waves is attributed to the early landslide dynamics. Figure A6 shows that the landslide thickness (resp. landslide velocity) from the simulations with $δ_1 = 3°, δ_1 = 11°$ and $Φ = 0.6$ differ from the rest of the simulations in the first 50 s (resp. 20 s). The impact of the friction angle and solid volume fraction is also observed on the gauge data in Figures 18 and 19. More gauge data related to landslide parameters are given in Appendix D and the supplementary material. The inundation pattern looks similar for all the landslide parameters except when comparing $δ_1 = 3°$ and $δ_1 = 11°$ (Figures 16, A7, A8b,d) and when comparing $Φ = 0.6$ and $Φ = 0.8$ (Figure 17b,d). In particular, for $δ_1 = 11°$, tsunami waves inundate less of the airport zone, which sees a maximum water height of around 0.4 m as opposed
to 0.7 m for $\delta_1 = 3^\circ$. The effect of the solid volume fraction is non-negligible on the inundation pattern as we observe less inundation for $\Phi = 0.6$.

## 7.2 Landslide/water and water/bottom friction analysis

We study here the effect of the Manning coefficient by comparing to a simulation with a constant manning coefficient ('Cst Manning' case). We also set the mass/water interface friction to 0 ('$m_f = 0$' case) and compare with the reference case ($m_f = 2 \times 10^{-5} \text{ m.s}^{-1}$). Values are specified in Table 1. The mass/water coefficient have a negligible effect on the water waves (Figure S4c in the supplementary material). The constant Manning coefficient leads to high differences compared to the spatialized Manning (reference case) on the generated waves: up to 90% differences are observed. With no surprise, these differences are located inside the lagoon (Figure S4a in the supplementary material) since a higher Manning value is applied at the reef location for the reference (Figure A2).

## 7.3 Tide effect

When adding 1 m to the water level, we observe that the waveform is similar but the water wave arrive slightly earlier (30 s) at the DEALM gauge n° 2 (Figure S8 in the supplementary material). This is consistent with a higher wave velocity $\sqrt{gH}$ when the water depth $H$ increases by 1 m. Note that on these figures, the sea surface elevation has been shifted (by -1 m) to match the hydrostatic zero of other simulations for comparison purposes. Unsurprisingly, the west coast of Petite Terre is largely impacted by a 1 m elevation of the initial water level: a larger part of west Petite Terre is under water and the water depth inland reaches 1 to 1.5 m in Figure 21b, especially south west around the aipport. The reference (+0 m water) only gives around 0.5 m water depth inland.

# 8 CONCLUSION

In this work, we perform a sensitivity analysis using the multilayer code HySEA on a island-type bathymetry characterized by large relief variations. We quantify the impact of the different model
parameters on the generated waves. Mayotte island is a relevant example due to the presence of the
submarine shelf surrounding Petite and Grande Terre. In the studied cases, the tsunami waves are
induced by a potential landslide of 200 millions of m³ initiated 2.5 km offshore Mayotte island.

The main result of this paper is the strong impact of the vertical discretisation (through mul-
tiple layers in the water column) when the bathymetry presents obstacles and steep slopes. In the
Mayotte configuration, the coral reef presents steep slopes with typical water depths ranging from
30 m (west of the reef) to 1000 m (east of the reef). The 1-layer (i.e. depth-averaged) model is
much more sensitive to this bathymetric barrier. The 1-layer model thus predicts waves that are
more attenuated inside the lagoon compared to multilayer models. On the west coast of Petite Terre
including the crucial zone occupied by the airport, the 1-layer model underestimates the impact
of the tsunami by a factor 2 compared to the 3-layer model. Adding more layers still increases
the capacity of the tsunami waves to overcome the coral reef, even though the model roughly
converges for 6 and 10 layers. The influence of increasing the number of layers is also visible in
the landslide area. Above the landslide, the wave pattern simulated with the 1-layer model clearly
reflects the landslide contours while multilayer models yield to a much smoother pattern, owing
to the higher degree of freedom for the water motion in the vertical direction. We show that the
numerical cost increases approximately linearly with the number of layers. In the particular case
of Mayotte, involving strong topography variations, a good balance between accuracy and compu-
tational cost seems to be around 4 layers. We expect that this could be the case for other locations
where bathymetry presents large variations and obstacles. Adding more layers allow to better cap-
ture the dispersion effects as detailed in Macías et al. (2021a), who recommend using 3 layers
in their simulations over much simpler topographies. Our results illustrate the major importance
of the dispersive effects (non-hydrostatic pressure and multilayer description) and the need to in-
clude them in the model, despite the increased computational cost, as recommended in Kirby et al.
(2022). Interestingly, our multilayer simulations capture the seiche of period T~s~ ≈ 600 s observed
at the airport gauge n°1, seemingly related to the water waves going back and forth between the
coasts of Grande and Petite Terre as observed and simulated in the Greenland fjords (Svennevig
et al., 2024).
Our results show that for the finest 10 m grid used here, the water waves overcome the reef barrier more easily than on rougher grids (30 m grids). A smoothing of the bathymetry (representing low quality bathymetric data) leads to a complete different result at the airport location for example: the airport area is completely submerged by a few meters of water as opposed to all other simulations. This points out the crucial need for more oceanographic campaigns to acquire accurate and high resolution bathymetric data in areas at risk. The major influence of the bathymetric data on landslide-tsunami simulations has also been highlighted by Svennevig et al. (2024).

The avalanche parameters also impact the amplitude and waveform of the tsunami, even locally, at Petite Terre for example. They obviously also change the landslide final deposit and mass distribution. While the form of the different friction laws investigated here and the grain diameter give quite similar water waves, the friction angles involved in these laws and the solid volume fraction have a large effect on both the landslide dynamics and tsunami waves. The impact on Petite Terre coastline is clearly visible: simulated wave amplitude doubles at the airport for $\delta_1 = 3°$ compared to $\delta_1 = 11°$. The landslide dynamic indeed affects the water waves generation and propagation.

While the mass/water friction coefficient is of small influence on the water waves, significant differences are observed when comparing the simulations with a constant and a variable Manning coefficient. Finally, a +1 m initial water level, potentially due to climate change, induces higher velocities, wider inundation distance, higher inundation depth and an earlier arrival time of around 30 s at gauge n°2 compared to the reference simulation.

Overall, the performed simulations allow us to roughly classify the parameters' influence on the maximum sea surface elevation and inundation pattern in the case of Mayotte from the more influential parameters to the least: 1) dispersion properties (non-hydrostatic versus hydrostatic and number of layers); 2) quality of bathymetric data, grid resolution, friction angles $\delta_s$ in the rheological laws and the solid volume fraction $\Phi$; 3) initial water level; 4) Manning coefficient; 5) type of rheological laws ($\mu(I)$ versus Coulomb); 6) mass/water friction and grain size related parameter $d$. The same classification applies locally at the Ponton gauge n°5 (Figure A9).

Those considerations should be useful to develop modeling strategy in case of near field con-
figuration for landslides induced tsunamis in a complex bathymetry context. Future works may include a deeper insight into the influence of the parameters or features targeted by this study using statistical analysis. As pointed out in this paper, a proper convergence with even finer grids and/or an increasing number of layers would be interesting to carry out. Additional features related to a multilayer model, like vertical viscosity terms and breaking wave mechanisms, would also be worth investigating. This sensitivity analysis shows that the parameter uncertainties and model hypotheses may lead to a factor of approximately 2 on the simulated wave amplitude and velocity. Consequently, we observe about the same error on the penetration length and the inundation depth. Increasing the level of accuracy of landslide-tsunami models would have a significant impact on hazard assessment by providing actionable insights for evacuation planning, such as defining more precise safe zones and optimizing warning times.

ACKNOWLEDGMENTS

A. M., A. M. and F. L. are part of the DT-GEO project funded by the European Union's Horizon Europe research and innovation programme under grant agreement n°101058129. M. J.C. has been partially supported and funded by MCIN/AEI/10.13039/501100011033, by the "European Union NextGenerationEU/PRTR" through the grant PDC2022-133663-C21, by MCIN/AEI/10.13039/501100011033, by "ERDF A way of making Europe", by the "European Union" through the grant PID2022-137637NB-C21. E. F.-N. has been partially supported by the European Union - NextGeneration EU program, by grant PID2022-137637NB-C22 funded by MCIN/AEI/10.13039/501100011033, the project "ERDF A way of making Europe", and by Junta de Andalucía research project ProyEx-cel_00525. A. L.F. benefits from funds managed by the Agence Nationale de la Recherche under the France 2030 programme with the reference ANR-22-EXIR-0005, (PEPR Risques, IRiMa). A. L. acknowledges the support from the Programme National de Télédétection Spatiale (PNTS, Grant PNTS-2022) and the LabEx UnivEarthS (ANR-10-LABX-0023 and ANR-18-IDEX-0001). R.P. and A.L. benefits from fundings of the french Ministry for the Ecological Transition and BRGM.
DATA AVAILABILITY

The simulations are performed with the code multilayer HySEA developed by the EDANYA group of Malaga University in Spain. The continuous equations and the associated numerical discretization are described in this paper. The reader is referred to the references cited in this paper and to the website available at https://edanya.uma.es/hysea/ for more details about the HySEA numerical codes. The input parameters used for the simulations are all mentioned in the paper. Bathymetric data are sourced and they may be distributed on demand. The maps are produced with the QGis open-source software (https://qgis.org) and homemade Python scripts for gauge data. The authors are available to provide any information needed to reproduce the results.

References

Abadie, S., Morichon, D., Grilli, S., & Glockner, S., 2010. Numerical simulation of waves generated by landslides using a multiple-fluid navier–stokes model, *Coastal Engineering*, **57**, 779–794, doi: 10.1016/J.COASTALENG.2010.03.003.

Allgeyer, S., Bristeau, M. O., Froger, D., Hamouda, R., Jauzein, V., Mangeney, A., Sainte-Marie, J., Souillé, F., & Vallée, M., 2019. Numerical approximation of the 3d hydrostatic navier–stokes system with free surface, *ESAIM: Mathematical Modelling and Numerical Analysis*, **53**, 1981–2024, doi: 10.1051/M2AN/2019044.

Andreotti, B., Pouliquen, O., & Forterre, Y., 2013. *Granular Media: Between Fluid and Solid*, Cambridge University Press, Cambridge.

Audusse, E., Bouchut, F., Bristeau, M. O., Klein, R., & Perthame, B., 2006. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows, https://doi.org/10.1137/S1064827503431090, **25**, 2050–2065, doi: 10.1137/S1064827503431090.

Aïssiouene, N., Bristeau, M. O., Godlewski, E., Mangeney, A., Madroñal, C. P., & Sainte-Marie, J., 2020. A two-dimensional method for a family of dispersive shallow water models, *SMAI Journal of Computational Mathematics*, **6**, 187–226, doi: 10.5802/SMAI-JCM.66/.

Behrens, J. & Dias, F., 2015. New computational methods in tsunami science, *Philosophical*
*Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **373**,
doi: 10.1098/RSTA.2014.0382.

Bonaventura, L., Fernández-Nieto, E. D., Garres-Díaz, J., & Narbona-Reina, G., 2018. Multilayer shallow water models with locally variable number of layers and semi-implicit time discretization, *Journal of Computational Physics*, **364**, 209–234, doi: 10.1016/J.JCP.2018.03.017.

Bonilauri, E. M., Aaron, C., Cerminara, M., Paris, R., Ongaro, T. E., Calusi, B., Mangione, D., & Harris, A. J. L., 2024. Inundation and evacuation of shoreline populations during landslide-triggered tsunami: An integrated numerical and statistical hazard assessment, doi: 10.5194/egusphere-2024-221.

Bouchut, F., Fernández-Nieto, E. D., Mangeney, A., & Narbona-Reina, G., 2016. A two-phase two-layer model for fluidized granular flows with dilatancy effects, *Journal of Fluid Mechanics*, **801**, 166–221, doi: 10.1017/jfm.2016.417.

Bouchut, F., Drach, E., Fernández-Nieto, E. D., Mangeney, A., & Narbona-Reina, G., 2025. A series of two-phase models for grain-fluid flows with dilatancy, *Journal of Fluid Mechanics*, **1008**, A43, doi: 10.1017/jfm.2025.131.

Brunet, M., Moretti, L., Friant, A. L., Mangeney, A., Nieto, E. D. F., & Bouchut, F., 2017. Numerical simulation of the 30–45 ka debris avalanche flow of montagne pelée volcano, martinique: from volcano flank collapse to submarine emplacement, *Natural Hazards*, **87**, 1189–1222, doi: 10.1007/S11069-017-2815-5.

Bunya, S., Dietrich, J. C., Westerink, J. J., Ebersole, B. A., Smith, J. M., Atkinson, J. H., Jensen, R., Resio, D. T., Luettich, R. A., Dawson, C., Cardone, V. J., Cox, A. T., Powell, M. D., Westerink, H. J., & Roberts, H. J., 2010. A high-resolution coupled riverine flow, tide, wind, wind wave, and storm surge model for southern louisiana and mississippi. part i: Model development and validation, *Monthly Weather Review*, **138**, 345–377, doi: 10.1175/2009MWR2906.1.

Cannat, M., Mangeney, A., Ondréas, H., Fouquet, Y., & Normand, A., 2013. High-resolution bathymetry reveals contrasting landslide activity shaping the walls of the mid-atlantic ridge axial valley, *Geochemistry, Geophysics, Geosystems*, **14**, 996–1011, doi: 10.1002/GGGE.20056.

Castro, M. J. & Parés, C., 2020. Well-balanced high-order methods for systems of balance laws,
*Continuum Mechanics, Applied Mathematics and Scientific Computing: Godunov's Legacy: A Liber Amicorum to Professor Godunov*, pp. 69–75, doi: 10.1007/978-3-030-38870-6_10.

Castro, M. J., Ferreiro, A. M. F., García-Rodríguez, J. A., González-Vida, J. M., MacíAs, J., Parés, C., & Vázquez-Cendón, M. E., 2005. The numerical treatment of wet/dry fronts in shallow flows: application to one-layer and two-layer systems, *Mathematical and Computer Modelling*, **42**, 419–439, doi: 10.1016/j.mcm.2004.01.016.

Castro, M. J., Ortega, S., de la Asunción, M., Mantas, J. M., & Gallardo, J. M., 2011. Gpu computing for shallow water flow simulation based on finite volume schemes, *Comptes Rendus Mécanique*, **339**, 165–184, doi: 10.1016/j.crme.2010.12.004.

Castro-Diaz, M. J. & Fernandez-Nieto, E. D., 2012. A class of computationally fast first order finite volume solvers: Pvm methods, https://doi.org/10.1137/100795280, **34**, doi: 10.1137/100795280.

Chen, Q., 2006. Fully nonlinear boussinesq-type equations for waves and currents over porous beds, *Journal of Engineering Mechanics*, **132**, 220–230, doi: 10.1061/(ASCE)0733-9399(2006)132:2(220).

Crosta, G. B., Imposimato, S., & Roddeman, D., 2016. Landslide spreading, impulse water waves and modelling of the vajont rockslide, *Rock Mechanics and Rock Engineering*, **49**, 2413–2436, doi: 10.1007/S00603-015-0769-Z/METRICS.

Delannay, R., Valance, A., Mangeney, A., Roche, O., & Richard, P., 2017. Granular and particle-laden flows: from laboratory experiments to field observations, *Journal of Physics D: Applied Physics*, **50**, 053001, doi: 10.1088/1361-6463/50/5/053001.

Delgado-Sánchez, J. M., Bouchut, F., Fernández-Nieto, E. D., Mangeney, A., & Narbona-Reina, G., 2020. A two-layer shallow flow model with two axes of integration, well-balanced discretization and application to submarine avalanches, *Journal of Computational Physics*, **406**, 109186, doi: 10.1016/j.jcp.2019.109186.

Escalante, C., de Luna, T. M., & Castro, M. J., 2018. Non-hydrostatic pressure shallow flows: Gpu implementation using finite volume and finite difference scheme, *Applied Mathematics and Computation*, **338**, 631–659, doi: 10.1016/j.amc.2018.06.035.
Escalante, C., Dumbser, M., & Castro, M. J., 2019. An efficient hyperbolic relaxation system for dispersive non-hydrostatic water waves and its solution with high order discontinuous galerkin schemes, *Journal of Computational Physics*, **394**, 385–416, doi: 10.1016/J.JCP.2019.05.035.

Escalante, C., Fernández-Nieto, E. D., Garres-Díaz, J., & Mangeney, A., 2023. Multilayer shallow model for dry granular flows with a weakly non-hydrostatic pressure, *Journal of Scientific Computing*, **96**, 1–50, doi: 10.1007/S10915-023-02299-Y/METRICS.

Esposti Ongaro, T., Cerminara, M., de' Michieli Vitturi, M., Tadini, A., Trolese, M., Fornaciai, A., Nannipieri, L., Calusi, B., Macías, J., Castro, M. J., Escalante, C., Ortega, S., González-Vida, J. M., & Rodríguez-Gálvez, J. F., 2025. Modeling and simulation of volcanic mass movements and induced tsunamis at Stromboli volcano (Aeolian archipelago, Tyrrhenian sea, Italy), *Frontiers in Earth Science*, **13**, 1548961, doi: 10.3389/feart.2025.1548961.

Fernández-Nieto, E. D., Bouchut, F., Bresch, D., Díaz, M. J. C., & Mangeney, A., 2008. A new Savage–Hutter type model for submarine avalanches and generated tsunami, *Journal of Computational Physics*, **227**, 7720–7754, doi: 10.1016/J.JCP.2008.04.039.

Fernández-Nieto, E. D., Garres-Díaz, J., Mangeney, A., & Narbona-Reina, G., 2016. A multilayer shallow model for dry granular flows with the rheology: application to granular collapse on erodible beds, *Journal of Fluid Mechanics*, **798**, 643–681, doi: 10.1017/JFM.2016.333.

Feuillet, N., Jorry, S., Crawford, W. C., Deplus, C., Thinon, I., Jacques, E., Saurel, J. M., Lemoine, A., Paquet, F., Satriano, C., Aiken, C., Foix, O., Kowalski, P., Laurent, A., Rinnert, E., Cathalot, C., Donval, J. P., Guyader, V., Gaillot, A., Scalabrin, C., Moreira, M., Peltier, A., Beauducel, F., Grandin, R., Ballu, V., Daniel, R., Pelleau, P., Gomez, J., Besançon, S., Geli, L., Bernard, P., Bachelery, P., Fouquet, Y., Bertil, D., Lemarchand, A., & der Woerd, J. V., 2021. Birth of a large volcanic edifice offshore mayotte via lithosphere-scale dyke intrusion, *Nature Geoscience* 2021 14:10, **14**, 787–795, doi: 10.1038/s41561-021-00809-x.

Garres-Díaz, J., Bouchut, F., Fernández-Nieto, E. D., Mangeney, A., & Narbona-Reina, G., 2020. Multilayer models for shallow two-phase debris flows with dilatancy effects, *Journal of Computational Physics*, **419**, 109699, doi: 10.1016/J.JCP.2020.109699.

Garres-Díaz, J., Fernández-Nieto, E. D., Mangeney, A., & de Luna, T. M., 2021. A weakly non-
hydrostatic shallow model for dry granular flows, *Journal of Scientific Computing*, **86**, 1–35,
doi: 10.1007/S10915-020-01377-9/METRICS.

Gonzalez-Vida, J. M., Macías, J., Castro, M. J., Sanchez-Linares, C., Asuncion, M. D. L., Ortega-Acosta, S., & Arcas, D., 2019. The lituya bay landslide-generated mega-tsunami - numerical simulation and sensitivity analysis, *Natural Hazards and Earth System Sciences*, **19**, 369–388,
doi: 10.5194/NHESS-19-369-2019.

Greve, R. & Hutter, K., 1993. Motion of a granular avalanche in a convex and concave
curved chute: experiments and theoretical predictions, *Philosophical Transactions of the
Royal Society of London. Series A: Physical and Engineering Sciences*, **342**, 573–600, doi:
10.1098/RSTA.1993.0033.

Grezio, A., Babeyko, A., Baptista, M. A., Behrens, J., Costa, A., Davies, G., Geist, E. L., Glims-
dal, S., González, F. I., Griffin, J., Harbitz, C. B., LeVeque, R. J., Lorito, S., Løvholt, F., Omira,
R., Mueller, C., Paris, R., Parsons, T., Polet, J., Power, W., Selva, J., Sørensen, M. B., & Thio,
H. K., 2017. Probabilistic tsunami hazard analysis: Multiple sources and global applications,
*Reviews of Geophysics*, **55**, 1158–1198, doi: 10.1002/2017RG000579.

Grilli, S. T., Vogelmann, S., & Watts, P., 2002. Development of a 3d numerical wave tank for modeling tsunami generation by underwater landslides, *Engineering Analysis with Boundary Elements*, **26**, 301–313, doi: 10.1016/S0955-7997(01)00113-8.

Grilli, S. T., Shelby, M., Kimmoun, O., Dupont, G., Nicolsky, D., Ma, G., Kirby, J. T., & Shi, F., 2017. Modeling coastal tsunami hazard from submarine mass failures: effect of slide rheology, experimental validation, and case studies off the us east coast, *Natural Hazards*, **86**, 353–391, doi: 10.1007/S11069-016-2692-3/METRICS.

Guan, X. & Shi, H., 2023. Translational momentum of deformable submarine landslides off a slope, *Journal of Fluid Mechanics*, **960**, doi: 10.1017/JFM.2023.177.

Gueugneau, V., Charbonnier, S., Ongaro, T. E., de' Michieli Vitturi, M., Peruzzetto, M., Man-
geney, A., Bouchut, F., Patra, A., & Kelfoun, K., 2021. Synthetic benchmarking of concentrated
pyroclastic current models, *Bulletin of Volcanology*, **83**, 75, doi: 10.1007/S00445-021-01491-Y.
Gylfadóttir, S. S., Kim, J., Helgason, J. K., Brynjólfsson, S., Ármann Höskuldsson, Jóhannesson,
T., Harbitz, C. B., & Løvholt, F., 2017. The 2014 lake askja rockslide-induced tsunami: Optimization of numerical tsunami model using observed data, *Journal of Geophysical Research: Oceans*, **122**, 4110–4122, doi: 10.1002/2016JC012496.

Harbitz, C. B., 1992. Model simulations of tsunamis generated by the storegga slides, *Marine Geology*, **105**, 1–21, doi: 10.1016/0025-3227(92)90178-K.

Harbitz, C. B., Løvholt, F., & Bungum, H., 2014. Submarine landslide tsunamis: How extreme and how likely?, *Natural Hazards*, **72**, 1341–1374, doi: 10.1007/S11069-013-0681-3/METRICS.

Hare, R., Eakins, B., & Amante, C., 2011. Modelling bathymetric uncertainty, *The International Hydrographic Review*.

Heinrich, P., Guibourg, S., Mangeney, A., & Roche, R., 1999. Numerical modeling of a landslide-generated tsunami following a potential explosion of the montserrat volcano, *Phgs. Chem. Earrh (A)*, **24**, 163–168.

Intergovernmental Panel on Climate Change (IPCC), 2023. Climate change 2022 – impacts, adaptation and vulnerability: Working group ii contribution to the sixth assessment report of the intergovernmental panel on climate change, in *Climate Change 2022 – Impacts, Adaptation and Vulnerability*, Cambridge University Press, doi: 10.1017/9781009325844.

Iverson, R. M. & George, D. L., 2014. A depth-averaged debris-flow model that includes the effects of evolving dilatancy. I. Physical basis, *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **470**(2170), 20130819, doi: 10.1098/rspa.2013.0819.

Jop, P., Forterre, Y., & Pouliquen, O., 2006. A constitutive law for dense granular flows, *Nature* 2006 **441**:7094, **441**, 727–730, doi: 10.1038/nature04801.

Kennedy, A. B., Kirby, J. T., Chen, Q., & Dalrymple, R. A., 2001. Boussinesq-type equations with improved nonlinear performance, *Wave Motion*, **33**, 225–243, doi: 10.1016/S0165-2125(00)00071-8.

Kirby, J. T., Wei, G., Chen, Q., Kennedy, A. B., & Dalrymple, R. A., 1998. Funwave 1.0: Fully nonlinear boussinesq wave model - documentation and user's manual, *research report NO. CACR-98-06*.
Kirby, J. T., Grilli, S. T., Horrillo, J., Liu, P. L., Nicolsky, D., Abadie, S., Ataie-Ashtiani, B., Castro, M. J., Clous, L., Escalante, C., Fine, I., González-Vida, J. M., Løvholt, F., Lynett, P., Ma, G., Macías, J., Ortega, S., Shi, F., Yavari-Ramshe, S., & Zhang, C., 2022. Validation and inter-comparison of models for landslide tsunami generation, *Ocean Modelling*, **170**, 101943, doi: 10.1016/J.OCEMOD.2021.101943.

Lane, E. M., Mountjoy, J. J., Power, W. L., & Mueller, C., 2016. Probabilistic hazard of tsunamis generated by submarine landslides in the cook strait canyon (new zealand), pp. 3757–3774, doi: 10.1007/978-3-319-55480-8_6.

Lemoine, A., Briole, P., Bertil, D., Roullé, A., Foumelis, M., Thinon, I., Raucoules, D., de Michele, M., Valty, P., & Colomer, R. H., 2020a. The 2018–2019 seismo-volcanic crisis east of mayotte, comoros islands: seismicity and ground deformation markers of an exceptional submarine eruption, *Geophysical Journal International*, **223**, 22–44, doi: 10.1093/GJI/GGAA273.

Lemoine, A., Pedreros, R., Filippini, A. G., Paquet, F., Thinon, I., Legendre, Y., Mompelat, P. J.-M., Arpaia, L., & Boulhaya, F., 2020b. Scénarios d’impact de tsunamis pour mayotte, *Final report RP-69869-FR*.

Leone, F., Lavigne, F., Paris, R., Denain, J. C., & Vinet, F., 2011. A spatial analysis of the december 26th, 2004 tsunami-induced damages: Lessons learned for a better risk assessment integrating buildings vulnerability, *Applied Geography*, **31**, 363–375, doi: 10.1016/J.APGEOG.2010.07.009.

Linares-Sánchez, C., Castro-Diaz, M. J., & Gonzalez-Vida, J. M., 2011. Simulación numérica de tsunamis generados por avalanchas submarinas: aplicación al caso de lituya-bay.

Liu, P. L.-F., Wu, T.-R., Raichlen, F., Synolakis, C. E., & Borrero, J. C., 2005. Runup and rundown generated by three-dimensional sliding masses, *Journal of Fluid Mechanics*, **536**, 107–144, doi: 10.1017/S0022112005004799.

Lucas, A., Mangeney, A., & Ampuero, J. P., 2014. Frictional velocity-weakening in land-slides on earth and on other planetary bodies, *Nature Communications* 2014 5:1, **5**, 1–9, doi: 10.1038/ncomms4417.

Lynett, P. & Liu, P. L., 2002. A numerical study of submarine–landslide–generated waves and
run-up, Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, **458**, 2885–2910, doi: 10.1098/RSPA.2002.0973.

Løvholt, F., Glimsdal, S., & Harbitz, C. B., 2020. On the landslide tsunami uncertainty and hazard, *Landslides*, **17**, 2301–2315, doi: 10.1007/S10346-020-01429-Z/FIGURES/9.

Løvholt, F., Glimsdal, S., & Harbitz, C. B., 2025. Tsunami hazard from subaerial landslides, *Probabilistic Tsunami Hazard and Risk Analysis: Towards Disaster Risk Reduction and Resilience*, pp. 289–301, doi: 10.1016/B978-0-443-18987-6.00010-5.

Ma, G., Kirby, J. T., Hsu, T. J., & Shi, F., 2015. A two-layer granular landslide model for tsunami wave generation: Theory and computation, *Ocean Modelling*, **93**, 40–55, doi: 10.1016/J.OCEMOD.2015.07.012.

Macías, J., Vázquez, J. T., Fernández-Salas, L. M., González-Vida, J. M., Bárcenas, P., Castro, M. J., del Río, V. D., & Alonso, B., 2015. The al-borani submarine landslide and associated tsunami. a modelling approach, *Marine Geology*, **361**, 79–95, doi: 10.1016/J.MARGEO.2014.12.006.

Macías, J., Escalante, C., Castro, M. J., González-Vida, J. M., & Ortega, S., 2017. Hysea model. landslide benchmarking results.

Macías, J., Castro, M. J., Ortega, S., & González-Vida, J. M., 2020. Performance assessment of tsunami-hysea model for nthmp tsunami currents benchmarking. field cases, *Ocean Modelling*, **152**, 101645, doi: 10.1016/J.OCEMOD.2020.101645.

Macías, J., Escalante, C., & Castro, M. J., 2021a. Multilayer-hysea model validation for landslide-generated tsunamis-part 1: Rigid slides, *Natural Hazards and Earth System Sciences*, **21**, doi: 10.5194/nhess-21-775-2021.

Macías, J., Escalante, C., & Castro, M. J., 2021b. Multilayer-hysea model validation for landslide-generated tsunamis-part 2: Granular slides, *Natural Hazards and Earth System Sciences*, **21**, doi: 10.5194/nhess-21-791-2021.

Magdalena, I., Rif'atin, H. Q., & Matin, A. M. A., 2020. Analytical and numerical studies for seiches in a closed basin with bottom friction, *Theoretical and Applied Mechanics Letters*, **10**, 429–437, doi: 10.1016/J.TAML.2020.01.057.
Mangeney, A., Heinrich, P., Rachel, R., Boudon, G., & Cheminée, J. L., 2000. Modeling of debris avalanche and generated water waves: Application to real and potential events in montserrat, *Phys. Chem. Earth (A)*, **25**, 741–745.

Mangeney, A., Bouchut, F., Thomas, N., Vilotte, J. P., & Bristeau, M. O., 2007. Numerical modeling of self-channeling granular flows and of their levee-channel deposits, *Journal of Geophysical Research: Earth Surface*, **112**, 2017, doi: 10.1029/2006JF000469.

Marras, S. & Mandli, K. T., 2021. Modeling and Simulation of Tsunami Impact: A Short Review of Recent Advances and Future Challenges, *Geosciences*, **11**(1), 5, doi: 10.3390/geosciences11010005.

Martin, H. A., Peruzzetto, M., Viroulet, S., Mangeney, A., Lagrée, P. Y., Popinet, S., Maury, B., Lefebvre-Lepot, A., Maday, Y., & Bouchut, F., 2023. Numerical simulations of granular dam break: Comparison between discrete element, navier-stokes, and thin-layer models, *Physical Review E*, **108**, 054902, doi: 10.1103/PHYSREVE.108.054902/FIGURES/17/MEDIUM.

Masson, D. G., Harbitz, C. B., Wynn, R. B., Pedersen, G., & Løvholt, F., 2006. Submarine landslides: processes, triggers and hazard prediction, *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **364**, 2009–2039, doi: 10.1098/RSTA.2006.1810.

Mathevet, T., Moine, N. L., Andréassian, V., Gupta, H., & Oudin, L., 2023. Multi-objective assessment of hydrological model performances using Nash–Sutcliffe and Kling–Gupta efficiencies on a worldwide large sample of watersheds, *Comptes Rendus. Géoscience*, **355**(S1), 117–141, doi: 10.5802/crgeos.189.

Mercury, N., Lemoine, A., Bertil, D., Vanderwoerd, J., Doubre, C., Hoste-Colomer, R., Mercury, N., Lemoine, A., Bertil, D., Vanderwoerd, J., Doubre, C., & Hoste-Colomer, R., 2020. The 2018-2020 seismo-volcanic crisis, east of mayotte, comoros islands: in-depth study of poorly instrumented first months of crisis, *AGUFM*, **2020**, V040–0007.

Mercury, N., Lemoine, A., Doubre, C., Bertil, D., Woerd, J. V. D., Hoste-Colomer, R., & Battaglia, J., 2022. Onset of a submarine eruption east of mayotte, comoros archipelago: The first ten months seismicity of the seismo-volcanic sequence (2018–2019), *Comptes Rend-
dus - Geoscience, **354**, 105-136, doi: 10.5802/CRGEOS.191/FILE/ATTACH/CRGEOS-191-SUPPL.PDF.

Mokhtarzadeh, G., Basirat, S., Bazargan, J., & Delavari, E., 2021. Impulse wave generation: A comparison of landslides of block and granular masses by coupled Lagrangian tracking using VOF over a set mesh, *Water Supply*, **22**(1), 510–526, doi: 10.2166/ws.2021.254.

Mulia, I. E., Ueda, N., Miyoshi, T., Gusman, A. R., & Satake, K., 2022. Machine learning-based tsunami inundation prediction derived from offshore observations, *Nature Communications* 2022 13:1, **13**, 1–14, doi: 10.1038/s41467-022-33253-5.

Neumann, B., Vafeidis, A. T., Zimmermann, J., & Nicholls, R. J., 2015. Future coastal population growth and exposure to sea-level rise and coastal flooding - a global assessment, *PLOS ONE*, **10**, e0118571, doi: 10.1371/JOURNAL.PONE.0118571.

Ongaro, T. E., de' Michieli Vitturi, M., Cerminara, M., Fornaciai, A., Nannipieri, L., Favalli, M., Calusi, B., Macías, J., Castro, M. J., Ortega, S., González-Vida, J. M., & Escalante, C., 2021. Modeling tsunamis generated by submarine landslides at stromboli volcano (aeolian islands, italy): A numerical benchmark study, *Frontiers in Earth Science*, **9**, 628652, doi: 10.3389/FEART.2021.628652/BIBTEX.

Paris, A., Okal, E. A., Guérin, C., Heinrich, P., Schindelé, F., & Hébert, H., 2019. Numerical modeling of the june 17, 2017 landslide and tsunami events in karrat fjord, west greenland, *Pure and Applied Geophysics*, **176**, 3035–3057, doi: 10.1007/S00024-019-02123-5/METRICS.

Pedrosa-González, M. T., González-Vida, J. M., Galindo-Záldivar, J., Ortega, S., Castro, M. J., Casas, D., & Ercilla, G., 2022. Simulation of tsunami induced by a submarine landslide in a glaciomarine margin: The case of storfjorden ls-1 (southwestern svalbard islands), *Natural Hazards and Earth System Sciences*, **22**, 3839–3858, doi: 10.5194/NHESS-22-3839-2022.

Popinet, S., 2003. Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries, *Journal of Computational Physics*, **190**, 572–600, doi: 10.1016/S0021-9991(03)00298-5.

Poulain, P., Friant, A. L., Pedreros, R., Mangeney, A., Filippini, A. G., Grandjean, G., Lemoine, A., Fernández-Nieto, E. D., Díaz, M. J. C., & Peruzzetto, M., 2022. Numerical simulation
of submarine landslides and generated tsunamis: application to the on-going mayotte seismo-volcanic crisis, *Comptes Rendus. Géoscience*, **354**, 361-390, doi: 10.5802/crgeos.138.

Poulain, P., Friant, A. L., Mangeney, A., Viroulet, S., Fernandez-Nieto, E., Diaz, M. C., Pe-ruzzetto, M., Grandjean, G., Bouchut, F., Pedreros, R., & Komorowski, J. C., 2023. Performance and limits of a shallow-water model for landslide-generated tsunamis: from laboratory exper-iments to simulations of flank collapses at montagne pelée (martinique), *Geophysical Journal International*, **233**, 796-825, doi: 10.1093/GJI/GGAC482.

Pouliquen, O. & Forterre, Y., 2002. Friction law for dense granular flows: application to the motion of a mass down a rough inclined plane, *Journal of Fluid Mechanics*, **453**, 133-151, doi: 10.1017/S0022112001006796.

Pörtner, H.-O., Roberts, D., Masson-Delmotte, V., Zhai, P., Tignor, M., Poloczanska, E., Minten-beck, K., Alegría, A., Nicolai, M., Okem, A., Petzold, J., Rama, B., & Weyer, N., 2019. Ipcc special report on the ocean and cryosphere in a changing climate, *IPCC*.

Radjai, F. & Dubois, F., 2011. Discrete-element modeling of granular materials, p. 425.

Rauter, M., Hoße, L., Mulligan, R. P., Take, W. A., & Løvholt, F., 2021. Numerical simulation of impulse wave generation by idealized landslides with openfoam, *Coastal Engineering*, **165**, 103815, doi: 10.1016/J.COASTALENG.2020.103815.

Rauter, M., Viroulet, S., Gylfadóttir, S. S., Fellin, W., & Løvholt, F., 2022. Granular porous landslide tsunami modelling – the 2014 lake askja flank collapse, *Nature Communications* 2022 **13:1**, **13**, 1–13, doi: 10.1038/s41467-022-28296-7.

Roeber, V., Cheung, K. F., & Kobayashi, M. H., 2010. Shock-capturing boussinesq-type model for nearshore wave processes, *Coastal Engineering*, **57**, 407-423, doi: 10.1016/J.COASTALENG.2009.11.007.

Roger, J., 2019. Potential tsunami hazard related to the seismic activity east of mayotte island, comoros archipelago., *Science of Tsunami Hazards*, **38**.

Roger, J. H., Bull, S., Watson, S. J., Mueller, C., Hillman, J. I., Wolter, A., Lamarche, G., Power, W., Lane, E., Woelz, S., & Davidson, S., 2024. A review of approaches for submarine landslide-tsunami hazard identification and assessment, *Marine and Petroleum Geology*, **162**, 106729,
doi: 10.1016/J.MARPETGEO.2024.106729.

Ruffini, G., Heller, V., & Briganti, R., 2021. Numerical characterisation and efficient prediction of landslide-tsunami propagation over a wide range of idealised bathymetries, *Coastal Engineering*, **167**, 103854, doi: 10.1016/j.coastaleng.2021.103854.

Sabeti, R. & Heidarzadeh, M., 2022. Numerical simulations of tsunami wave generation by submarine landslides: Validation and sensitivity analysis to landslide parameters, *Journal of Waterway, Port, Coastal, and Ocean Engineering*, **148**, 05021016, doi: 10.1061/(ASCE)WW.1943-5460.0000694/ASSET/3F278BB6-8B66-496A-9773-F990510ACB0F/ASSETS/IMAGES/LARGE/FIGURE13.JPG.

Saurel, J. M., Jacques, E., Aiken, C., Lemoine, A., Retailleau, L., Lavayssière, A., Foix, O., Dofal, A., Laurent, A., Mercury, N., Crawford, W., Lemarchand, A., Daniel, R., Pelleau, P., Berc, M. B. D., Dectot, G., Bertil, D., Roullé, A., Broucke, C., Colombain, A., Jund, H., Besançon, S., Guyavarch, P., Kowalski, P., Roudaut, M., Apprioual, R., Battaglia, J., Bodihar, S., Boissier, P., Bouin, M. P., Brunet, C., Canjamale, K., Catherine, P., Desfete, N., Doubre, C., Dretzen, R., Dumouche, T., Fernagu, P., Ferrazzini, V., Fontaine, F. R., Gaillot, A., Géli, L., Griot, C., Grunberg, M., Guzel, E. C., Hoste-Colomer, R., Lambotte, S., Lauret, F., Léger, F., Maros, E., Peltier, A., Vergne, J., Satriano, C., Tronel, F., Woerd, J. V. D., Fouquet, Y., Jorry, S. J., Rinnert, E., Thinon, I., & Feuillet, N., 2021. Mayotte seismic crisis: building knowledge in near real-time by combining land and ocean-bottom seismometers, first results, *Geophysical Journal International*, **228**, 1281–1293, doi: 10.1093/GJI/GGAB392.

Savage, S. B. & Hutter, K., 1989. The motion of a finite mass of granular material down a rough incline, *Journal of Fluid Mechanics*, **199**, 177–215, doi: 10.1017/S0022112089000340.

Scala, A., Lorito, S., Sánchez, C. E., Romano, F., Festa, G., Abbate, A., Bayraktar, H. B., Castro, M. J., Macías, J., & Gonzalez-Vida, J. M., 2024. Assessing the optimal tsunami inundation modeling strategy for large earthquakes in subduction zones, *Journal of Geophysical Research: Oceans*, **129**, e2024JC020941, doi: 10.1029/2024JC020941.

Shi, F., Kirby, J. T., Harris, J. C., Geiman, J. D., & Grilli, S. T., 2012. A high-order adaptive time-stepping tvd solver for boussinesq modeling of breaking waves and coastal inundation,
*Ocean Modelling*, **43-44**, 36–51, doi: 10.1016/J.OCEMOD.2011.12.004.

Shu, C. W. & Osher, S., 1988. Efficient implementation of essentially non-oscillatory shock-capturing schemes, *Journal of Computational Physics*, **77**, 439–471, doi: 10.1016/0021-9991(88)90177-5.

Sultan, N., Jouet, G., Riboulot, V., Terzariol, M., Garziglia, S., Cattaneo, A., Giraudeau, J., & Jorry, S. J., 2023. Sea-level fluctuations control the distribution of highly liquefaction-prone layers on volcanic-carbonate slopes, *Geology*, **51**, 402–407, doi: 10.1130/G50785.1.

Svennevig, K., Hicks, S. P., Forbriger, T., Lecocq, T., Widmer-Schnidrig, R., Mangeney, A., Hibert, C., Korsgaard, N. J., Lucas, A., Satriano, C., Anthony, R. E., Mordret, A., Schippkus, S., Rysgaard, S., Boone, W., Gibbons, S. J., Cook, K. L., Glimsdal, S., Løvholt, F., Noten, K. V., Assink, J. D., Marboeuf, A., Lomax, A., Vanneste, K., Taira, T., Spagnolo, M., Plaen, R. D., Koelemeijer, P., Ebeling, C., Cannata, A., Harcourt, W. D., Cornwell, D. G., Caudron, C., Poli, P., Bernard, P., Larose, E., Stutzmann, E., Voss, P. H., Lund, B., Cannavo, F., Castro-Díaz, M. J., Chaves, E., Dahl-Jensen, T., Dias, N. D. P., Déprez, A., Develter, R., Dreger, D., Evers, L. G., Fernández-Nieto, E. D., Ferreira, A. M. G., Funning, G., Gabriel, A.-A., Hendrickx, M., Kafka, A. L., Keiding, M., Kerby, J., Khan, S. A., Dideriksen, A. K., Lamb, O. D., Larsen, T. B., Lipovsky, B., Magdalena, I., Malet, J.-P., Myrup, M., Rivera, L., Ruiz-Castillo, E., Wetter, S., & Wirtz, B., 2024. A rockslide-generated tsunami in a greenland fjord rang earth for 9 days, *Science*, doi: 10.1126/SCIENCE.ADM9247.

Sánchez-Linares, C., de la Asunción, M., Castro, M. J., Mishra, S., & Šukys, J., 2015. Multi-level monte carlo finite volume method for shallow water equations with uncertain parameters applied to landslides-generated tsunamis, *Applied Mathematical Modelling*, **39**, 7211–7226, doi: 10.1016/J.APM.2015.03.011.

Thinon, I., Lemoine, A., Leroy, S., Paquet, F., Berthod, C., Zaragosi, S., Famin, V., Feuillet, N., Boymond, P., Masquelet, C., Mercury, N., Rusquet, A., Scalabrin, C., der Woerd, J. V., Bernard, J., Bignon, J., Clouard, V., Doubre, C., Jacques, E., Jorry, S. J., Rolandone, F., Chamot-Rooke, N., Delescluse, M., Franke, D., Watremez, L., Bachelery, P., Michon, L., Sauter, D., Bujan, S., Canva, A., Dassie, E., Roche, V., Ali, S., Allaouia, A. H. S., De-
plus, C., Rad, S., & Sadeski, L., 2022. Volcanism and tectonics unveiled in the comoros archipelago between Africa and Madagascar, *Comptes Rendus - Geoscience*, **354**, 7–34, doi: 10.5802/CRGEOS.159/FILE/ATTACH/CRGEOS-159-SUPPL.PDF.

Tinti, S., Bortolucci, E., & Romagnoli, C., 2000. Computer simulations of tsunamis due to sector collapse at stromboli, Italy, *Journal of Volcanology and Geothermal Research*, **96**, 103–128, doi: 10.1016/S0377-0273(99)00138-9.

Tozato, K., Moriguchi, S., Takase, S., Otake, Y., Motley, M. R., Suppasri, A., & Terada, K., 2023. Optimal probabilistic placement of facilities using a surrogate model for 3d tsunami simulations, *Natural Hazards and Earth System Sciences*, **23**, 1891–1909, doi: 10.5194/NHESS-23-1891-2023.

Watts, P., Grilli, S. T., Kirby, J. T., Fryer, G. J., & Tappin, D. R., 2003. Landslide tsunami case studies using a boussinesq model and a fully nonlinear tsunami generation model, *Natural Hazards and Earth System Sciences*, **3**, 391–402, doi: 10.5194/NHESS-3-391-2003.

Wei, G., Kirby, J. T., Grilli, S. T., & Subramanya, R., 1995. A fully nonlinear boussinesq model for surface waves. part 1. highly nonlinear unsteady waves, *Journal of Fluid Mechanics*, **294**, 71–92, doi: 10.1017/S0022112095002813.

Windows-Yule, C. R., Tunuguntla, D. R., & Parker, D. J., 2015. Numerical modelling of granular flows: a reality check, *Computational Particle Mechanics 2015* 3:3, **3**, 311–332, doi: 10.1007/S40571-015-0083-2.

Yavari-Ramshe, S. & Ataie-Ashtiani, B., 2016. Numerical modeling of subaerial and submarine landslide-generated tsunami waves—recent advances and future challenges, *Landslides 2016* 13:6, **13**, 1325–1368, doi: 10.1007/S10346-016-0734-2.

Zhao, H., Amann, F., & Kowalski, J., 2021. Emulator-based global sensitivity analysis for flow-like landslide run-out models, *Landslides*, **18**, 3299–3314, doi: 10.1007/S10346-021-01690-W/FIGURES/10.
<figure><img src="image_26.png" /><figcaption>**Figure A1.** Bathymetric data and isolines at 20 m, 150 m, 300 m, 450 m, 600 m, 750 m, 900 m and 1050 m depth. Cross section of the bathymetric data between the airport and east airport gauges n°1 and 6 (gauges locations described in Subsection 2.3) for the 50 m, 30 m (reference case), and 10 m grids as well as for the smoothed 30 m bathymetry. The dashed black line on the left panel indicates the cross section location. The cross section is built using the QGis open-source software.</figcaption></figure>

## APPENDIX A: SPECIFICITY OF THE MAYOTTE'S CONFIGURATION AND MANNING COEFFICIENT

Isolines of the Mayotte's bathymetry and a cross section across the coral reef are shown in Figure A1 to highlight the specificity of the Mayotte configuration. The spatialized Manning coefficient is presented in Figure A2.

## APPENDIX B: ADDITIONAL RESULTS ON THE MULTILAYER STRUCTURE AND BATHYMETRY EFFECTS

This Appendix presents the relative difference between cases related to the number of layers and bathymetric data are shown in Figures A3 and A4. The relative difference in the maximum sea surface elevation $l_{ij}$ between two simulations called Simu 1 and Simu 2 is computed by

$$ \Delta_{ij} = 100 \times \frac{l_{ij}^{\text{Simu 2}} - l_{ij}^{\text{Simu 1}}}{l_{ij}^{\text{Simu 1}}}. \quad (B.1) $$

The indices i and j stand for the grid points. A comparison of the landslide dynamics with the reference case is also shown in Figure S2 (supplementary material) for the Smoothed and the 10 m
<figure><img src="image_27.png" /><figcaption>**Figure A2.** Values of the spatialized Manning coefficient used in this study (Lemoine et al., 2020b).</figcaption></figure>

grid simulations. The maximum sea surface velocity is displayed in Figure S1 in the supplementary material.

<figure><img src="image_28.png" /><figcaption>**Figure A3.** Relative signed difference $\Delta_{ij}$ in the maximum sea surface elevation $\ell_{ij}$ (Equation (B.1)) given in percentage: a) 1 Layer (Simu 1) versus 3 Layers (Simu 2, Ref of Table 1); b) 6 Layers (Simu 2) versus 3 Layers (Simu 1, Ref of Table 1); c) 10 Layers (Simu 2) versus 6 Layers (Simu 1). Minimum and maximum differences are: a) -77% up to 500%; b) -96% up to 312%; c) -53% up to 75%.</figcaption></figure>
<figure><img src="image_29.png" /></figure>

**Figure A4.** Relative signed difference $\Delta_{ij}$ in the maximum sea surface elevation $\nu_{ij}$ given in percentage (Equation (B.1)): a) 30 m grid (Simu 2, Ref of Table 1) versus 50 m grid (Simu 1); b) 10 m grid (Simu 2) versus 30 m grid (Simu 1, Ref of Table 1); e) Smoothed (Simu 2) versus 30 m grid (Simu 1, Ref of Table 1). Minimum and maximum differences are: a) -40% up to 140%; b) -40% up to 180%; c) -50% up to 200%.

## APPENDIX C: ADDITIONAL RESULTS ON EFFECTS OF LANDSLIDE RHEOLOGY AND FRICTION PARAMETERS

Landslide maximum velocities are shown in Figures S3 (reference case) and S4 (landslide parameters' simulations) in the supplementary material. Relative differences on the maximum sea surface elevation are displayed in Figure A5 for the different landslide parameters. Figures A7 and A8 provide the maximum sea surface elevation and the inundation pattern around Petite Terre related to the grain size parameter $d$ and the rheological law $\mu$. The results related to the friction angles have been kept in the main text (Figure 16) since it presents the higher differences with the reference case among all the tested landslide parameters. The maximum landslide thickness and the mean landslide velocity for the first 50 s are presented in Figure A6 for all tested landslide parameters. Table A1 collects the runout, area and maximum thicknesses for the simulations related to the landslide parameters.
<figure><img src="image_30.png" /></figure>

**Figure A5.** Relative signed difference $\Delta_{ij}$ in the maximum sea surface elevation (Equation (B.1)) given in percentage: a) 3° friction coefficient (Simu 2) versus 11° friction coefficient (Simu 1); b) 0.1 m grain diameter (Simu 2) versus 1 m grain diameter (Simu 1, Ref of Table 1); c) Coulomb friction law (Simu 2) versus $\mu(I)$ rheology (Simu 1, Ref of Table 1); and d) $\Phi = 0.8$ (Simu 2, Ref of Table 1) versus $\Phi = 0.6$ (Simu 1). Minimum and maximum differences are: a) -97% up to 270%; b) -29% up to 51%; c) -63% up to 78%; and d) -13% up to 170%.

## APPENDIX D: ADDITIONAL NUMERICAL GAUGES DATA

Appendix D presents additional gauge data.

### D1 Tables

Table A2 collects the coordinates and depth of each numerical gauge.
<figure><img src="image_31.png" /><figcaption>**Figure A6.** Maximum landslide thickness and velocity for the first 50 s.</figcaption></figure>

## D2 Gauge sea surface elevation

The data for the hydrostatic / non-hydrostatic feature, for the grain size related parameter *d*, and for the different rheological laws have been moved to the supplementary material (Figures S6, S7 and S8). The gauge n°5 results are all gathered in Figure A9. It allows a comparison of all our simulations from Table 1. From Figure A9 and at this specific location, the parameters are classified from the most influential to the least: 1) dispersion properties (non-hydrostatic versus
<figure><img src="image_32.png" /><figcaption>**Figure A7.** a): maximum sea surface elevation $t_{ij}$ (Equation (11)) for a 0.1 m grain diameter. b): maximum inundation depth for a 0.1 m grain diameter. c): maximum sea surface elevation $t_{ij}$ (Equation (11)) for a 1 m grain diameter (reference case). d): maximum inundation depth for a 1 m grain diameter (reference case).</figcaption></figure>

hydrostatic and the layer number); 2) bathymetric data and grid resolution; 3) friction angles $\delta s$; 4) initial water level; 5) Manning coefficient; 6) rheological law type and the solid volume fraction $\Phi$; 7) mass/water friction and grain size related parameter $d$. The gauge data for friction parameters and initial water level is also presented in the supplementary material (Figure S8). The spectral amplitude at the gauge n° 1 (airport) complete the supplementary material (Figure S9).
<figure><img src="image_33.png" /><figcaption>**Figure A8.** a): maximum sea surface elevation $\iota_{ij}$ (Equation (11)) for the Coulomb friction law. b): maximum inundation depth for the Coulomb friction law. c): maximum sea surface elevation $\iota_{ij}$ (Equation (11)) for the $\mu(I)$ rheology (reference case). d): maximum inundation depth for the $\mu(I)$ rheology (reference case).</figcaption></figure>

## APPENDIX E: COMPUTATION TIME

Figure A10 presents the computation time against the number of layers and the mean error of Equation (E.1)

$$ \text{mean}(\iota_{ij}^{10 \text{ Layers}} - \iota_{ij}^{\text{x Layers}}). \qquad (\text{E.1}) $$
<table><thead><tr><th>Landslide data</th><th>Runout (km)</th><th>Area ×10<sup>7</sup> (m<sup>2</sup>)</th><th>Max final thickness (m)</th></tr></thead><tbody><tr><td>Ref</td><td>8.14</td><td>2.05</td><td>59.04</td></tr><tr><td>0.1 m grain</td><td>10.43</td><td>3.32</td><td>48.87</td></tr><tr><td>0.5 m grain</td><td>8.87</td><td>2.38</td><td>55.48</td></tr><tr><td>3° friction</td><td>12.79</td><td>3.92</td><td>51.63</td></tr><tr><td>11° friction</td><td>5.36</td><td>1.19</td><td>85.54</td></tr><tr><td>Coulomb</td><td>11.46</td><td>4.21</td><td>49.00</td></tr><tr><td>Pouliquen</td><td>8.37</td><td>2.15</td><td>54.57</td></tr><tr><td>Φ = 0.6</td><td>8.18</td><td>2.01</td><td>59.29</td></tr></tbody></table>

**Table A1.** Avalanche parameters sensitivity analysis: recap of landslide data and comparison with the reference case.

Table A3 sums up the computation time and the number of used GPUs for a few relevant simulations listed in Table 1. The other simulations take the same time as the reference. The machine used for all the calculation possesses 4 GPUs NVIDIA Tesla V100 SXM2, each GPU attached with 16Go memory.

<table><thead><tr><th></th><th>X (m)</th><th>Y (m)</th><th>Depth (m)</th></tr></thead><tbody><tr><td>Gauge 1: Airport</td><td>531228</td><td>8582886</td><td>3.2</td></tr><tr><td>Gauge 2: DEALM</td><td>524120</td><td>8585883</td><td>2.5</td></tr><tr><td>Gauge 3: Piton 200</td><td>533524</td><td>8588823</td><td>76.5</td></tr><tr><td>Gauge 4: Offshore</td><td>538926</td><td>8588820</td><td>1106</td></tr><tr><td>Gauge 5: Ponton</td><td>528057</td><td>8586964</td><td>9.8</td></tr><tr><td>Gauge 6: East Airport</td><td>533521</td><td>8582886</td><td>466</td></tr></tbody></table>

**Table A2.** Gauge locations and depth. Coordinates are given in the system EPSG:4471-RGM04/UTM zone 38S as mentioned in Subsection 2.1.
<figure><img src="image_34.png" /></figure>

**Figure A9.** Sea surface elevation at the gauge n° 5: Ponton for the different simulations. The sea surface elevation has been shifted by -1 m for the '+1 m water' case to match the hydrostatic zero of other simulations for comparison purposes.

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

**Figure A10.** Mean error compared to the 10 layers case and computational real time for simulations with different numbers of layers: a) Real computation time on 1 GPU (h); b) average error (m) of Equation (E.1).
<table><caption>**Table A3.** Computational real time for a few simulations performed in this paper.</caption><thead><tr><th></th><th>Computation time</th><th>GPU number</th></tr></thead><tbody><tr><td>Ref</td><td>8 h 40 min</td><td>1</td></tr><tr><td>Hydro</td><td>20 min</td><td>1</td></tr><tr><td>1 Layer</td><td>3 h 20 min</td><td>1</td></tr><tr><td>6 Layers</td><td>20 h</td><td>1</td></tr><tr><td>10m grid</td><td>5 days 19 h</td><td>2</td></tr><tr><td>50m grid</td><td>3 h 36 min</td><td>1</td></tr></tbody></table>

## APPENDIX F: NUMERICAL APPROXIMATION OF THE NON-HYDROSTATIC MULTILAYER MODEL

Systems (9) and (2) are solved simultaneously with a second-order HLL (Harten-Lax-van Leer), positivity-preserving, well-balanced, path-conservative finite volume numerical scheme (Macías et al., 2021b; Castro-Diaz & Fernandez-Nieto, 2012). The time step is identical for the two systems and obey the usual CFL (Courant-Friedrichs-Lewy) derived for the complete system

$$ \Delta t \max \left( \frac{1}{\Delta x} \max(|u| + \sqrt{gh}, |u_s| + \sqrt{g(1-r)h_s}), \frac{1}{\Delta y} \max(|v| + \sqrt{gh}, |v_s| + \sqrt{g(1-r)h_s}) \right) \le \frac{\text{CFL}}{2} $$

with CFL usually between 0.1 and 0.99. The CFL number is set to CFL = 0.55 for all simulations in this paper. It should be stressed that the non-hydrostatic corrections $K_{NH}$ are computed separately since they require the discretization of an elliptic operator. It is achieved through standard second-order central finite differences. Special treatments are applied in regions with small water depth (Macías et al., 2021a; Castro et al., 2005). All the lengths are dimensionalized by $h_c = 1000$ m corresponding to characteristic depth and distance for all the Mayotte's simulations.

We consider a rectangular domain $[x_a, x_b] \times [y_a, y_b]$, which is discretized by a set of rectangular cells. Let us consider a partition $\{x_i\}_{i=1}^{N_x}$ of $[x_a, x_b]$, defined by $x_i = x_a + (i - \frac{1}{2})\Delta x$, with $\Delta x = (x_b - x_a)/N_x$, being $N_x$ the number of control volumes in the X-direction. Analogously, we consider a partition $\{y_j\}_{j=1}^{N_y}$ of $[y_a, y_b]$, defined by $y_j = y_a + (j - \frac{1}{2})\Delta y$, with $\Delta y = (y_b - y_a)/N_y$, being $N_y$ the number of control volumes in the Y-direction. Then, we consider the following
definition of 2D control volumes,

$$V_{i,j} = [x_{i-1/2}, x_{i+1/2}] \times [y_{j-1/2}, y_{j+1/2}], \quad \begin{cases} i = 1, \dots, N_x \\ j = 1, \dots, N_y \end{cases}$$

being {$V_{ij}$}$_{i,j=1}^{N_x,N_y}$ a partition of the domain. Let us denote by $(\mathbf{U}_f)^n_{i,j}$ the approximation of the average of $\mathbf{U}_f(x, y, t)$ on the control volume $V_{i,j}$ at time $t^n$

$$(\mathbf{U}_f)^n_{i,j} = \begin{bmatrix} h_{i,j}^n \\ \begin{bmatrix} (h u)^n_{\alpha,i,j} \\ (h v)^n_{\alpha,i,j} \\ (h w)^n_{\alpha,i,j} \end{bmatrix}_{\{\alpha=1,\dots,L\}} \end{bmatrix}.$$

The non-hydrostatic pressure unknowns are {$p^n_{\alpha+1/2,i+1/2,j+1/2}$}$_{i,j=0}^{N_x,N_y}$, approximations of $p(x, y, z, t)$ at $z = z_{\alpha+1/2}(x, y, t)$ for $(x, y) = (x_{i+1/2}, y_{j+1/2})$, the corners of the control volumes. Note that for $\alpha = L$, it corresponds to the non-hydrostatic pressure at the free surface, set to zero.

A splitting semi-implicit projection method is considered. For the explicit part, a two-steps second order TVD method is used (Shu & Osher, 1988). The second order is achieved through reconstruction of states. For the sake of brevity, let us describe only the first order approximation of the model (9). The second order extension is done following Castro & Parés (2020). We solve the following discrete system

$$\left\{
\begin{aligned}
& \frac{(\mathbf{U}_f)_{i,j}^{n+1} - (\mathbf{U}_f)_i^n}{\Delta t} + \frac{(\mathbf{F}_f)_{i+1/2,j}^n - (\mathbf{F}_f)_{i-1/2,j}^n}{\Delta x} \\
& + \frac{(\mathbf{G}_f)_{i,j+1/2}^n - (\mathbf{G}_f)_{i,j-1/2}^n}{\Delta y} + (\mathbf{B}_f)_i^n = \\
& (\mathbf{H}_f)_i^n + (\mathbf{K}_{KH})_{i,j}^{n,n+1} - (\mathbf{S}_f)_i^{n,n+1} \\
& (\tilde{\nabla}_{NH})_{i+1/2,j+1/2}^{n,n+1} \cdot \mathbf{U}_f = \mathbf{0}
\end{aligned}
\right.
\quad (F.1)$$

Let us denote by $(\mathbf{U}_f)^{n, \pm}_{i+1/2, j}$ the hydrostatic reconstruction defined in terms of $(\mathbf{U}_f)^n_{i,j}$ and $(\mathbf{U}_f)^n_{i+1, j}$ (Audusse et al., 2006), and analogously for $(\mathbf{U}_f)^{n, \pm}_{i,j+1/2}$ in terms of $(\mathbf{U}_f)^n_{i,j}$ and $(\mathbf{U}_f)^n_{i,j+1}$. Then,
the terms that appear in the previous system are defined as follows:

$$ (\mathbf{F}_f)_{i+1/2,j}^n = \frac{1}{2}[\mathbf{F}_f((\mathbf{U}_f)_{i+1/2,j}^{n,-}) + \mathbf{F}_f((\mathbf{U}_f)_{i+1/2,j}^{n,+})] - \frac{1}{2}\left(\alpha_{0,i+1/2,j}^n[(\mathbf{U}_f)_{i+1/2,j}^{n,+} - (\mathbf{U}_f)_{i+1/2,j}^{n,-}] + \alpha_{1,i+1/2,j}^n[\mathbf{F}_f((\mathbf{U}_f)_{i+1/2,j}^{n,+}) - \mathbf{F}_f((\mathbf{U}_f)_{i+1/2,j}^{n,-}) - (\mathbf{H}_f)_{i+1/2,j}^n]\right) $$

$$ (\mathbf{G}_f)_{i,j+1/2}^n = \frac{1}{2}[\mathbf{G}_f((\mathbf{U}_f)_{i,j+1/2}^{n,-}) + \mathbf{G}_f((\mathbf{U}_f)_{i,j+1/2}^{n,+})] - \frac{1}{2}\left(\alpha_{0,i,j+1/2}^n[(\mathbf{U}_f)_{i,j+1/2}^{n,+} - (\mathbf{U}_f)_{i,j+1/2}^{n,-}] + \alpha_{1,i,j+1/2}^n[\mathbf{G}_f((\mathbf{U}_f)_{i,j+1/2}^{n,+}) - \mathbf{G}_f((\mathbf{U}_f)_{i,j+1/2}^{n,-}) - (\mathbf{H}_f)_{i,j+1/2}^n]\right). $$

The coefficients $\alpha_0$ and $\alpha_1$ defining the HLL method (Castro-Diaz & Fernandez-Nieto, 2012), are given by

$$ (\alpha_0)_{i+1/2,j}^n = \frac{S_R |S_L| - S_L |S_R|}{S_R - S_L}, \quad (\alpha_1)_{i+1/2,j}^n = \frac{|S_R| - |S_L|}{S_R - S_L}, $$

where

$$ S_L = \min\left( \begin{array}{l} \bar{u}_{i+1/2,j}^n - \sqrt{gh_{i+1/2,j}^n}, \\ (u_s)_{i+1/2,j}^n - \sqrt{g(1-r)(h_s)_{i+1/2,j}^n} \end{array} \right). $$

$$ S_R = \max\left( \begin{array}{l} \bar{u}_{i+1/2,j}^n + \sqrt{gh_{i+1/2,j}^n}, \\ (u_s)_{i+1/2,j}^n + \sqrt{g(1-r)(h_s)_{i+1/2,j}^n} \end{array} \right). $$

Coefficients $(\alpha_0)_{i,j+1/2}$ and $(\alpha_1)_{i,j+1/2}$ are defined analogously, in terms of $(v_s)_{i,j+1/2}^n$, $\bar{v}_{i,j+1/2}$, $(h_s)_{i,j+1/2}^n$ and $h_{i,j+1/2}^n$.

The approximation of the term $\mathbf{H}_f$ is the following

$$ (\mathbf{H}_f)_{i,j}^n = \frac{(\mathbf{H}_f)_{i+1/2,j}^n + (\mathbf{H}_f)_{i-1/2,j}^n}{2\Delta x} + \frac{(\mathbf{H}_f)_{i,j+1/2}^n + (\mathbf{H}_f)_{i,j-1/2}^n}{2\Delta y} $$

where

$$ (\mathbf{H}_f)_{i,j+1/2}^n = \left[ \begin{array}{c} 0 \\ 0 \\ -g \frac{h_{i,j+1/2}^{n,+} + h_{i,j+1/2}^{n,-}}{2} (h_{i,j+1/2}^{n,+} - h_{i,j+1/2}^{n,-}) \\ 0 \end{array} \right]_{\{\alpha\}}. $$

for $\alpha = 1, \dots, L$. The approximation of $\mathbf{B}_f$ is defined as follows:

$$ (\mathbf{B}_f)_{i,j}^n = \frac{(\mathbf{B}_f)_{i+1/2,j}^n + (\mathbf{B}_f)_{i-1/2,j}^n}{2\Delta x} + \frac{(\mathbf{B}_f)_{i,j+1/2}^n + (\mathbf{B}_f)_{i,j-1/2}^n}{2\Delta y} $$
with

$$ (\mathbf{B}_f)_{i+1/2,j}^n = \begin{bmatrix} 0 \\ \begin{bmatrix} (u_{\alpha+1/2,i+1/2,j}^n \Gamma_{\alpha+1/2,i+1/2,j}^n \\ -u_{\alpha-1/2,i+1/2,j}^n \Gamma_{\alpha-1/2,i+1/2,j}^n \\ (v_{\alpha+1/2,i+1/2,j}^n \Gamma_{\alpha+1/2,i+1/2,j}^n \\ -v_{\alpha-1/2,i+1/2,j}^n \Gamma_{\alpha-1/2,i+1/2,j}^n \\ (w_{\alpha+1/2,i+1/2,j}^n \Gamma_{\alpha+1/2,i+1/2,j}^n \\ -w_{\alpha-1/2,i+1/2,j}^n \Gamma_{\alpha-1/2,i+1/2,j}^n \end{bmatrix} \end{bmatrix}_{\{\alpha\}} $$

and

$$ (\mathbf{B}_f)_{i,j+1/2}^n = \begin{bmatrix} 0 \\ \begin{bmatrix} (u_{\alpha+1/2,i,j+1/2}^n \Gamma_{\alpha+1/2,i,j+1/2}^n \\ -u_{\alpha-1/2,i,j+1/2}^n \Gamma_{\alpha-1/2,i,j+1/2}^n) / l_\alpha \\ (v_{\alpha+1/2,i,j+1/2}^n \Gamma_{\alpha+1/2,i,j+1/2}^n \\ -v_{\alpha-1/2,i,j+1/2}^n \Gamma_{\alpha-1/2,i,j+1/2}^n) / l_\alpha \\ (w_{\alpha+1/2,i,j+1/2}^n \Gamma_{\alpha+1/2,i,j+1/2}^n \\ -w_{\alpha-1/2,i,j+1/2}^n \Gamma_{\alpha-1/2,i,j+1/2}^n) / l_\alpha \end{bmatrix} \end{bmatrix}_{\{\alpha\}} $$

for $\alpha = 1, \dots, L$. The terms in $(\mathbf{B}_f)_{i+1/2,j}^n$ and $(\mathbf{B}_f)_{i,j+1/2}^n$ are given by

$$ \Gamma_{1/2,i+1/2,j}^n = \Gamma_{L+1/2,i+1/2,j}^n = 0, $$

$$ \Gamma_{1/2,i,j+1/2}^n = \Gamma_{L+1/2,i,j+1/2}^n = 0, $$

$$ \Gamma_{\alpha+1/2,i+1/2,j}^n = \sum_{\beta=1}^{\alpha} l_{\beta} \left( h_{i+1/2,j}^{n,+} (\bar{u}_{i+1/2,j}^{n,+} - u_{\beta,i+1/2,j}^{n,+}) - h_{i+1/2,j}^{n,-} (\bar{u}_{i+1/2,j}^{n,-} - u_{\beta,i+1/2,j}^{n,-}) \right), $$

$$ \Gamma_{\alpha+1/2,i,j+1/2}^n = \sum_{\beta=1}^{\alpha} l_{\beta} \left( h_{i,j+1/2}^{n,+} (\bar{v}_{i,j+1/2}^{n,+} - v_{\beta,i,j+1/2}^{n,+}) - h_{i,j+1/2}^{n,-} (\bar{v}_{i,j+1/2}^{n,-} - v_{\beta,i,j+1/2}^{n,-}) \right), $$

for $\alpha = 1, \dots, L-1$, with $u_{\alpha+1/2,i+1/2,j}^n = (u_{\alpha,i+1/2,j}^{n,+} + u_{\alpha,i+1/2,j}^{n,-} + u_{\alpha+1,i+1/2,j}^{n,+} + u_{\alpha+1,i+1/2,j}^{n,-})/4$

and analogous definitions for $v_{\alpha+1/2,i+1/2,j}^n, w_{\alpha+1/2,i+1/2,j}^n, v_{\alpha+1/2,i,j+1/2}^n$ and $w_{\alpha+1/2,i,j+1/2}^n$.
The discretization of the non-hydrostatic gradient pressure terms is given by

$$ (\mathbf{K}_{NH})_{i,j}^{n,n+1} = \begin{bmatrix} 0 \\ (K_x)_{\alpha,i,j}^{n,n+1} \\ (K_y)_{\alpha,i,j}^{n,n+1} \\ -(p_{\alpha+1/2,i,j}^{n+1} - p_{\alpha-1/2,i,j}^{n+1})/l_{\alpha} \end{bmatrix}_{\{\alpha\}} $$

for $\alpha = 1, \dots, L$, where

$$ (K_x)_{\alpha,i,j}^{n,n+1} = \left( -h_{i+1/2,j}^n p_{\alpha,i+1/2,j}^{n+1} + h_{i-1/2,j}^n p_{\alpha,i-1/2,j}^{n+1} \right. \\ \left. \qquad + p_{\alpha+1/2,i,j}^{n+1} (z_{\alpha+1/2,i+1/2,j} - z_{\alpha+1/2,i-1/2,j})/l_{\alpha} \right. \\ \left. \qquad - p_{\alpha-1/2,i,j}^{n+1} (z_{\alpha-1/2,i+1/2,j} - z_{\alpha-1/2,i-1/2,j})/l_{\alpha} \right) / (\Delta x), $$

$$ (K_y)_{i,j}^{\alpha,n,n+1} = \left( -h_{i,j+1/2}^n p_{\alpha,i,j+1/2}^{n+1} + h_{i,j-1/2}^n p_{\alpha,i,j-1/2}^{n+1} \right. \\ \left. \qquad + p_{\alpha+1/2,i,j}^{n+1} (z_{\alpha+1/2,i,j+1/2} - z_{\alpha+1/2,i,j-1/2})/l_{\alpha} \right. \\ \left. \qquad - p_{\alpha-1/2,i,j}^{n+1} (z_{\alpha-1/2,i,j+1/2} - z_{\alpha-1/2,i,j-1/2})/l_{\alpha} \right) / (\Delta y), $$

$$ p_{\alpha,i+1/2,j}^{n+1} = (p_{\alpha+1/2,i+1/2,j} + p_{\alpha-1/2,i+1/2,j})/2 \quad \text{and} \quad p_{\alpha+1/2,i,j} = (p_{\alpha+1/2,i-1/2,j-1/2} + p_{\alpha+1/2,i-1/2,j+1/2} + p_{\alpha+1/2,i+1/2,j-1/2} + p_{\alpha+1/2,i+1/2,j+1/2})/4. $$

ORIGINAL UNEDITED MANUSCRIPT
The friction and breaking terms, defined in $S_f$, are discretized semi-implicitly

$$ (S_f)_{i,j}^{n,n+1} = \left[ \begin{array}{l} 0 \\ g h_{i,j}^n \frac{n^2}{l_1 (h_{i,j}^n)^{4/3}} u_{1,i,j}^{n+1} \sqrt{(u_{1,i,j}^n)^2 + (v_{1,i,j}^n)^2} \\ -m_f (u_{s,i,j}^{n+1} - u_{1,i,j}^n) \\ \\ g h_{i,j}^n \frac{n^2}{l_1 (h_{i,j}^n)^{4/3}} v_{1,i,j}^{n+1} \sqrt{(u_{1,i,j}^n)^2 + (v_{1,i,j}^n)^2} \\ -m_f (v_{s,i,j}^{n+1} - v_{1,i,j}^n) \\ \\ m_1^{w,n,n+1} \\ \left[ \begin{array}{l} 0 \\ 0 \\ m_\alpha^{w,n,n+1} \end{array} \right]_{\{\alpha=2,...,L\}} \end{array} \right], $$

with

$$ m_{\alpha}^{w,n,n+1} = \frac{C_{i,j}^{n}}{l_{\alpha}^{2}} G_{\alpha,i,j}^{n} w_{\alpha,i,j}^{n+1}, $$

$$ G_{\alpha,i,j}^{n} = \left( \left( \frac{(hu)^{n}_{\alpha,i+1,j} - (hu)^{n}_{\alpha,i-1,j}}{2\Delta x} \right)^{2} + \left( \frac{(hv)^{n}_{\alpha,\bar{u},j+1} - (hv)^{n}_{\alpha,\alpha,j-1}}{2\Delta x} \right)^{2} \right)^{1/2} $$

$$ \text{for } \alpha = 1, \dots, N, \text{ and } C_{i,j}^{n} = 35 \max \left( \frac{\sqrt{(\bar{u}_{i,j}^{n})^2 + (\bar{v}_{i,j}^{n})^2}}{0.4 \sqrt{gh_{i,j}^{n}}} - 1, 0 \right). $$

Finally, a projection method is considered, which is based on a semi-implicit discretization of the divergence operator at the corners of the control volumes $x_{i+1/2,j+1/2}$ for $i = 0, \dots, N_x$ and $j = 0, \dots, N_y$. The discretization of the divergence type operator associated to the multilayer
system is

$$h_{i+1/2,j+1/2}^n (\tilde{\nabla}_{NH})_{i+1/2,j+1/2}^{n,n+1} \cdot U_f = \\ \left[ \begin{array}{l} l_1 h_{i+1/2,j+1/2}^n ((hu)_{1,i+1,j+1/2}^{n+1} - (hu)_{1,i,j+1/2}^{n+1}) / \Delta x - l_1 (hu)_{i+1/2,j+1/2}^{n+1} (h_{i+1,j+1/2}^n - h_{i,j+1/2}^n) / \Delta x \\ \\ + l_1 h_{i+1/2,j+1/2}^n ((hv)_{1,i+1/2,j+1}^{n+1} - (hv)_{1,i+1/2,j}^{n+1}) / \Delta y - l_1 (hv)_{i+1/2,j+1/2}^{n+1} (h_{i+1/2,j+1}^n - h_{i+1/2,j}^n) / \Delta y \\ \\ + 2((hw)_{1,i+1/2,j+1/2}^{n+1} - (hw)_{0,i+1/2,j+1/2}^{n+1}) \\ \\ \left[ \begin{array}{l} ((hu)_{\alpha+1,i+1/2,j+1/2} - (hu)_{\alpha,i+1/2,j+1/2}) (z_{\alpha+1/2,i+1,j+1/2} - z_{\alpha+1/2,i,j+1/2}) / \Delta x \\ \\ + ((hv)_{\alpha+1,i+1/2,j+1/2} - (hv)_{\alpha,i+1/2,j+1/2}) (z_{\alpha+1/2,i+1/2,j+1} - z_{\alpha+1/2,i+1/2,j}) / \Delta y \\ \\ - ((hw)_{\alpha+1,i+1/2,j+1/2}^{n+1} - (hw)_{\alpha,i+1/2,j+1/2}^{n+1}) \end{array} \right]_{\{\alpha\}} \end{array} \right]$$

for $\alpha = 1, \dots, L-1$, where $w_0$ is defined as an approximation of $\partial_t(-H + h_s)$. In the previous equation, the following notation has been considered: $(hu)_{\alpha,i+1/2,j+1/2} = ((hu)_{\alpha,i,j} + (hu)_{\alpha,i+1,j} + (hu)_{\alpha,i,j+1} + (hu)_{\alpha,i+1,j+1})/4$, $(hu)_{\alpha,i,j+1/2} = ((hu)_{\alpha,i,j} + (hu)_{\alpha,i,j+1})/2$, and analogously for $h$, $(hv)$ and $(hw)$. The non-hydrostatic pressure approximation are obtained by solving a linear system. It is deduced by substituting the expression of $U_f^{n+1}$ into the previous equation associated to the discretization of the velocity divergence.

This paper has been produced using the Blackwell Scientific Publications GJI LATEX2e class file.