Multilayer-HySEA model validation for landslide-generated tsunamis – Part 2: Granular slides

The final aim of the present work is to propose a NTHMP-benchmarked numerical tool for landslide-generated tsunami hazard assessment. To achieve this, the novel Multilayer-HySEA model is validated using laboratory experiment data for landslide-generated tsunamis. In particular, this second part of the work deals with granular slides, while the first part, in a companion paper, considers rigid slides. The experimental data used have been proposed by the US National Tsunami Hazard and Mitigation Program (NTHMP) and were established for the NTHMP Landslide Benchmark Workshop, held in January 2017 at Galveston (Texas). Three of the seven benchmark problems proposed in that workshop dealt with tsunamis generated by rigid slides and are collected in the companion paper (Macías et al., 2021). Another three benchmarks considered tsunamis generated by granular slides. They are the subject of the present study. The seventh benchmark problem proposed the field case of Port Valdez, Alaska, 1964 and can be found in Macías et al. (2017). In order to reproduce the laboratory experiments dealing with granular slides, two models need to be coupled: one for the granular slide and a second one for the water dynamics. The coupled model used consists of a new and efficient hybrid finite-volume–finite-difference implementation on GPU architectures of a non-hydrostatic multilayer model coupled with a Savage–Hutter model. To introduce the multilayer model more fluidly, we first present the equations of the one-layer model, Landslide-HySEA, with both strong and weak couplings between the fluid layer and the granular slide. Then, a brief description of the multilayer model equations and the numerical scheme used is included. The dispersive properties of the multilayer model can be found in the companion paper. Then, results for the three NTHMP benchmark problems dealing with tsunamis generated by granular slides are presented with a description of each benchmark problem.

1 Introduction

Following the introduction of the companion paper Macías et al. (2021), a landslide tsunami model benchmarking and validation workshop was held on 9–11 January 2017 in Galveston, TX. This workshop was organized on behalf of the NOAA National Weather Service's National Tsunami Hazard Mitigation Program (NTHMP) Mapping and Modeling Subcommittee (MMS) with the expected outcome being to (i) develop a set of community-accepted benchmark tests for validating models used for landslide tsunami generation and propagation in NTHMP inundation mapping work; (ii) develop workshop documentation and a web-based repository, for benchmark data, model results, and workshop documentation, results, and conclusions; and (iii) provide recommendations as a basis for developing best-practice guidelines for landslide tsunami modeling in NTHMP work.

A set of seven benchmark tests was selected (Kirby et al., 2018). The selected benchmarks were taken from a subset of available laboratory data sets for solid slide experiments (three of them) and deformable slide experiments (another three) that included both submarine and subaerial slides. Finally, a benchmark based on a historic field event (Valdez, AK, 1964) closed the list of proposed benchmarks. The EDANYA group (https://www.uma.es/edanya, last access: 21 February 2021) from the University of Malaga participated in the aforementioned workshop, and the numerical codes Multilayer-HySEA and Landslide-HySEA were used to produce our modeled results. We presented numerical results for six out of the seven benchmark problems proposed, including the field case (Macías et al., 2017). The sole benchmark we did not perform at the time was BP6, for which numerical results are included here.

The present work aims at showing the numerical results obtained with the Multilayer-HySEA model in the framework of the validation effort described above for the case of granular-slide-generated tsunamis for the complete set of the three benchmark problems proposed by the NTHMP. However, the ultimate goal of the present work is to provide the tsunami community with a numerical tool, tested and validated meeting the defined criteria for the NTHMP, for landslide-generated tsunami hazard assessment. This NTHMP acceptance has already been achieved by the Tsunami-HySEA model for the case of earthquake-generated tsunamis (Macías et al., 2017, 2020a, c).

Fifteen years ago, at the beginning of the century, solid block landslide modeling challenged researchers and was undertaken by a number of authors (see companion paper, Macías et al., 2021, for references), and laboratory experiments were developed for those cases and for tsunami model benchmarking. In contrast, some early models (e.g., Heinrich, 1992; Harbitz et al., 1993; Rzadkiewicz et al., 1997; Fine et al., 1998) and a number of more recent models have simulated tsunami generation by deformable slides, based either on depth-integrated two-layer model equations or on solving more complete sets of equations in terms of featured physics (dispersive, non-hydrostatic, Navier–Stokes). Examples include solutions of 2D or 3D Navier–Stokes equations to simulate subaerial or submarine slides modeled as dense Newtonian or non-Newtonian fluids (Ataie-Ashtiani and Shobeyri, 2008; Weiss et al., 2009; Abadie et al., 2010, 2012; Horrillo et al., 2013), flows induced by sediment concentration (Ma et al., 2013), or fluid or granular flow layers penetrating or failing underneath a 3D water domain – for example, the two-layer models of Macías et al. (2015) or González-Vida et al. (2019), where a fully coupled non-hydrostatic SW/Savage–Hutter model is used, or the model used in Ma et al. (2015) and Kirby et al. (2016), in which the upper water layer is modeled with the non-hydrostatic σ-coordinate 3D model NHWAVE (Ma et al., 2012). For a more comprehensive review of recent modeling work, see Yavari-Ramshe and Ataie-Ashtiani (2016). A number of recent laboratory experiments have modeled tsunamis generated by subaerial landslides composed of gravel (Fritz et al., 2004; Ataie-Ashtiani and Najafi-Jilani, 2008; Heller and Hager, 2010; and Mohammed and Fritz, 2012) or glass beads (Viroulet et al., 2014). For deforming underwater landslides and related tsunami generation, 2D experiments were performed by Rzadkiewicz et al. (1997), who used sand, and Ataie-Ashtiani and Najafi-Jilani (2008), who used granular material. Well-controlled 2D glass bead experiments were reported and modeled by Grilli et al. (2017) using the model of Kirby et al. (2016).

The benchmark problems performed in the present work are based on the laboratory experiments of Kimmoun and Dupont (see Grilli et al., 2017) for BP4, Viroulet et al. (2014) for BP5, and Mohammed and Fritz (2012) for BP6. The basic reference for these three benchmarks, but also the three related to solid slides and the Alaska field case, all of them proposed by the NTHMP, is Kirby et al. (2018). That is a key reference for readers interested in the benchmarking initiative which the present work is based on.

2 The Landslide-HySEA model for granular slides

First we consider the Landslide-HySEA model, applied in Macías et al. (2015) and González-Vida et al. (2019), which for the case of one-dimensional domains reads

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}{\partial }_{t}h+{\partial }_x\left(hu\right)=\mathrm{0},\\ {\partial }_t\left(hu\right)+{\partial }_x\left(h{u}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\right)-gh{\partial }_x\left(H-z_s\right)\\ =n_a\left(u_s-u\right),\\ {\partial }_t{z}_s+{\partial }_x\left(z_s{u}_s\right)=\mathrm{0},\\ {\partial }_t\left(z_s{u}_s\right)+{\partial }_x\left(z_s{u}_s^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g(\mathrm{1}-r)z_s^{\mathrm{2}}\right)\\ =g{z}_s{\partial }_x\left((\mathrm{1}-r)H-r\mathit{\eta }\right)-r{n}_a\left(u_s-u\right)+{\mathit{\tau }}_P,\end{array}\right.\end{array}$$

