Multilayer-HySEA model validation for landslide 1 generated tsunamis. Part I Rigid slides

6 This paper is devoted to benchmarking the Multilayer-HySEA model using lab- 7 oratory experimental data for landslide-generated tsunamis. This article deals 8 with rigid slides and the second part, in a companion paper, addresses granular 9 slides. The US National Tsunami Hazard and Mitigation Program (NTHMP) 10 has proposed the experimental data used and established for the NTHMP Land- 11 slide Benchmark Workshop, held in January 2017 at Galveston (Texas). The 12 ﬁrst three benchmark problems proposed in this workshop deal with rigid slides. 13 Rigid slides must be simulated as a moving bottom topography and, therefore, 14 they must be modelled as a prescribed boundary condition. These three bench- 15 marks are used here to validate the Multilayer-HySEA model. This new HySEA 16 model consists of an eﬃcient hybrid ﬁnite-volume/ﬁnite-diﬀerence implementa- 17 tion on GPU architectures of a non-hydrostatic multilayer model. A brief de- 18 scription of model equations, dispersive properties, and the numerical scheme 19 is included. The benchmarks are described and the numerical results compared 20 against the lab-measured data for each of them. The speciﬁc aim is to validate 21 this new code for tsunamis generated by rigid slides. Nevertheless, the over- 22 all objective of the current benchmarking eﬀort is to produce a ready-to-use 23 numerical tool for real-world landslide generated tsunami hazard assessment. 24 This tool has already been used to reproduce the Port Valdez Alaska 1964 and 25 Stromboli Italy 2002 events. 26

(3) Total depth, h, is split along the vertical axis into L ≥ 1 layers and ∆s = 1/L 143 (see Figure 1). The variables u α and w α are the depth-averaged velocities in 144 the x and z directions, respectively, t is time and g is gravitational acceleration.
The non-hydrostatic pressure at the interface z α+1/2 is denoted by p α+1/2 . The 146 water surface elevation measured from the still-water level is η = h − H, where 147 H is the water depth when the water is at rest. Finally, τ is a friction law term, 148 and the terms Γ α+1/2 account for the mass transfer across interfaces and are  To obtain such properties, the system (1) is linearised around the water at 164 rest steady-state solution. After that, a Stokes-type Fourier analysis is carried 165 out looking for first-order planar wave solutions. This method constitutes a 166 standard procedure to study systems that model dispersive water waves (see 167 velocities as well as the linear shoaling gradient are, respectively, defined as: where ω denotes the angular frequency, k the local wave-number and H the 171 typical depth.

172
The measured quantities C, G and γ are solely functions of the local wave- one can also obtain the corresponding linear dispersion relations that state the 176 linear theory for the considered quantities (see Schäffer and Madsen (1995) for 177 the Airy reference formulae). For example, the expression for the phase velocity 178 from the Airy's theory is The expressions of the phase velocity for the system (1) are given in Table 1 180 for the non-linear hydrostatic shallow water system (SWE) and the Multilayer-

181
HySEA (non-hydrostatic) system with j ≥ 1 layers (NH-jL to 3 (still a small value) the system (1) is in an excellent agreement with the 203 Airy theory for kH up to 15. For the phase celerity, the percentage error is less 204 than 0.62%, and for the group velocity is less than 1% for kH smaller than 10 205 (see Figure 2). Linear shoaling is also well reproduced in this same range.

206
The Multilayer-HySEA model presents enhanced dispersive properties. In 207 order to have similar dispersive results as the ones obtained here using a three-  reads for each α ∈ {1, 2, . . . , L} as: where ς = z α−1/2 z α−1/2 ∂ z ν is the eddy viscosity. In this work, as in Escalante et al.
where B is an empirical parameter related to a breaking criteria to switch on where denote the flow speed at the onset and termination of the wave-breaking process which gives the exact value of u α and w α for h ≥ , and gives a smooth 254 transition of u α and w α to zero when h tends to zero without truncation.

255
In this work we set = 10 −3 for the numerical tests.

Numerical Solution Method
We describe now the discretization of system (1) that follows the natural ex-

267
First, we shall solve the non-conservative hyperbolic underlying system (1) 268 given by the compact equation where the following compact notation has been used: and B SW is a matrix such B SW ∂ x U contains the non-conservative products 271 related to the mass transfer across interfaces appearing at the momentum equa-272 tions.

273
Then, in a second step, non-hydrostatic terms given by the pressure vector 274 correction term as well as the divergence constraints at each layer will be taken into account.
Non-hydrostatic terms will be approximated by second order compact finite-285 differences.

286
Let us detail de time stepping procedure followed. Assume given time steps ∆t n , and denote t n = k≤n ∆t k . To obtain second order accuracy in time, the two-stage second-order TVD Runge-Kutta scheme is adopted. At the kth stage, k ∈ {1, 2}, the two-step projection-correction method is given by Observe that, equations (20b-20c) requires, at each stage of the calculation 292 respectively, to solve a Poisson-like equation for each one of the variables con-  Regarding the wave breaking model, the breaking mechanism described in Sub-321 section 3.2 was implemented, adopting B 1 = 0.6, B 2 = 0.15 for all the bench-322 mark problems, and K = 10 for the third benchmark and K = 2 for the rest.

323
The description of all these benchmarks can be found at LTMBW (2017) and  Table 3-).
where S 0 = u 2 t /a 0 = 2.110 m, t 0 = u t /a 0 = 1.677 s, a 0 = 0.75 m/s 2 and u t = 346 1.258 m/s is the terminal velocity. Figure 4 shows the prescribed acceleration, The benchmark here consists of using the above information on slide shape,        , y) in mm, as shown in Figure 6.
presented. An excellent agreement can be observed between these time series.

407
The comparisons for the second required case (d = 120 mm) in the 3 gauges 408 with data provided (gauge g 3 was not available) are shown in Figure 8. Good 409 agreement can also be observed in this case. Finally, Figure 9 shows the com-       of kH (see Table 1 and Figure 2), this explains the aforementioned excellent 435 agreement between the computed time series and the measured lab data.

436
Finally, although the phase velocity for the one-layer system shows an error 437 bounded by only 3.02 % for kH ∈ [0, 5] (see Table 1), it can be seen in Figure 10 438 that the one-layer non-hydrostatic pressure system cannot represent the waves 439 correctly. In contrast, the one-layer system tends to amplify waves. This be-440 haviour can be explained by observing the shoaling gradient for this model (see 441 Figure 2). The shoaling gradient verifies the ODE where A denotes the amplitude. Then, it can be stated by inspecting the so-   g 4 2.815 3.093 3.218 3.512 3.14 3.218 3.512 Table 8: Computed kH values from the measured wave period (see Table 7) and the Airy dispersion relation 2π/T = gk tanh(kH) for test cases d = 61, 80, 100, 120, 140, 149 and 189.

463
The block movement was provided by means of a polynomial fitting to mea-464 sured data, giving the horizontal distance as: with β = arctan(1/2) and x (0,t=0) = −2∆. The polynomial coefficients for the 466 two cases proposed are given in Table 9.       cost.

518
Regarding model results, they show a good agreement with the experimental 519 data for the three benchmark problems studied in the present work. In partic-520 ular, for BP2, but this also occurs for the other two benchmark problems, we 521 have shown that a one layer, hydrostatic or non-hydrostatic, model is not able 522 to reproduce the complexity in the observed lab data considered in the pro-   The authors declare that they have no conflict of interest.