This paper is devoted to benchmarking the Multilayer-HySEA model using laboratory experimental data for landslide-generated tsunamis. This article deals with rigid slides, and the second part, in a companion paper, addresses granular slides. The US National Tsunami Hazard and Mitigation Program (NTHMP) has proposed the experimental data used and established for the NTHMP Landslide Benchmark Workshop, held in January 2017 at Galveston (Texas). The first three benchmark problems proposed in this workshop deal with rigid slides. Rigid slides must be simulated as a moving bottom topography, and, therefore, they must be modeled as a prescribed boundary condition. These three benchmarks are used here to validate the Multilayer-HySEA model. This new HySEA model consists of an efficient hybrid finite-volume–finite-difference implementation on GPU architectures of a non-hydrostatic multilayer model. A brief description of model equations, dispersive properties, and the numerical scheme is included. The benchmarks are described and the numerical results compared against the lab-measured data for each of them. The specific aim is to validate this new code for tsunamis generated by rigid slides. Nevertheless, the overall objective of the current benchmarking effort is to produce a ready-to-use numerical tool for real-world landslide-generated tsunami hazard assessment. This tool has already been used to reproduce the Port Valdez, Alaska, 1964 and Stromboli, Italy, 2002 events.

Model development and benchmarking for earthquake-induced tsunamis is a task addressed
in the past and to which much effort and time has been dedicated.
In particular, just to mention a couple of NTHMP efforts, the 2011 Galveston benchmarking workshop

In its 2019 strategic plan,
the NTHMP required that all numerical tsunami inundation models to be used in hazard assessment
studies in the United States be verified as accurate and consistent through a model benchmarking process.
This mandate was fulfilled in 2011 but only for seismic tsunami sources and to a limited extent for idealized
solid underwater landslides.
However, recent work by various NTHMP states has shown that landslide tsunami hazard may in fact be greater than
seismically induced hazard and may be also the dominant risk along significant parts of the US coastline