where g is the gravity acceleration (g=9.81 m s⁻²); H(x) is the non-erodible (does not evolve in time) bathymetry measured from a given reference level; z_s(x, t) represents the thickness of the layer of granular material at each point x at time t; h(x, t) is the total water depth; η(x, t) denotes the free surface (measured form the same fixed reference level used for the bathymetry, for example, the mean sea surface) and is given by $\mathit{\eta }=h+z_s-H$; u(x, t) and u_s(x, t) are the averaged horizontal velocity for the water and for the granular material, respectively; and $r=\frac{{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{2}}}$ is the ratio of densities between the ambient fluid and the granular material. The term n_a(u_s−u) parameterizes the friction between the fluid and the granular layer. Finally, the term τ_P(x, t) represents the friction between the granular slide and the non-erodible bottom surface. It is parameterized as in Pouliquen and Forterre (2002), and it is described in the next section.

System Eq. (1) presents the difficulty of considering the complete coupling between sediment and water, including the corresponding coupled pressure terms. That makes its numerical approximation more complex. Moreover, it also makes it difficult to consider its natural extension to non-hydrostatic flows.

Now, if ∂_xη is neglected in the momentum equation of the granular material, that is, the fluctuation of pressure due to the variations of the free surface is neglected in the momentum equation of the granular material, then the following weakly coupled system could be obtained:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{S}\text{-}\mathrm{W}\mathrm{system}\left\{\begin{array}{l}{\partial }_{t}h+{\partial }_x\left(hu\right)=\mathrm{0},\\ {\partial }_t\left(hu\right)+{\partial }_x\left(h{u}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\right)\\ -gh{\partial }_x\left(H-z_s\right)=n_a\left(u_s-u\right),\end{array}\right.\text{(3)}& {\displaystyle }\mathrm{S}\text{-}\mathrm{H}\mathrm{system}\left\{\begin{array}{l}{\partial }_t{z}_s+{\partial }_x\left(z_s{u}_s\right)=\mathrm{0},\\ {\partial }_t\left(z_s{u}_s\right)+{\partial }_x\left(z_s{u}_s^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g(\mathrm{1}-r)z_s^{\mathrm{2}}\right)\\ -g(\mathrm{1}-r)z_s{\partial }_{x}H=-r{n}_a\left(u_s-u\right)+{\mathit{\tau }}_P,\end{array}\right.\end{array}$$

where the first system is the standard one-layer shallow-water system and the second one is the one-layer reduced-gravity Savage–Hutter model (Savage and Hutter, 1989) that takes into account that the granular landslide is underwater. Note that the previous system could be also adapted to simulate subaerial/submarine landslides by a suitable treatment of the variation of the gravity terms. Under this formulation, it is now straightforward to improve the numerical model for the fluid phase by including non-hydrostatic effects.

In the present study, the governing equations of the landslide motion are derived in Cartesian coordinates. In some cases where steep slopes are involved, landslide models based on local coordinates allow representing the slide motion better. However, when general topographies are considered and not only simple geometries, landslide models based on local coordinates also introduce some difficulties on the final numerical model and on its implementation compromising, at the same time, the computational efficiency of the numerical model. Here, we focus on the hydrodynamic component of the system, and that is one of the reasons for choosing a simple landslide model based on Cartesian coordinates. Of course, the strategies presented here can also be adapted for more sophisticated landslide models. For example, in Garres-Díaz et al. (2021) a non-hydrostatic model for the hydrodynamic part that is similar to the one presented here for the case of a single layer was introduced. In the work mentioned above, the authors study the influence of coupling the hydrodynamic model with a granular model that is derived in both reference systems: Cartesian and local coordinates. The front positions calculated with the Cartesian model progress faster, and, after some time, they are slightly ahead compared with the local coordinate model solution (see, for instance, Fig. 4 in Garres-Díaz et al., 2021). This is due to the fact that the Cartesian model uses the horizontal velocity instead of the velocity tangent to the topography. In any case, the differences between the two models are not very noticeable. A granular slide model based on local coordinates might give better results. However, when introducing a non-hydrostatic pressure, the model is closer to a 3D solver. In such a case, the influence on the reference coordinate system barely exists. That is the reason why in Garres-Díaz et al. (2021) both non-hydrostatic models based on different coordinate systems show similar results. In any case, although on the present work we focus on the hydrodynamic part, it can be observed on the benchmark tests that the numerical results are in very good agreement with the laboratory-measured data, despite the simple landslide model chosen here.

3 The Multilayer-HySEA model

The Multilayer-HySEA model implements a two-phase model intended to reproduce the interaction between the slide granular material (submarine or subaerial) and the fluid. In the present work, a multilayer non-hydrostatic shallow-water model is considered for modeling the evolution of the ambient water (see Fernández-Nieto et al., 2018), and for simulating the kinematics of the submarine/subaerial landslide the Savage–Hutter model (Eq. 3) is used. The coupling between these two models is performed through the boundary conditions at their interface. The parameter r represents the ratio of densities between the ambient fluid and the granular material (slide liquefaction parameter). Usually

$$\begin{array}{}\text{(4)}& r={\displaystyle \frac{{\mathit{\rho }}_f}{{\mathit{\rho }}_b}},{\mathit{\rho }}_b=(\mathrm{1}-\mathit{\phi }){\mathit{\rho }}_s+\mathit{\phi }{\mathit{\rho }}_f,\end{array}$$

where ρ_s stands for the typical density of the granular material, ρ_f is the density of the fluid (ρ_s>ρ_f), both constant, and φ represents the porosity ($\mathrm{0}\le \mathit{\phi }<\mathrm{1}$). In the current work, the porosity, φ, is supposed to be constant in space and time, and, therefore, the ratio r is also constant. This ratio ranges from 0 to 1 (i.e., $\mathrm{0}<r<\mathrm{1}$) and, even on a uniform material, is difficult to estimate as it depends on the porosity (and ρ_f and ρ_s are also supposed to be constant). Typical values for r are in the interval   [0.3, 0.8].

Figure 1 Schematic diagram describing the multilayer system.

The fluid model

The ambient fluid is modeled by a multilayer non-hydrostatic shallow-water system (Fernández-Nieto et al., 2018) to account for dispersive water waves. The model considered, that is obtained by a process of depth-averaging of the Euler equations, can be interpreted as a semi-discretization with respect to the vertical coordinate. In order to take into account dispersive effects, the total pressure is decomposed into the sum of hydrostatic and non-hydrostatic components. In this process, the horizontal and vertical velocities are supposed to have constant vertical profiles. The resulting multilayer model admits an exact energy balance, and when the number of layers increases, the linear dispersion relation of the linear model converges to that of Airy's theory. Finally, the model proposed in Fernández-Nieto et al. (2018) can be written in compact form as

$$\begin{array}{}\text{(5)}& \left\{\begin{array}{l}{\partial }_{t}h+{\partial }_x\left(hu\right)=\mathrm{0},\\ {\partial }_t\left(h{u}_{\mathit{\alpha }}\right)+{\partial }_x\left(h{u}_{\mathit{\alpha }}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\right)-gh{\partial }_x\left(H-z_s\right)\\ +u_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}-u_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}\\ =-h\left({\partial }_x{p}_{\mathit{\alpha }}+{\mathit{\sigma }}_{\mathit{\alpha }}{\partial }_z{p}_{\mathit{\alpha }}\right)-{\mathit{\tau }}_{\mathit{\alpha }}\\ {\partial }_t\left(h{w}_{\mathit{\alpha }}\right)+{\partial }_x\left(h{u}_{\mathit{\alpha }}{w}_{\mathit{\alpha }}\right)+w_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}\\ -w_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=-h{\partial }_z{p}_{\mathit{\alpha }},\\ {\partial }_x{u}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}+{\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}{\partial }_z{u}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}+{\partial }_z{w}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}=\mathrm{0},\end{array}\right.\end{array}$$

for α∈{1, 2, …, L}, with L the number of layers and where the following notation has been used:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{f}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(f_{\mathit{\alpha }+\mathrm{1}}+f_{\mathit{\alpha }}\right),\\ {\displaystyle }& {\displaystyle }{\partial }_z{f}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{h\mathrm{\Delta }s}}\left(f_{\mathit{\alpha }+\mathrm{1}}-f_{\mathit{\alpha }}\right),\end{array}$$

where f denotes one of the generic variables of the system, i.e., u, w, and p; $\mathrm{\Delta }s=\mathrm{1}/L$ and, finally,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathit{\alpha }}={\partial }_x\left(H-z_s-h\mathrm{\Delta }s(\mathit{\alpha }-\mathrm{1}/\mathrm{2})\right),\\ {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_{\mathit{\alpha }-\mathrm{1}/\mathrm{2}}={\partial }_x\left(H-z_s-h\mathrm{\Delta }s(\mathit{\alpha }-\mathrm{1})\right).\end{array}$$

Figure 1 shows a schematic picture of the model configuration, where the total water height h is decomposed along the vertical axis into L≥1 layers. The depth-averaged velocities in the x and z directions are written as u_α and w_α, respectively. The non-hydrostatic pressure at the interface $z_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}$ is denoted by $p_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}$. The free-surface elevation measured from a fixed reference level (for example the still-water level) is written as η and $\mathit{\eta }=h-H+z_s$, where again H(x) is the unchanged non-erodible bathymetry measured from the same fixed reference level. τ_α=0 for α>1, and τ₁ is given by

$${\mathit{\tau }}_{\mathrm{1}}={\mathit{\tau }}_b-n_a\left(u_s-u_{\mathrm{1}}\right),$$

where τ_b stands for an classical Manning-type parameterization for the bottom shear stress and, in our case, is given by

$${\mathit{\tau }}_b=gh{\displaystyle \frac{n^{\mathrm{2}}}{h^{\mathrm{4}/\mathrm{3}}}}{u}_{\mathrm{1}}\left|u_{\mathrm{1}}\right|,$$

and n_a(u_s−u₁) accounts for the friction between the fluid and the granular layer. The latest two terms are only present at the lowest layer (α=1). Finally, for α=1, …, L−1, ${\mathrm{\Gamma }}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}$ parameterizes the mass transfer across interfaces, and those terms are defined by

$${\mathrm{\Gamma }}_{\mathit{\alpha }+\mathrm{1}/\mathrm{2}}=\sum _{\mathit{\beta }=\mathit{\alpha }+\mathrm{1}}^L{\partial }_x\left(h\mathrm{\Delta }s\left(u_{\mathit{\beta }}-\stackrel{\mathrm{‾}}{u}\right)\right),\stackrel{\mathrm{‾}}{u}=\sum _{\mathit{\alpha }=\mathrm{1}}^L\mathrm{\Delta }s{u}_{\mathit{\alpha }}.$$

Here we suppose that ${\mathrm{\Gamma }}_{\mathrm{1}/\mathrm{2}}={\mathrm{\Gamma }}_{L+\mathrm{1}/\mathrm{2}}=\mathrm{0}$; this means that there is no mass transfer through the sea floor or the water free surface. In order to close the system, the boundary condition

$$p_{L+\mathrm{1}/\mathrm{2}}=\mathrm{0}$$

is imposed at the free-surface and the boundary conditions

$$u_{\mathrm{0}}=\mathrm{0},w_{\mathrm{0}}=-{\partial }_t\left(H-z_s\right)$$

are imposed at the bottom. The last two conditions enter into the incompressibility relation for the lowest layer (α=1), given by

$${\partial }_x{u}_{\mathrm{1}/\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{1}/\mathrm{2}}{\partial }_z{u}_{\mathrm{1}/\mathrm{2}}+{\partial }_z{w}_{\mathrm{1}/\mathrm{2}}=\mathrm{0}.$$

It should be noted that both models, the hydrodynamic model described here and the morphodynamic model described in the next subsection, are coupled through the unknown z_s, which, in the case of the model described here, is present in the equations and in the boundary condition ($w_{\mathrm{0}}=-{\partial }_t(H-z_s)$).

Some dispersive properties of system Eq. (5) were originally studied in Fernández-Nieto et al. (2018). Moreover, for a better-detailed study on the dispersion relation (such as phase velocity, group velocity, and linear shoaling) the reader is referred to the companion paper Macías et al. (2021).