As a result of this demonstrated gap, a set of candidate benchmarks was proposed to perform the required validation process.
The selected benchmarks are based on a subset of available laboratory data sets for solid slide experiments and deformable slide experiments and include both submarine and subaerial slides.
In order to complete this list of laboratory data, a benchmark based on a historic field event (Valdez, AK, 1964)
was also chosen.
The EDANYA group (

Twenty years ago, at the beginning of the century, the challenge of solid block landslide modeling was taken by a number of researchers

The HySEA (Hyperbolic Systems and Efficient Algorithms) software consists
of a family of geophysical codes based on either single-layer, two-layer
stratified systems, or multilayer shallow-water models.
HySEA codes

The waves generated in the laboratory tests proposed in the NTHMP selected benchmarks are high frequency
and dispersive, and the generated flows have a complex vertical structure.
Therefore, the numerical model used must be able to reproduce such effects.
This makes the two-layer Landslide-HySEA model unsuitable for reproducing these experimental results
as non-hydrostatic effects play an important role and a richer vertical structure
is required.
To address these requirements, the Multilayer-HySEA model was very recently implemented, considering a
stratified structure in the simulated fluid and including non-hydrostatic terms. A multilayer model is able to better approximate the vertical structure
of a complex flow than a standard one-layer depth-averaged model.
In particular, by increasing the number of layers, the linear dispersion relation of the model converges towards the exact dispersion relation from the Stokes linear theory (see

The Multilayer-HySEA model implements one of the multilayer non-hydrostatic models
of the family introduced and described in

An alternative deduction for this system is performed in

Schematic diagram describing the multilayer system.

Total depth,

In order to close the system of equations, the following boundary conditions are considered:

Note that the motion of the bottom surface can be taken into account as a boundary condition, imposing

Phase velocity expressions and maximum of the relative error

Some dispersive properties of the system (Eq.

To obtain such properties, system Eq. (

The measured quantities

The expressions of the phase velocity for the system (Eq.

The main goal when deriving dispersive shallow-water systems is to get the most accurate
dispersive relations as possible, compared with the Airy wave theory, without highly increasing
the complexity of the system.
See

The Multilayer-HySEA model presents enhanced dispersive properties.
In order to have similar dispersive results as the ones obtained here
using a three-layer system, at least five layers are required for other similar multilayer models
as the one presented in

Relative error for the phase velocities

In shallow areas the breaking of waves can be observed near the coast. As pointed out in

Finally, a natural and simple extension of the criterion proposed by

For the computation of variables in areas of small water depth, a wet–dry treatment adapting the ideas described in

The hydrostatic pressure terms

To overcome the difficulties due to large round-off
errors in computing the velocities

We describe now the discretization of system Eq. (

First, we shall solve the non-conservative hyperbolic underlying shallow-water (SW) system (Eq.

where the following compact
notation has been used:

BP1. Sketch of main parameters and variables for wave generation by 2D rigid slide (modified from Grilli and Watts, 2005).

BP1. Prescribed acceleration, velocity, and displacement of the solid slide.

Then, in a second step, non-hydrostatic terms given by
the pressure vector correction term

System Eq. (

Let us detail the time stepping procedure implemented.
Assume given time steps

With respect to the breaking mechanism introduced in Sect.

In this section, the numerical results obtained with the Multilayer-HySEA model and the comparison with the
measured lab data for waves generated by the movement of a rigid bottom surface or of a solid block are presented.
In particular, BP1 deals with a 2D submarine solid slide, BP2 with a 3D submarine slide, and, finally,
BP3 consists of two 3D slides, one partially submerged and a second one representing a completely submarine slide.
In all these cases, a moving bottom condition has been used to model the solid block movement.
Regarding the wave-breaking model, the breaking mechanism described in Sect.

This benchmark problem is based on the 2D laboratory experiments of

BP1. Filtered data (in red) and numerically simulated (in blue) time series
at wave gauges

BP2. Sketch of the plan view (case 61 mm) (from

Values for variables defining setup configuration.

Gauge positions in dimensional and non-dimensional units.

In this benchmark, two items remained not completely determined in the original description provided:
the first one is related with the initialization of the numerical experiment, and
the second one is related with how and where the solid moving block must stop.
Other small issues related to the description of the benchmark were put forward in

The motion of the rigid slide was prescribed as a function of time as

The benchmark here consists of using the above information on slide shape, submergence, and kinematics, together with reproducing the experimental setup to simulate surface elevations measured at the four wave gauges (average of two replicates of experiments provided).

Then, in order to reproduce the lab experiment, the interval

In Fig.

Measured and curve-fitted slide and wave parameters for the seven experiments performed by

This second benchmark consists of a 3D extension of BP1. The longitudinal sketch of the experiment is the same as in Fig.

Test case

Test case

For the numerical simulations, the two-dimensional computational domain

The benchmark test proposed consists in reproducing the slide shape and complete experimental setup in
and using the information about submergence and kinematics to replicate numerically
Enet and Grilli's experiments for

Wave gauge locations

Execution times in seconds for SWE and non-hydrostatic GPU implementations. Ratios compared with SWE.

In Fig.

Comparison of data and numerical solution time series
at wave gauges (dashed) for the cases

In Table

Figure

Moreover, Table

Since multilayer models have good dispersion relation errors within this range of

Test case

Definition sketch for BP3 laboratory experiments –
here for a submerged (

Finally, although the phase velocity for the one-layer system shows an error bounded
by only 3.02

Measured wave period

Computed

This benchmark problem is based on the 3D laboratory experiment of

A plane slope

The block movement was provided by means of a polynomial fitting to
measured data, giving the horizontal distance as

Polynomial coefficients defining slide motion.

For each case, measured free-surface elevations are given for two wave gauges
placed at

The two-dimensional computational domain

The numerical results obtained for the subaerial test case are presented in
Figs.

Subaerial test case. Lab-measured water height (red) and numerical time series (blue)
at wave gauges

Subaerial test case. Lab-measured run-up (red) and numerical time series
(blue) at run-up gauges

Submerged test case. Lab-measured water height (red) and numerical time series (blue)
at wave gauges

Submerged test case. Lab-measured run-up (red) and numerical time series
(blue) at run-up gauges

Validation of numerical models is a first unavoidable step before their use as predictive tools.
This requirement is even more necessary when the developed models are going to be used for risk assessment in natural events where human lives are involved.
The present work is the first step in this task for the Multilayer-HySEA model, a novel dispersive multilayer model of the HySEA suite developed at the University of Malaga.
This model considers a stratified vertical structure and includes non-hydrostatic terms;
this is done in order to include the dispersive effects in the propagation of the waves
in a homogeneous, inviscid, and incompressible fluid.
The numerical scheme implemented combines a highly robust and efficient finite-volume path-conservative
scheme for the underlying hyperbolic system and finite differences for the discretization
of the non-hydrostatic terms.
In order to increase numerical efficiency, the numerical model is implemented to run in GPU architectures –
in particular in NVIDIA graphics cards and using CUDA language.
In the case of the traditional SW non-dispersive model, this kind of implementation produces an extremely
efficient and fast code

Regarding model results, they show a good agreement with the experimental data for the three benchmark
problems studied in the present work.
In particular, for BP2, but this also occurs for the other two benchmark problems, we have shown that
a one-layer, hydrostatic or non-hydrostatic, model is not able to reproduce the complexity
in the observed lab data considered in the proposed benchmarks.
The waves to be modeled in the test cases proposed here are high in frequency and dispersive.
Hence, it is at least necessary that a two-layer structure and non-hydrostatic terms in the model be used
in order to capture the dynamics of the generated waves.
As noted in

The numerical code used to perform the numerical simulations in this
paper is available at the HySEA codes web page at

All the data used in the present work and necessary to reproduce the setup of the numerical experiments as well as the laboratory-measured data to compare with can be downloaded from

JM led the HySEA codes benchmarking effort undertaken by the EDANYA group, wrote most of the paper, reviewed and edited it, and assisted in the numerical experiments and in their setup. CE implemented the numerical code and performed all the numerical experiments; he also contributed to writing the manuscript. JM and CE did all the figures. MC significantly contributed to the design and implementation of the numerical code.

The authors declare that they have no conflict of interest.

The authors are indebted to Diego Arcas (PMEL/NOAA) and Victor Huérfano (PRSN) for supporting our participation in the 2017 Galveston workshop and to the MMS of the NTHMP for kindly inviting us to that event.

This research has been supported by the Spanish Government–FEDER (MEGAFLOW project, grant no. RTI2018-096064-B-C21) and the Junta de Andalucía–FEDER (grant no. UMA18-Federja-161).

This paper was edited by Maria Ana Baptista and reviewed by two anonymous referees.