Along the derivation of the hydrodynamic model presented here, the rigid-lid assumption for the free surface of the ambient fluid is adopted. This means that pressure variations induced by the fluctuation on the free surface of the ambient fluid over the landslide are neglected.

The landslide model

The 1D Savage–Hutter model used and implemented in the present work is given by system Eq. (3). The friction law τ_P (Pouliquen and Forterre, 2002) is given by the expression

$${\mathit{\tau }}_P=-g(\mathrm{1}-r)\mathit{\mu }{z}_s{\displaystyle \frac{u_s^{\mathrm{2}}}{\left|u_s\right|}},$$

where μ is a constant friction coefficient with a key role, as it controls the movement of the landslide. Usually μ is given by the Coulomb friction law as the simpler parameterization that can be used in landslide models. However, it is well known that a constant friction coefficient does not allow us to reproduce steady uniform flows over rough beds observed in the laboratory for a range of inclination angles. To reproduce these flows, in Pouliquen and Forterre (2002), the authors introduce an empirical friction coefficient μ that depends on the norm of the mean velocity u_s, on the thickness z_s of the granular layer, and on the Froude number $Fr=\frac{u_s}{\sqrt{g{z}_s}}$. The friction law is given by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\mu }\left(z_s,u_s\right)=\\ {\displaystyle }& {\displaystyle }\left\{\begin{array}{ll}{\mathit{\mu }}_{\mathrm{start}}\left(z_s\right)+{\left(\frac{Fr}{\mathit{\beta }}\right)}^{\mathit{\gamma }}\left({\mathit{\mu }}_{\mathrm{stop}}\left(z_s\right)-{\mathit{\mu }}_{\mathrm{start}}\left(z_s\right)\right),& \mathrm{for}Fr<\mathit{\beta },\\ {\mathit{\mu }}_{\mathrm{stop}}\left(z_s\right),& \mathrm{for}\mathit{\beta }\le Fr,\end{array}\right.\end{array}$$

with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\mu }}_{\mathrm{start}}\left(z_s\right)=\tan \left({\mathit{\delta }}_{\mathrm{3}}\right)+\left(\tan \left({\mathit{\delta }}_{\mathrm{2}}\right)-\tan \left({\mathit{\delta }}_{\mathrm{1}}\right)\right)\exp \left(-{\displaystyle \frac{z_s}{d_s}}\right),\\ {\displaystyle }& {\displaystyle }{\mathit{\mu }}_{\mathrm{stop}}\left(z_s\right)=\tan \left({\mathit{\delta }}_{\mathrm{1}}\right)+\left(\tan \left({\mathit{\delta }}_{\mathrm{2}}\right)-\tan \left({\mathit{\delta }}_{\mathrm{1}}\right)\right)\exp \left(-{\displaystyle \frac{z_s\mathit{\beta }}{d_{s}Fr}}\right),\end{array}$$

where d_s represents the mean size of grains. β=0.136 and $\mathit{\gamma }={\mathrm{10}}^{-\mathrm{3}}$ are empirical parameters. tan (δ₁) and tan (δ₂) are the characteristic angles of the material, and tan (δ₃) is other friction angle related to the behavior when starting from rest. This law has been widely used in the literature (see, e.g., Brunet et al., 2017).

Note that the slide model can also be adapted to simulate subaerial landslides. The presence of the term (1−r) in the definition of the Pouliquen–Forterre friction law is due to the buoyancy effects, which must be taken into account only in the case that the granular material layer is submerged in the fluid. Otherwise, this term must be replaced by 1.

4 Numerical solution method

System Eq. (3) can be written in the following compact form:

$$\begin{array}{}\text{(6)}& {\partial }_t{U}_s+{\partial }_x{F}_s\left(U_s\right)=G_s\left(U_s\right){\partial }_{x}H-S_s\left(U_s\right),\end{array}$$

being

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{U}_s=\left[\begin{array}{c}{z}_s\\ u_s{z}_s\end{array}\right],F_s\left(U_s\right)=\left[\begin{array}{c}{z}_s{u}_s\\ z_s{u}_s^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g(\mathrm{1}-r)z_s^{\mathrm{2}}\end{array}\right],\\ {\displaystyle }& {\displaystyle }{G}_s\left(U_s\right)=\left[\begin{array}{c}\mathrm{0}\\ g(\mathrm{1}-r)z_s\end{array}\right],\\ {\displaystyle }& {\displaystyle }{S}_s\left(U_s\right)=\left[\begin{array}{c}\mathrm{0}\\ -r{n}_a\left(u_s-u\right)+{\mathit{\tau }}_P\end{array}\right].\end{array}$$

Analogously, the multilayer non-hydrostatic shallow-water system Eq. (5) can also be expressed in a similar way:

$$\begin{array}{}\text{(7)}& \left\{\begin{array}{l}{\partial }_t{U}_f+{\partial }_x{F}_f\left(U_f\right)+B_f\left(U_f\right){\partial }_x{U}_f\\ =G_f\left(U\right){\partial }_x\left(H-z_s\right)+{\mathcal{T}}_{\mathrm{NH}}-S_f\left(U_f\right),\\ B\left(U_f,{\left(U_f\right)}_x,H,H_x,z_s,{\left(z_s\right)}_x\right)=\mathrm{0},\end{array}\right.\end{array}$$

where

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{U}_f=\left[\begin{array}{c}h\\ h{u}_{\mathrm{1}}\\ \mathrm{⋮}\\ h{u}_L\\ h{w}_{\mathrm{1}}\\ \mathrm{⋮}\\ h{w}_L\end{array}\right],F_f\left(U_f\right)=\left[\begin{array}{c}h\stackrel{\mathrm{‾}}{u}\\ h{u}_{\mathrm{1}}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\\ \mathrm{⋮}\\ h{u}_L^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}g{h}^{\mathrm{2}}\\ h{u}_{\mathrm{1}}{w}_{\mathrm{1}}\\ \mathrm{⋮}\\ h{u}_L{w}_L\end{array}\right],\\ {\displaystyle }& {\displaystyle }{G}_f\left(U_f\right)=\left[\begin{array}{c}\mathrm{0}\\ gh\\ \mathrm{⋮}\\ gh\\ \mathrm{0}\\ \mathrm{⋮}\\ \mathrm{0}\end{array}\right],\end{array}$$

and B_f(U_f)∂_x(U_f) contains the non-conservative products involving the momentum transfer across the interfaces, and, finally, S_f(U_f) represents the friction terms,

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{B}_f\left(U_f\right){\partial }_x\left(U_f\right)=\left[\begin{array}{c}\mathrm{0}\\ u_{\mathrm{3}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{3}/\mathrm{2}}\\ u_{\mathrm{5}/\mathrm{3}}{\mathrm{\Gamma }}_{\mathrm{5}/\mathrm{2}}-u_{\mathrm{3}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{3}/\mathrm{2}}\\ \mathrm{⋮}\\ -u_{L-\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{L-\mathrm{1}/\mathrm{2}}\\ w_{\mathrm{3}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{3}/\mathrm{2}}\\ w_{\mathrm{5}/\mathrm{3}}{\mathrm{\Gamma }}_{\mathrm{5}/\mathrm{2}}-w_{\mathrm{3}/\mathrm{2}}{\mathrm{\Gamma }}_{\mathrm{3}/\mathrm{2}}\\ \mathrm{⋮}\\ -w_{L-\mathrm{1}/\mathrm{2}}{\mathrm{\Gamma }}_{L-\mathrm{1}/\mathrm{2}}\end{array}\right],\\ {\displaystyle }& {\displaystyle }{S}_f\left(U_f\right)=\left[\begin{array}{c}\mathrm{0}\\ {\mathit{\tau }}_b-n_a\left(u_s-u_{\mathrm{1}}\right)\\ \mathrm{0}\\ \mathrm{⋮}\\ \mathrm{0}\end{array}\right].\end{array}$$

The non-hydrostatic corrections in the momentum equations are given by

$$\begin{array}{ll}{\displaystyle }{\mathcal{T}}_{\mathrm{NH}}& {\displaystyle }={\mathcal{T}}_{\mathrm{NH}}\left(h,h_x,H,H_x,z_s,{\left(z_s\right)}_x,p,p_x\right)\\ {\displaystyle }& {\displaystyle }=-\left[\begin{array}{c}\mathrm{0}\\ h\left({\partial }_x{p}_{\mathrm{1}}+{\mathit{\sigma }}_{\mathrm{1}}{\partial }_z{p}_{\mathrm{1}}\right)\\ \mathrm{⋮}\\ h\left({\partial }_x{p}_L+{\mathit{\sigma }}_L{\partial }_z{p}_L\right)\\ h{\partial }_z{p}_{\mathrm{1}}\\ \mathrm{⋮}\\ h{\partial }_z{p}_L\end{array}\right],\end{array}$$

and, finally, the operator related with the incompressibility condition at each layer is given by

$$\begin{array}{ll}{\displaystyle }B& {\displaystyle }\left(U_f,{\left(U_f\right)}_x,H,H_x,z_s,{\left(z_s\right)}_x\right)\\ {\displaystyle }& {\displaystyle }=\left[\begin{array}{c}{\partial }_x{u}_{\mathrm{1}/\mathrm{2}}+{\mathit{\sigma }}_{\mathrm{1}/\mathrm{2}}{\partial }_z{u}_{\mathrm{1}/\mathrm{2}}+{\partial }_z{w}_{\mathrm{1}/\mathrm{2}}\\ \mathrm{⋮}\\ {\partial }_x{u}_{L-\mathrm{1}/\mathrm{2}}+{\mathit{\sigma }}_{L-\mathrm{1}/\mathrm{2}}{\partial }_z{u}_{L-\mathrm{1}/\mathrm{2}}+{\partial }_z{w}_{L-\mathrm{1}/\mathrm{2}}\end{array}\right].\end{array}$$

The discretization of systems (Eqs. 6 and 7) becomes difficult. In the present work, the natural extension of the numerical schemes proposed in Escalante et al. (2018a, b) is considered. These authors propose, describe, and use a splitting technique. Initially, the systems (Eqs. 6 and 7) are expressed as the following non-conservative hyperbolic system:

$$\begin{array}{}\text{(8)}& \left\{\begin{array}{l}{\partial }_t{U}_s+{\partial }_x{F}_s\left(U_s\right)=G_s\left(U_s\right){\partial }_{x}H,\\ {\partial }_t{U}_f+{\partial }_x{F}_f\left(U_f\right)+B_f\left(U_f\right){\partial }_x\left(U_f\right)=G_f\left(U_f\right){\partial }_x\left(H-z_s\right).\end{array}\right.\end{array}$$

Both equations are solved simultaneously using a second-order HLL (Harten–Lax–van Leer), positivity-preserving, and well-balanced, path-conservative finite-volume scheme (see Castro and Fernández-Nieto, 2012) and using the same time step. The synchronization of time steps is performed by taking into account the Courant–Friedrichs–Lewy (CFL) condition of the complete system Eq. (8). A first-order estimation of the maximum of the wave speed for system Eq. (8) is the following:

$${\mathit{\lambda }}_{\max }=\max \left(\left|u_s\right|+\sqrt{\left(g\right(\mathrm{1}-r)z_s},\left|\stackrel{\mathrm{‾}}{u}\right|+\sqrt{gh}\right).$$

Then, the non-hydrostatic pressure corrections $p_{\mathrm{1}/\mathrm{2}}$, ⋯, $p_{L-\mathrm{1}/\mathrm{2}}$ at the vertical interfaces are computed from

$$\left\{\begin{array}{l}{\partial }_t{U}_f={\mathcal{T}}_{\mathrm{NH}}\left(h,h_x,H,H_x,z_s,{\left(z_s\right)}_x,p,p_x\right),\\ B\left(U_f,{\left(U_f\right)}_x,H,H_x,z_s,{\left(z_s\right)}_x\right)=\mathrm{0},\end{array}\right.$$

which requires the discretization of an elliptic operator that is done using standard second-order central finite differences. This results in a linear system than in our case is solved using an iterative scheduled Jacobi method (see Adsuara et al., 2016). Finally, the computed non-hydrostatic corrections are used to update the horizontal and vertical momentum equations at each layer, and, at the same time, the frictions S_s(U_s) and S_f(U_f) are also discretized (see Escalante et al., 2018a, b). For the discretization of the Coulomb friction term, we refer the reader to Fernández-Nieto et al. (2008).

The resulting numerical scheme is well balanced for the water-at-rest stationary solution and is linearly L^∞ stable under the usual CFL condition related to the hydrostatic system. It is also worth mentioning that the numerical scheme is positive preserving and can deal with emerging topographies. Finally, its extension to 2D is straightforward. For dealing with numerical experiments in 2D regions, the computational domain must be decomposed into subsets with a simple geometry, called cells or finite volumes. The 2D numerical algorithm for the hydrodynamic hyperbolic component of the coupled system is well suited to be parallelized and implemented in GPU architectures, as is shown in Castro et al. (2011). Nevertheless, a standard treatment of the elliptic part of the system does not allow the parallelization of the algorithms. The method used here and proposed in Escalante et al. (2018a, b) makes it possible that the second step can also be implemented on GPUs, due to the compactness of the numerical stencil and the easy and massive parallelization of the Jacobi method The above-mentioned parallel GPU and multi-GPU implementation of the complete algorithm results in much shorter computational times.

5 Benchmark problem comparisons

This section presents the numerical results obtained with the Multilayer-HySEA model for the three benchmark problems dealing with granular slides and the comparison with the measured lab data for the generated water waves. In particular, BP4 deals with a 2D submarine granular slide, BP5 with a 2D subaerial slide, and BP6 with a 3D subaerial slide. The description of all these benchmarks can be found in LTMBW (2017) and Kirby et al. (2018). In this paper, all units, unless otherwise indicated, will be expressed in the SI (International System of Units) units.

The model parameters required in each simulation are

$$g,r,n_a,n_m,d_s,{\mathit{\delta }}_i,\mathit{\beta },\mathrm{and}\mathit{\gamma }.$$

The parameters g, r, n_m, and d_s are related to physical settings given in each experiment. β and γ are empirical parameters that were chosen as in the seminal paper of Pouliquen and Forterre (2002).

The friction angles δ₁ and δ₂ are characteristic angles of the material, and δ₃ is related to the behavior of the slide motion when starting from the rest. Thus, the values of these angles strongly depend on the granular material. They were adjusted within a range of feasible values according to the references (Brunet et al., 2017; Mangeney et al., 2007; Pouliquen and Forterre, 2002):

$${\mathit{\delta }}_{\mathrm{1}}\in [\mathrm{1},\mathrm{22}{}^{\circ }],{\mathit{\delta }}_{\mathrm{2}}\in [\mathrm{11},\mathrm{34}{}^{\circ }],{\mathit{\delta }}_{\mathrm{3}}\in [\mathrm{3},\mathrm{23}{}^{\circ }].$$

In the present paper we have used the values

$${\mathit{\delta }}_{\mathrm{1}}=\mathrm{6}{}^{\circ },{\mathit{\delta }}_{\mathrm{2}}\in [\mathrm{17},\mathrm{30}{}^{\circ }],{\mathit{\delta }}_{\mathrm{3}}=\mathrm{12}{}^{\circ }$$

for the three benchmark problems, which is consistent with the values found in the literature. As noted in Mangeney et al. (2007), in general for real problems involving complex rheologies, smaller values of these parameters δ_i should be employed.

With regard to the sensitivity of the model to parameter variation, an appropriate sensitivity analysis can be performed, as is done in González-Vida et al. (2019). However, the aim of the present work was to prove if the non-hydrostatic model coupled with the granular model was able to accurately reproduce the three benchmarks considered.

Regarding the parameter denoting the buoyancy effect, for field cases, r=0.5 is usually taken, and then the parameter is eventually adjusted based on available field data. In general, the complexity of the rheology introduces a difficulty that is always present on the modeling as well as on the adjustment of the parameters. Moreover, the more sophisticated the model is (considering, for example, the rheology of the material), more input data will be required.

We would like to stress the simplicity of the slide model used here as a great advantage regarding parameter setup. Although the end user has to adjust some input parameters of the model, within a range of acceptable value, the simplicity of the proposed numerical model makes this task remain simple – not representing an obstacle to run the model. On the other hand, the efficient GPU implementation of the model allows for performing uncertainty quantification (see Sánchez-Linares et al., 2016) on a few parameters and investigating the sensitivity to them, varying on small ranges (as in González-Vida et al., 2019). This will be the aim of future works. When field or experimental observations are available, a different approach is proposed in Ferreiro-Ferreiro et al. (2020), where an automatic data assimilation strategy for a similar landslide non-hydrostatic model is proposed. The same strategy can be adapted for the model used here.

5.1 Benchmark problem 4: two-dimensional submarine granular slide

The first proposed benchmark problem for granular slides, BP4 in the list, aims to reproduce the generation of tsunamis by submarine granular slides modeled in the laboratory experiment by means of glass beads. The corresponding 2D laboratory experiments were performed at the École Centrale de Marseille (see Grilli et al., 2017, for a description of the experiment). A set of 58 (29 with their corresponding replicate) experiments were performed at the IRPHE (Institut de Recherche sur les Phénomènes Hors Equilibre) precision tank. The experiments were performed using a triangular submarine cavity filled with glass beads that were released by lifting a sluice gate and then moving down a plane slope, all underwater. Figure 2 shows a schematic picture of the experiment setup.

Figure 2 BP4 sketch showing the longitudinal cross section of the IRPHE's precision tank. The figure shows the location of the plane slope, the sluice gate, and the four gages (WG1–WG4).

The one-dimensional domain [0, 6] is discretized with Δx=0.005 m, and wall boundary conditions are imposed. The simulated time is 10 s. The CFL number was set to 0.5, and model parameters take the following values:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }g=\mathrm{9.81},r=\mathrm{0.78},n_a=\mathrm{0.2},n_m={\mathrm{10}}^{-\mathrm{3}},\\ {\displaystyle }& {\displaystyle }{d}_s=\mathrm{7}\times {\mathrm{10}}^{-\mathrm{3}},{\mathit{\delta }}_{\mathrm{1}}=\mathrm{6}{}^{\circ },{\mathit{\delta }}_{\mathrm{2}}=\mathrm{17}{}^{\circ },{\mathit{\delta }}_{\mathrm{3}}=\mathrm{12}{}^{\circ },\\ {\displaystyle }& {\displaystyle }\mathit{\beta }=\mathrm{0.136},\mathit{\gamma }={\mathrm{10}}^{-\mathrm{3}}.\end{array}$$

Figure 3 depicts the modeled times series for the water height at the four wave gages and compares them with the lab-measured data. Note that the computed free surface matches well with the laboratory data for gauges WG2–WG4, both in amplitude and frequency. For gauge WG1, some mismatch is observed in amplitude, which could be explained by the simplicity of the landslide model and the absence of turbulent effects in the model.

Figure 3 Numerical time series for the simulated water surface (in blue) compared with lab-measured data (red) at wave gauges (a) WG1, (b) WG2, (c) WG3, and (d) WG4.

Figure 4 Modeled location of the granular material and water free surface elevation at times t=0, 0.3, 0.6, and 0.9 s.

Figure 4 shows the location and evolution of the granular material and water free surface at several times during the numerical simulation. In Grilli et al. (2017) some snapshots of the landslide evolution are shown at different time steps. Compared to those snapshots, it can be observed that the location of the landslide front is well captured by the model, but there is some mismatch in the landslide shape at the front, mainly due to the simplicity of the landslide model considered here. In particular, we consider that density remains constant in the landslide layer during the simulation, which is not true due to the water entrainment.

In the numerical experiments presented in this section, the number of layers was set up to five. Similar results were obtained with a lower number of layers (four or three) but slightly closer to measured data when considering five layers. This justifies our choice in the present test problem. A larger number of layers does not further improve the numerical results. This may indicate that to get better numerical results, it is no longer a question related to the dispersive properties of the model (which improve with the number of layers) but is more likely due to some missing physics.

5.2 Benchmark problem 5: two-dimensional subaerial granular slide

This benchmark is based on a series of 2D laboratory experiments performed by Viroulet et al. (2014) in a small tank at the École Centrale de Marseille, France. The simplified picture of the setup for these experiments can be found in Fig. 5.

Figure 5 BP5 sketch of the setup for the laboratory experiments.

The granular material was confined in triangular subaerial cavities and composed of dry glass beads of diameter d_s (which was varied) and density ρ_s=2500 km m⁻³. This was located on a plane with a 45^∘ slope just on top of the water surface. Then the slide was released by lifting a sluice gate and coming in contact with water right away. The experimental setup used by Viroulet et al. (2014) consisted of a wave tank 2.2 m long, 0.4 m high, and 0.2 m wide.

The granular material is initially retained by a vertical gate on the dry slope. The gate is suddenly lowered, and in the numerical experiments, it should be assumed that the gate release velocity is large enough to neglect the time it takes the gate to withdraw. The front face of the granular slide touches the water surface at t=0. The initial slide shape has a triangular cross section over the width of the tank, with down-tank length L and front face height B=L as the slope angle is 45^∘.

For the present benchmark, two cases are considered. Case 1 is defined by the following setup: d_s=1.5 mm, H=14.8 cm, and L=11 cm. Case 2 is given by d_s=10 mm, H=15 cm, and L=13.5 cm. The benchmark problem proposed consists of simulating the free-surface elevation evolution at the four gauges WG1 to WG4, where measured data are provided, for the two test cases described above.

The same model configuration as in the previous benchmark problem is used here. The vertical structure is reproduced using three layers in the present case. The one-dimensional domain is given by the interval [0, 2.2], and it is discretized using a step Δx=0.003 m. As boundary conditions, rigid walls were imposed. The simulation time is 2.5 s. The CFL number is set to 0.9, and the model parameters take the following values:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }g=\mathrm{9.81},r=\mathrm{0.6},n_a={\mathrm{10}}^{-\mathrm{2}},n_m=\mathrm{9}\times {\mathrm{10}}^{-\mathrm{2}},\\ {\displaystyle }& {\displaystyle }{\mathit{\delta }}_{\mathrm{1}}=\mathrm{6}{}^{\circ },{\mathit{\delta }}_{\mathrm{2}}=\mathrm{26}{}^{\circ },{\mathit{\delta }}_{\mathrm{3}}=\mathrm{12}{}^{\circ },\mathit{\beta }=\mathrm{0.136},\mathit{\gamma }={\mathrm{10}}^{-\mathrm{3}}.\end{array}$$

Finally d_s was set to $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{10}\times {\mathrm{10}}^{-\mathrm{3}}$ depending on the test case.

Figure 6 shows the comparison for Case 1. In this case, the numerical results show a very good agreement when compared with lab-measured data, and, in particular, the two leading waves are very well captured. Figure 7 shows the comparison for Case 2. In this case, the agreement is good, but larger differences between model and lab measurements can be observed.

Figure 6 Numerical time series for the simulated water surface (in blue) compared with lab-measured data (red). Case 1 at gauges (a) G1, (b) G2, (c) G3, and (d) G4.

Figure 7 Numerical time series for the simulated water surface (in blue) compared with lab-measured data (red). Case 2 at gauges (a) G1, (b) G2, (c) G3, and (d) G4.

Two things can be concluded from the observation of Figs. 6 and 7: (1) a much better agreement is obtained for Case 1 than for Case 2 and (2) the agreement is better for gauges located farther from the slide compared with closer to the slide gauges. Although paradoxical, this second differential behavior among gauges can be explained as a consequence of the hydrodynamic component being much better resolved and simulated than the morphodynamic component (the movement of the slide material), which is obviously much more difficult to reproduce. But, at the same time, this implies a correct transfer of energy at the initial stages of the interaction between the slide and fluid.

Finally, Fig. 8 shows the location of the granular material and the free-surface elevation at several times for the numerical simulation of Case 1. In Viroulet et al. (2014) some snapshots of the landslide evolution are shown at different time steps that can be compared with Fig. 8. As for the benchmark problem 4, it can be seen that the location of the landslide front is well captured, but there is some mismatch in the landslide shape at the front, mainly due to the simplicity of the landslide model considered here.

Figure 8 Modeled water free-surface elevation and granular slide location at times t=0, 0.2, 0.4, and 0.8 s for Case 1.

5.3 Benchmark problem 6: three-dimensional subaerial granular slide

This benchmark problem is based on the 3D laboratory experiment of Mohammed and Fritz (2012) and Mohammed (2010). Benchmark 6 simulates the rapid entry of a granular slide into a 3D water body. The landslide tsunami experiments were conducted at Oregon State University in Corvallis. The landslides are deployed off a plane with a 27.1^∘ slope, as shown in Fig. 9.

Figure 9 Schematic picture of the computational domain. Plan view in panel (a). Cross section at y=0 m in panel (b). The red dots represent the distribution of the wave gauge positions in the laboratory setup.

The landslide material is deployed using a box measuring 2.1 m × 1.2 m × 0.3 m, with a volume of 0.756 m³, and weighting approximately 1360 kg. The case selected by the NTHMP as a benchmarking test is the one with a still-water depth of H=0.6 m (see Fig. 9). The computational domain is the rectangle defined by [0, 48] × [−14, 14], and it is discretized with $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{0.06}$ m. At the boundaries, wall boundary conditions were imposed. The simulation time is 20 s and we set the CFL=0.5. According to Mohammed and Fritz (2012) and Mohammed (2010), the three-dimensional granular landslide parameters were set to

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }g=\mathrm{9.81},r=\mathrm{0.55},n_a=\mathrm{4},n_m=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{2}},\\ {\displaystyle }& {\displaystyle }{d}_s=\mathrm{13.7}\times {\mathrm{10}}^{-\mathrm{3}},{\mathit{\delta }}_{\mathrm{1}}=\mathrm{6}{}^{\circ },{\mathit{\delta }}_{\mathrm{2}}=\mathrm{30}{}^{\circ },{\mathit{\delta }}_{\mathrm{3}}=\mathrm{12}{}^{\circ },\\ {\displaystyle }& {\displaystyle }\mathit{\beta }=\mathrm{0.136},\mathit{\gamma }={\mathrm{10}}^{-\mathrm{3}}.\end{array}$$

The vertical structure of the fluid layer is modeled using three layers. Similar results were obtained with two layers.

Initially, the slide box is driven using four pneumatic pistons. Here we provide comparisons for the case where the pressure for the pneumatic pistons generating the slide is P=0.4 MPa (P=58 PSI). In Mohammed (2010), it is shown that for this test case, the landslide box reached a velocity of $v_b=\mathrm{2.3}\cdot \sqrt{g\cdot \mathrm{0.6}}=\mathrm{5.58}$ m s⁻¹. Thus, the initial condition for the water velocities is set to zero,

$$u_i=\mathrm{0},i=\mathrm{1},\mathrm{2},\mathrm{\dots },L,$$

and for the landslide velocity it is set to the above-mentioned constant value,

$$u_s=\mathrm{5.58},\mathrm{wherever}{z}_s>\mathrm{0},$$

for the x component. The y component of the landslide velocity was initially set to zero.

Figure 10 BP6 cross section at y=0 m. Panel (a) shows the location and velocity of the granular slide and the generated wave at time t=0.4 s from the triggering and (b) the final deposit location (at t=20 s).

Figure 11 Simulated (solid blue lines) time series compared with measured (dashed red lines) free-surface waves for the nine wave gauges considered at (R, θ^∘) positions.

The benchmark problem proposed consists of simulating the free-surface elevation at some wave gauges. In the present study, we include the comparison for the nine wave gauges displayed in Fig. 9 as red dots. A total number of 21 wave gauges composed the whole set of data plus five run-up gauges. The wave gauges in coordinates (R, θ^∘) are given more precisely in Table 1.

Table 1 Location of the nine waves gauges referenced to the toe's slope.

Before comparing time series, we first check the simulated landslide velocity at impact with the measured one. The slide impact velocity measured in the lab experiment is 5.72 m s⁻¹ at time t=0.44 s. The numerically computed slide impact velocity is slightly underestimated with a value of 5.365 m s⁻¹ at time t=0.4 s as it can be seen in the upper panel of Fig. 10. The final simulated granular deposit is located partially on the final part of the sloping floor and partially at the flat bottom closer to the point of change of slope as is shown in the lower panel of Fig. 10. This can be compared with the actual final location of the granular material in the experimental setup. The simulated deposits, being thinner, extend further. This is probably due to the fact that we are neglecting the friction that it is produced by the change in the slope at the transition area. In Ma et al. (2015) a similar result and discussion can be found.

Figure 11 presents the comparisons between the simulated and the measured waves at the nine gauges we have retained. Model results are in good agreement with measured time. Despite this, wave heights are overestimated at some stations, specially those closer to the shoreline (for example, the station with R=3.9 and θ=30^∘). This effect has been also observed and discussed in Ma et al. (2015). At some of the time series, it can be observed that the small free-surface oscillations at the final part of the time series are not well captured by the model. This is partially due to the relatively coarse horizontal grids used in the simulation. This same behavior can be also observed in Fig. 12 in this case for the comparisons between the simulated and measured run-up values at some measure locations situated at the shoreline (as for x=7.53).

Figure 12 Time series comparing numerical run-up (solid blue) at the four run-up gauges with the measured (dashed red) data at (x, y) positions.

Table 2 Wall-clock times in seconds for the hydrostatic shallow-water Savage–Hutter system (SWE-SH) and the non-hydrostatic GPU implementations. The ratios are with respect to the SWE-SH model implementation.

Table 2 shows the wall-clock times on an NVIDIA Tesla P100 GPU. It can be observed that including non-hydrostatic terms in the SWE-SH system results in an increase in the computational time of 2.9 times. If a richer vertical structure is considered, then larger computational times are required. As in the examples for the two- and three-layer systems, there is a 3.48 and 4.66 times increase in the computational effort.

6 Concluding remarks

Numerical models need to be validated prior to their use as predictive tools. This requirement becomes even more necessary when these models are going to be used for risk assessment in natural hazards where human lives are involved. The current work aims at benchmarking the novel Multilayer-HySEA model for landslide-generated tsunamis produced by granular slides, in order to provide the tsunami community with a robust, efficient, and reliable tool for landslide tsunami hazard assessment in the future.

The Multilayer-HySEA code implements a two-phase model to describe the interaction between landslides (aerial or subaerial) and water body. The upper phase describes the hydrodynamic component. This is done using a stratified vertical structure that includes non-hydrostatic terms in order to include dispersive effects in the propagation of simulated waves. The motion of the landslide is taken into account by the lower phase, consisting of a Savage–Hutter model. To reproduce these flows, the friction model given in Pouliquen and Forterre (2002) is considered here. The hydrodynamic and morphodynamic models are weakly coupled through the boundary condition at their interface.

The implemented numerical algorithm combines a finite-volume path-conservative scheme for the underlying hyperbolic system and finite differences for the discretization of the non-hydrostatic terms. The numerical model is implemented to be run in GPU architectures. The two-layer non-hydrostatic code coupled with the Savage–Hutter used here has been shown to run at very efficient computational times. To assess this, we compare with respect to the one-layer SWE/Savage–Hutter GPU code. For the numerical simulations performed here, the execution times for the non-hydrostatic model are always below 4.66 times the times for the SWE model for a number of layers up to three. We can conclude that the numerical scheme presented here is very robust, is extremely efficient, and can model dispersive effects generated by submarine/subaerial landslides at a low computational cost considering that dispersive effects and a vertical multilayer structure are included in the model. Model results show a good agreement with the experimental data for the three benchmark problems considered, in particular for BP5, but this also occurs for the other two benchmark problems. In general, a better agreement for the hydrodynamic component, compared with its morphodynamic counterpart, is shown, which is more challenging to reproduce.