A Numerical Landslide-Tsunami Hazard Assessment Technique Applied on Hypothetical Scenarios at Es Vedrà, Offshore Ibiza

This study presents a numerical landslide-tsunami hazard assessment technique for applications in reservoirs, lakes, fjords, and the sea. This technique is illustrated with hypothetical scenarios at Es Vedrà, offshore Ibiza, although currently no evidence suggests that this island may become unstable. The two selected scenarios include two particularly vulnerable locations, namely: (i) Cala d’Hort on Ibiza (3 km away from Es Vedrà) and (ii) Marina de Formentera (23 km away from Es Vedrà). The violent wave generation process is modelled with the meshless Lagrangian method smoothed particle hydrodynamics. Further offshore, the simulations are continued with the less computational expensive code SWASH (Simulating WAves till SHore), which is based on the non-hydrostatic non-linear shallow water equations that are capable of considering bottom friction and frequency dispersion. The up to 133-m high tsunamis decay relatively fast with distance from Es Vedrà; the wave height 5 m offshore Cala d’Hort is 14.2 m, reaching a maximum run-up height of over 21.5 m, whilst the offshore wave height (2.7 m) and maximum inundation depth at Marina de Formentera (1.2 m) are significantly smaller. This study illustrates that landslide-tsunami hazard assessment can nowadays readily be conducted under consideration of site-specific details such as the bathymetry and topography, and intends to support future investigations of real landslide-tsunami cases.


Introduction
Landslides are a main source of impulse waves in reservoirs [1,2], lakes [3,4], fjords [5], and the sea [6]. Such impulse waves are better known as landslide tsunamis if they occur in the open sea, and for simplicity, this term is adopted herein irrespectively of the type of water body in question. Landslide tsunamis involve several complex processes such as slide propagation and impact into a water body, tsunami generation, propagation, and run-up, which often results in casualties and destroyed infrastructure. Landslide tsunamis from nearshore to offshore involve multiple time and length scales such that more than one numerical code is required in order to reliably deal with all of these scales, e.g., DualSPHysics for the nearshore [7,8] combined with SWASH (Simulating WAves till SHore) [9,10] for the offshore region.
Landslide tsunamis occurred for example in Lake Askja, on Iceland, on 21 July 2014. An approximately 20 million m 3 large landslide initiated a 50-m large tsunami resulting in a run-up SWASH is thus a good alternative to the commonly recommended Boussinesq-type models for landslide-tsunami propagation.
Examples relying on a coupling approach include Narayanaswamy et al. [33], who proposed a hybrid SPHysics-FUNWAVE model, and Altomare et al. [34], who proposed a hybrid DualSPHysics-SWASH model to study wave propagation from offshore to nearshore. Viroulet et al. [35] used SPHysics for landslide-tsunami generation, and Gerris, which is a two-phase finite volume approach, for tsunami propagation. However, these simulations [33][34][35] were conducted in two dimensions. Three-dimensional (3D) simulations, as presented by Abadie et al. [12], are much rarer. They combined a Navier-Stokes model (THETIS) and FUNWAVE-TVD to study the 3D tsunami generation and propagation from the potential collapse of the Cumbre Vieja Volcano, La Palma, on the Canary Islands.
In this work, a 3D landslide-tsunami hazard assessment technique is illustrated by taking advantage of the two distinct models DualSPHysics (for slide impact and wave generation) and SWASH (for wave propagation). The feasibility of the proposed coupling method is demonstrated by investigating hypothetical landslide-tsunami scenarios originating from Es Vedrà, offshore Ibiza, although currently no evidence suggests that this island may become unstable. Es Vedrà has been selected because of its steep flanks exceeding typical basal friction angles >30° and for its proximity to other islands, allowing for run-up and inundation investigations. That landslide tsunamis at Es Vedrà have not been investigated previously added also to the motivation to select this particular case. More importantly, the developed techniques in this work are transferable to other cases in reservoirs, lakes, fjords, and the sea. Figure 1 shows the plan view of Ibiza and Formentera with Es Vedrà, including the paths of two critical landslide-tsunami scenarios investigated herein, and a picture of Es Vedrà taken from Ibiza. These two scenarios were selected because inhabitants, tourists, and infrastructure at these low-lying locations are particularly exposed to potential tsunami run-up in contrast to more elevated areas elsewhere on Ibiza and Formentera. These two scenarios involve tsunamis propagating towards (i) Cala d'Hort (3 km away from Es Vedrà) and (ii) Marina de Formentera (23 km away from Es Vedrà).  [36] and topography [37] in proximity of Es Vedrà with the range covered by the picture in (c) marked with blue dashed lines. (c) Picture of Es Vedrà taken by V. Heller from Ibiza.  [36] and topography [37] in proximity of Es Vedrà with the range covered by the picture in (c) marked with blue dashed lines. (c) Picture of Es Vedrà taken by V. Heller from Ibiza.
The remaining sections of this article are organised as follows. Section 2 presents the basic features of DualSPHysics and SWASH. The results of the convergence tests, as well as the wave generation and propagation in the two scenarios, are described in Section 3. In Section 4, the results are compared with another study, the implications of this work on landslide-tsunami hazard assessment are discussed, and the limitations are highlighted. The main conclusions are finally presented in Section 5.

Methods
The slide movement and landslide-tsunami generation were conducted with DualSPHysics v4.0, University of Vigo, weakly coupled with SWASH v4.01, Delft University of Technology, which was used for the wave propagation. The numerical backgrounds of these two open source codes are introduced in Sections 2.1 and 2.2.

Basic Principles
The open source code DualSPHysics v4.0 was used herein [30], which is based on weakly compressible SPH (WCSPH). The fluid phase in DualSPHysics is governed by the Navier-Stokes equations with the partial differential equations reduced to ordinary differential equations in a Lagrangian framework. The conservation of mass and momentum are expressed in Equations (1) and (2) [31,38]: where ρ is the density, v is the velocity vector, p is the pressure, g is the gravitational acceleration vector, and Γ is the dissipative term. DualSPHysics provides two options for dissipative terms, namely: artificial viscosity and laminar viscosity with sub-particle scale turbulence [30]. These two viscosity treatments provided similar results for the dam break case [39]. Thus, the artificial viscosity is widely used in landslide-tsunami modelling, resulting in the reasonable agreement of the numerical wave profiles with experimental results [8,40]. Based on this, the artificial viscosity was applied in the present work. In SPH, the fluid is discretised into a set of particles carrying properties such as density, pressure, velocity, position, etc. In general, two steps are required to transform Equations (1) and (2) into the SPH formalism, i.e., a kernel approximation and a particle approximation [41]. In the first step, any field variable f associated with particle a (located at x a ) can be represented by an integral at location x in the form of Equation (3): In Equation (3), Ω is the computation domain, W is the weighting function or smoothing kernel, which monotonically decreases with distance, and h p is the smoothing length determining the size of the kernel support. The kernels cubic spline and Wendland are available in DualSPHysics, and in this work, the former has been used.
In the second step, the integral in Equation (3) is approximated by interpolating the characteristics of the surrounding particles resulting in Equation (4): where the summation is over all the particles in the kernel support. W ab is short for W x a − x b , h p , and m b and ρ b are the mass and density, respectively, of particle b (located at x b ). Similarly, the derivative of the field variable f can be expressed as shown in Equation (5): where ∇ a is the derivative with respect to the coordinates of particle a.

Governing Equations
The governing equations in the SPH formalism are briefly discussed with details given by Crespo et al. [30]. The WCSPH method, although it is able to resolve the fluid kinematics well, suffers from unphysical pressure oscillations. In WCSPH, the pressure is linked to the density by means of an equation of state. Hence, remedies are proposed to stabilise the density field, and thus ensuring that the pressure field is noise-free. This is achieved with a density filter [38] and a density diffusion correction [42]. In the density diffusion correction, which is also referred to as the delta SPH algorithm, a diffusion term is added into the continuity equation to eliminate the numerical noise on the pressure field.
The continuity equation including the density diffusion correction is written in Equation (6): where subscript a and b denote fluid particles a and b, respectively, ρ a is the density of the fluid particle a, v ab = v a − v b , v a and v b are the velocity vectors of fluid particles a and b, respectively, δ is the delta SPH coefficient, c ab = 0.5(c a + c b ), c a and c b are the speed of sound at the locations of fluid particles a and b, respectively, x ab = x a − x b , x a and x b are the position vector of fluid particles a and b, respectively, and 0.01h 2 p is added to prevent singularities. The second term at the right-hand side of Equation (6) is the density diffusion term.
The momentum equation is written in the form of Equation (7) [30,38]: where p a and p b are the pressure of fluid particles a and b, respectively, and Π ab is the artificial viscosity, which accounts for the effects of dissipation, and is given by Equation (8): In Equation (8), κ is the artificial viscosity coefficient, µ ab = h p v ab ·x ab / x 2 ab + 0.01h 2 p and ρ ab = 0.5(ρ a + ρ b ).
In contrast to incompressible SPH where the pressure is solved by the Poisson equation, the pressure in WCSPH is determined via an equation of state. This Equation (9) relates the pressure to the density of the fluid to close the governing equation system [30,38]: In Equation (9), γ is typically selected as 7, ρ 0 = 1000 kg/m 3 is the reference density, and c 0 = .c(ρ 0 ) = ∂p a /∂ρ a | ρ 0 is the speed of sound at the reference density. The relative density fluctuation can be related to c 0 and the maximum velocity of the SPH particles v max , namely |∆ρ|/ρ ∼ v 2 max /c 2 0 [43]. Therefore, c 0 should be set to be at least 10 times larger than the maximum flow velocity such that the relative density fluctuation can be constrained to less than 1%, satisfying the condition that the fluid is nearly incompressible. Equations (6), (7), and (9) are complete, and result in the density, velocity, and pressure of each fluid particle by using an integration scheme. In DualSPHysics, the two integration schemes Verlet and Symplectic are provided [30], whilst this work adopted the Verlet scheme.
Each fluid particle position is then updated according to its velocity as shown in Equation (10): where ε = 0.5. The second term on the right-hand side of Equation (10) is the XSPH variant [44], which is aimed at moving each fluid particle with a velocity close to the average of its neighbourhood. The representation of solid boundaries in SPH is challenging primarily due to the truncation of the kernel support near a boundary. In general, solid boundaries are represented by particles, and three kinds of boundary particles are widely used in the technical literature, namely: (i) repulsive, (ii) ghost, and (iii) dynamic particles.
In terms of (i), repulsive boundary conditions, a repulsion force, normally in the Lennard-Jones form, is felt by those fluid particles in proximity of the boundary to prevent its motion beyond the domain of interest [45]. However, the repulsive particles have no contribution to the density of the fluid particles, which may further contaminate the pressure field. (ii) Ghost particles, which can either be dynamically generated or predefined, are placed beyond the boundary to fill the truncated domain of the kernel [46]. The field variables of the ghost particles are then mirrored from the inner fluid particles. Method (ii) has proven to be effective for simple boundaries, but encounters difficulties in dealing with more complex boundaries [31]. In (iii), dynamic boundary conditions [47], the dynamic boundary particles are forced to satisfy the same equations as the fluid particles, but are fixed in position, providing a sufficient repelling force to the fluid particles nearby. In this work, dynamic boundary conditions were used, given that this is currently the only option in DualSPHysics [30].
The motion of the landslide is fully resolved with no initial velocity imposed. The solid phase (landslide) is assumed to be rigid [13] and treated as a floating object [30,38]. The translational and rotational motion was obtained according to the force and momentum exerted on the floating object, i.e., by boundaries or ambient fluid particles. All of the relevant parameters that are used in DualSPHysics are summarised in Table 1. This includes the particle spacing, number of particles, ratio of smoothing length to particle spacing, delta-SPH coefficient, artificial viscosity coefficient, reference density ρ 0 , dimensionless parameter γ, ε, coefficient of speed of sound, number of time steps applied to Eulerian equations, the Courant-Friedrichs-Lewy (CFL) number, and physical time. Table 1. Parameters used in DualSPHysics. The values marked with * were only used in the convergence tests. CFL: Courant-Friedrichs-Lewy, SPH: smoothed particle hydrodynamics.

Scenario 1 Scenario 2
Particle spacing dp (mm) 10. The 3D numerical models of scenario 1 (Cala d'Hort) and scenario 2 (Marina de Formentera) in DualSPHysics were reconstructed in a similar manner. First, the bathymetric data available from the European Marine Observation and Data Network [36] and topographic data from the Spanish National Center of Geographic Information (CNIG) [37] were visualised and superimposed with each other. Second, the superimposed bathymetry and topography were interpolated such that the resolution was 5 m × 5 m. Third, the 4000 m × 4000 m large domain of interest was selected, and the slide plane was defined on Es Vedrà. Fourth, the slide profile and the interpolated bathymetry and topography were scaled at 1:500 before importing them into DualSPHysics. At this reduced scale, the simulations were run with the same particle spacing and similar SPH parameters as in Heller et al. [8], i.e., some experience was available. Last, the water body was added into the scaled model.
In DualSPHysics, a series of wave probes were placed at x = 200 m, 300 m, 500 m, 750 m, 1000 m, etc., to record wave profiles along the slide axis, with the coordinate origin for x placed at the intersection of the still water level (SWL) with the slide axis (white circles in Figures 1b and 2). The coupling location was selected by inspecting the numerical wave profiles, and the first location was selected where both the wave crest and trough are fully developed such that SWASH can cope with the input time series. This is the case at x = 750 m in scenario 1 and at x = 1000 m in scenario 2. Alternative coupling criteria, e.g., based on the location of the maximum wave amplitude defined in Heller and Hager [14], are discussed in Ruffini et al. [48].
One wave probe was placed on the slide axis, and 25 were placed perpendicularly on either side on a straight line. This resulted in a total of 51 wave probes separated by 30 m. The wave profiles and wave kinematics outputs from DualSPHysics were used as input in SWASH using a time series for each of the 51 grid points and depth-averaged velocities over both layers. The total length of the wave generation boundary in both scenarios was 1500 m.

Numerical Background
SWASH v4.01 [9,10,49] was used in this study to simulate the wave propagation in both scenarios. SWASH is a numerical model based on the depth-averaged non-hydrostatic NLSWEs given in Equations (11)-(13): ∂η ∂t ∂v ∂t + u ∂v ∂x + v ∂v ∂y + g ∂η ∂y where t is time, x, y, and z are the coordinates located at the mean SWL, g is the gravitational acceleration, h(x, y) is the still water depth, η(x, y, t) is the water surface elevation from the SWL, and d = h + η is the total water depth. u and v are the depth-averaged flow velocities in the two main directions. τ xx = 2ν t ∂u/∂x, τ xy = τ yx = ν t (∂v/∂x + ∂u/∂y), and τ yy = 2ν t ∂v/∂y are the horizontal turbulent stresses where ν t (x, y, t) is the horizontal eddy viscosity defined in Zijlema et al. [10]. c f is the bottom friction coefficient defined by Manning's formula, and q(x, y, z, t) is the non-hydrostatic pressure. Equations (11)-(13) were expanded for the multi-layer case by Stelling and Zijlema [48] using a discretisation method based on the Keller-box scheme.
The non-hydrostatic pressure is defined as a term of the total pressure in Equation (14) [50]: where p h is the hydrostatic term. The hydrostatic balance is given by Equation (15): The computation of the integral of the non-hydrostatic pressure gradient in Equations (12) and (13) is introduced in Zijlema et al. [10] where the free surface boundary condition of the non-hydrostatic pressure is q z=η = 0 and at the bottom the non-hydrostatic pressure is defined by applying the Keller-box method. Then, the vertical velocities at the free surface w s and at the bottom w b are considered with their respective momentum equations. Here, the vertical acceleration is instantaneously determined from the non-hydrostatic pressure. Finally, by combining the vertical momentum equations with the non-hydrostatic pressure equation at the bottom and using the kinematic bottom boundary condition w b = −u∂h/∂x − v∂h/∂y, the conservation of local mass yields Equation (16): Equation (16) closes the equation system and enables, together with the boundary conditions, to solve Equations (11) to (13).

Numerical Model Setup and Boundary Conditions
The use of SWASH v4.01 for the wave propagation allowed reducing the computational time significantly in both scenarios, and helped to avoid the wave decay problem of DualSPHysics (Section 4.2). Two different regular grids were created, one for each scenario, both with a grid resolution of ∆x = ∆y = 30 m, allowing at least 30 grid points per wavelength, which is sufficient to ensure the convergence of the solution [48,49]. The bathymetry was retrieved from the European Marine Observation and Data Network [36] with a resolution of approximatively 200 m. The topography was taken from the CNIG [37] with digital terrain model data with an accuracy of 25 m to match the grid resolution closely. The created domains counted 94 × 94 grid points for scenario 1 (Cala d'Hort) and 900 × 400 grid points for scenario 2 (Marina de Formentera). The Delft3D QUICKIN v5.0, Deltares, Delft, has been used to create the grids by processing the data with a triangular interpolation. Furthermore, smoothing of the interpolation has been applied on the west side of the Cap Llentrisca topography to avoid instabilities in the simulations from the large slopes at this location.
The code has been compiled for the use with multiple CPUs. All of the simulations were performed on the Nottingham High Performance Computing (HPC) cluster Minerva. For the biggest domain investigated, a real time of 25 min took 8.25 h of simulation time using four CPU cores and 10 GB of RAM. To solve the equations using multiple cores, the model divided the computational domain into subdomains. A stripwise decomposition method along the propagation direction of the tsunami was chosen for these purposes.
All of the simulations were performed using two layers, which generates a maximum error of 1% in the phase velocity for waves with kd ≤ 3 [10], where k specifies the wave number. An accurate simulation of the frequency dispersion can be very important for landslide tsunamis, as they can be highly dispersive [32]. In addition, the breaking criterion with default parameters has been applied to ensure that the energy dissipation due to wave breaking is accurately simulated.
Both water surface and a depth-averaged velocity time series were used as input for each point, as described in Section 2.1. The water surface time series was assigned using a weakly reflective boundary condition [51], and the depth-averaged velocity was assigned as a time series using velocity components that were directed perpendicularly to the coupling line. To achieve this, the velocity output from DualSPHysics was decomposed, and only the required component was used. However, it is expected that this did not affect the results significantly, as the grids were oriented along the main tsunami propagation direction. Furthermore, all of the open sides of the domain were defined using a Sommerfeld boundary condition [52], which allows waves to cross the outflow boundary without reflections, and is particularly suited for long waves. This condition, for a boundary oriented parallel to the y-axis, is given in Equation (17): The Manning's roughness coefficient n was used to include bottom friction via Equation (18): This study applies the default value n = 0.019 s/m 1/3 as recommended for wave simulation on sandy beaches [53]. Further, a value of 5 mm is chosen as the threshold for the minimum computed water depth.
Finally, an explicit time discretisation method is used. This method uses the CFL condition to adjust the time step during the simulation accordingly. The Courant number C r for the performed simulations is defined in Equation (19): where ∆x and ∆y are the distances between two grid points in the x and y directions, respectively.
To calculate whether to increase or reduce the time step, a minimum and maximum C r threshold can be applied in the simulation in order to control the convergence of the solution accurately. C r,min = 0.1 and C r,max = 0.5 are used herein as suggested for large, nonlinear waves and wave interaction with structures and steep slopes [53]. Figure 2 shows the two-dimensional (2D) slide profiles along the main wave propagation directions (Figure 1a). The origins are defined at the intersections of the island, the slide axis and the SWL (Figure 1b). Es Vedrà is mainly composed of Miocene evaporates [54], which is a special type of limestone. The Miocene limestone has a density of 1210 kg/m 3 to 2510 kg/m 3 [55]. In this work, the slide density is assumed to be 1600 kg/m 3 . The slide profile is obtained by cutting the island at an angle of 30 • , assuming that the sliding surface is planar. Smaller and larger slope angles in 5 • intervals were also investigated. Smaller slope angles result in larger slide volumes, but in smaller slide impact velocities. No slide movement was observed for very flat angles. In scenario 1, 30 • resulted in the largest tsunamis. In scenario 2, both slope angles 25 • and 30 • resulted essentially in the same wave heights such that again, 30 • was selected to be consistent with scenario 1. The slide volume in scenario 1 was 4.33 × 10 6 m 3 , corresponding to a mass of 6.93 × 10 9 kg and in scenario 2, the volume was 6.92 × 10 6 m 3 and the mass 11.08 × 10 9 kg.

Slide Scenarios
This study applies the default value = 0.019 s/m / as recommended for wave simulation on sandy beaches [53]. Further, a value of 5 mm is chosen as the threshold for the minimum computed water depth.
Finally, an explicit time discretisation method is used. This method uses the CFL condition to adjust the time step during the simulation accordingly. The Courant number for the performed simulations is defined in Equation (19): where Δ and Δ are the distances between two grid points in the and directions, respectively. To calculate whether to increase or reduce the time step, a minimum and maximum threshold can be applied in the simulation in order to control the convergence of the solution accurately. , = 0.1 and , = 0.5 are used herein as suggested for large, nonlinear waves and wave interaction with structures and steep slopes [53].

Wave Generation
3.1.1. Slide Scenarios Figure 2 shows the two-dimensional (2D) slide profiles along the main wave propagation directions ( Figure 1a). The origins are defined at the intersections of the island, the slide axis and the SWL (Figure 1b). Es Vedrà is mainly composed of Miocene evaporates [54], which is a special type of limestone. The Miocene limestone has a density of 1210 kg/m 3 to 2510 kg/m 3 [55]. In this work, the slide density is assumed to be 1600 kg/m 3 . The slide profile is obtained by cutting the island at an angle of 30°, assuming that the sliding surface is planar. Smaller and larger slope angles in 5° intervals were also investigated. Smaller slope angles result in larger slide volumes, but in smaller slide impact velocities. No slide movement was observed for very flat angles. In scenario 1, 30° resulted in the largest tsunamis. In scenario 2, both slope angles 25° and 30° resulted essentially in the same wave heights such that again, 30° was selected to be consistent with scenario 1. The slide volume in scenario 1 was 4.33 × 10 6 m 3 , corresponding to a mass of 6.93 × 10 9 kg and in scenario 2, the volume was 6.92 × 10 6 m 3 and the mass 11.08 × 10 9 kg.

Convergence Tests
Convergence tests have been conducted for scenario 2. Five wave probes were placed along the slide axis at 0.40 m, 0.60 m, 1.00 m, 1.50 m, and 2.00 m. Four particle resolutions were selected with the resulting wave profiles at 2.00 m, as shown in Figure 3, thereby focusing on the primary wave. The results in Figure 3 are shown at scale 1:500, while all of the other results in this article are upscaled to the nature scale with Froude scaling [56,57]. Generally speaking, the wave profiles converge with finer particle resolution towards a larger wave amplitude, except for dp = 7.5 mm. The wave heights for dp = 7.5 mm, 10 mm, 15 mm, and 20 mm are 0.137 m, 0.142 m, 0.119 m, and 0.101 m, respectively. In relation to the wave height observed at dp = 7.5 mm, this results in differences of 3.1%, −13.3%, and −26.5% for dp = 10 mm, 15 mm, and 20 mm, respectively. The discrepancy of the wave height between dp = 7.5 mm and 10 mm is relatively small, indicating convergence. Based on this, the particle resolution 10 mm was selected as a compromise between accuracy and computational cost. This resolution was then applied to both scenarios, which resulted in approximately 8.40 million particles in scenario 1 and 11.74 million particles in scenario 2 ( Table 1) to the nature scale with Froude scaling [56,57]. Generally speaking, the wave profiles converge with finer particle resolution towards a larger wave amplitude, except for dp = 7.5 mm. The wave heights for dp = 7.5 mm, 10 mm, 15 mm, and 20 mm are 0.137 m, 0.142 m, 0.119 m, and 0.101 m, respectively. In relation to the wave height observed at dp = 7.5 mm, this results in differences of 3.1%, −13.3%, and −26.5% for dp = 10 mm, 15 mm, and 20 mm, respectively. The discrepancy of the wave height between dp = 7.5 mm and 10 mm is relatively small, indicating convergence. Based on this, the particle resolution 10 mm was selected as a compromise between accuracy and computational cost. This resolution was then applied to both scenarios, which resulted in approximately 8.40 million particles in scenario 1 and 11.74 million particles in scenario 2 ( Table 1)  Convergence test: wave profiles at 2.00 m along the slide axis for different particle spacings dp in scenario 2. Figure 4 shows the slide kinematics in both scenarios upscaled to the nature scale. The slide front impact velocity in scenario 1 is 49.48 m/s, which is close to the maximum slide velocity. In scenario 2, the slide front impact velocity is 13.02 m/s, with the slide velocity further increasing by approximately a factor of four. It is emphasised that the slide velocity in DualSPHysics may overestimate the real slide velocity [8]. This is because the slide block was treated as a floating object, and the basal friction was not fully modelled (Section 4.2). However, the maximum slide velocities shown in Figure 4 are not unrealistic for landslides, e.g., the slide impact velocity of the 1958 Lituya Bay rockslide was estimated at 92 m/s [14], and of the 2007 Chehalis landslide was estimated at 60 m/s [58]. The velocities in Figure 4 are considered to be extreme case scenarios where the initial friction coefficients are reduced, which is a phenomenon known as hypermobility, resulting sometimes in anomalously high velocities and long runout distances in nature for slide volumes >10 6 m 3 [59]. The slide position time histories, which were directly derived from the slide velocity time histories, are also shown in Figure  4. The slide positions from initiation to deposition in both scenarios exceed 1000 m, which further indicates hypermobility features of the herein investigated landslides.  Convergence test: wave profiles at 2.00 m along the slide axis for different particle spacings dp in scenario 2. Figure 4 shows the slide kinematics in both scenarios upscaled to the nature scale. The slide front impact velocity in scenario 1 is 49.48 m/s, which is close to the maximum slide velocity. In scenario 2, the slide front impact velocity is 13.02 m/s, with the slide velocity further increasing by approximately a factor of four. It is emphasised that the slide velocity in DualSPHysics may overestimate the real slide velocity [8]. This is because the slide block was treated as a floating object, and the basal friction was not fully modelled (Section 4.2). However, the maximum slide velocities shown in Figure 4 are not unrealistic for landslides, e.g., the slide impact velocity of the 1958 Lituya Bay rockslide was estimated at 92 m/s [14], and of the 2007 Chehalis landslide was estimated at 60 m/s [58]. The velocities in Figure 4 are considered to be extreme case scenarios where the initial friction coefficients are reduced, which is a phenomenon known as hypermobility, resulting sometimes in anomalously high velocities and long runout distances in nature for slide volumes >10 6 m 3 [59]. The slide position time histories, which were directly derived from the slide velocity time histories, are also shown in Figure 4. The slide positions from initiation to deposition in both scenarios exceed 1000 m, which further indicates hypermobility features of the herein investigated landslides.

Analysis of the Results
The wave propagations at the nature scale are depicted in Figure 5 (scenario 1) and Figure 6 (scenario 2). In both scenarios, the water is fully displaced at the wake of the landslide. This is an artifact of the discrete SPH particles, and is not expected to occur to such a degree in reality. In scenario 1, before t = 16.9 s, the slide moves below the primary wave crest, indicating that the slide motion is faster than the wave. This is also observed in scenario 2, where the primary wave finally overtakes the landslide at t = 20.1 s, and continues to propagate outwards (Figure 6d). Both Figures 5 and 6 confirm that the landslide tsunamis are largest on the slide axis and significantly smaller at the peripheries [29,48]. estimated at 92 m/s [14], and of the 2007 Chehalis landslide was estimated at 60 m/s [58]. The velocities in Figure 4 are considered to be extreme case scenarios where the initial friction coefficients are reduced, which is a phenomenon known as hypermobility, resulting sometimes in anomalously high velocities and long runout distances in nature for slide volumes >10 6 m 3 [59]. The slide position time histories, which were directly derived from the slide velocity time histories, are also shown in Figure  4. The slide positions from initiation to deposition in both scenarios exceed 1000 m, which further indicates hypermobility features of the herein investigated landslides.  The wave propagations at the nature scale are depicted in Figure 5 (scenario 1) and Figure 6 (scenario 2). In both scenarios, the water is fully displaced at the wake of the landslide. This is an artifact of the discrete SPH particles, and is not expected to occur to such a degree in reality. In scenario 1, before t = 16.9 s, the slide moves below the primary wave crest, indicating that the slide motion is faster than the wave. This is also observed in scenario 2, where the primary wave finally overtakes the landslide at t = 20.1 s, and continues to propagate outwards (Figure 6d). Both Figures 5  and 6 confirm that the landslide tsunamis are largest on the slide axis and significantly smaller at the peripheries [29,48].  Four wave probes in scenario 1 and five wave probes in scenario 2 were placed along the slide axes ( Figure 1) to record the landslide tsunamis, as shown in Figure 7. The flat parts in Figure 7 indicate that the water at these wave probe locations is fully displaced, as discussed in the last paragraph. The maximum wave height in scenario 1 of 133.0 m is larger than the 75.4 m that is observed in scenario 2 ( Figure 7). However, the landslide mass is with 6.93 × 10 9 kg smaller in scenario 1 than in scenario 2, in which it is 11.08 × 10 9 kg. This seemingly contradiction can be explained with the larger slide Froude number v/(gh) 1/2 in scenario 1, which has a significantly more dominant effect on the tsunami magnitude than the slide mass for subaerial landslides, according to the impulse product parameter [14].
The maximum wave amplitude in scenario 1 was more than twice as large as the water depth (Figure 2a). This is not unusual in the slide impact zone for large slide impact velocities (see e.g., Figure 7a,b in Huang et al. [60] or Figure 5a in Heller and Hager [14], where the relative maximum wave amplitudes were similarly large). Such large "waves" are the consequences of the impact crater, water splash, and the slide trapped below the primary wave, and does not represent a stable wave. Figure 7 also reveals that the wave decays faster in scenario 1; the water column is elevated by the landslide located below the primary wave at the impact zone, and the sudden drop in wave elevation may be related to the disappearance of this effect once the wave travels faster than the slide. The wave propagations at the nature scale are depicted in Figure 5 (scenario 1) and Figure 6 (scenario 2). In both scenarios, the water is fully displaced at the wake of the landslide. This is an artifact of the discrete SPH particles, and is not expected to occur to such a degree in reality. In scenario 1, before t = 16.9 s, the slide moves below the primary wave crest, indicating that the slide motion is faster than the wave. This is also observed in scenario 2, where the primary wave finally overtakes the landslide at t = 20.1 s, and continues to propagate outwards (Figure 6d). Both Figures 5  and 6 confirm that the landslide tsunamis are largest on the slide axis and significantly smaller at the peripheries [29,48].  Four wave probes in scenario 1 and five wave probes in scenario 2 were placed along the slide axes ( Figure 1) to record the landslide tsunamis, as shown in Figure 7. The flat parts in Figure 7 indicate that the water at these wave probe locations is fully displaced, as discussed in the last paragraph. The maximum wave height in scenario 1 of 133.0 m is larger than the 75.4 m that is observed in scenario 2 ( Figure 7). However, the landslide mass is with 6.93 × 10 9 kg smaller in scenario 1 than in scenario 2, in which it is 11.08 × 10 9 kg. This seemingly contradiction can be explained with the larger slide Froude number v/(gh) 1/2 in scenario 1, which has a significantly more dominant effect on the tsunami magnitude than the slide mass for subaerial landslides, according to the impulse product parameter [14].
The maximum wave amplitude in scenario 1 was more than twice as large as the water depth (Figure 2a). This is not unusual in the slide impact zone for large slide impact velocities (see e.g., Figure 7a,b in Huang et al. [60] or Figure 5a in Heller and Hager [14], where the relative maximum wave amplitudes were similarly large). Such large "waves" are the consequences of the impact crater, water splash, and the slide trapped below the primary wave, and does not represent a stable wave. Figure 7 also reveals that the wave decays faster in scenario 1; the water column is elevated by the landslide located below the primary wave at the impact zone, and the sudden drop in wave elevation may be related to the disappearance of this effect once the wave travels faster than the slide. Four wave probes in scenario 1 and five wave probes in scenario 2 were placed along the slide axes ( Figure 1) to record the landslide tsunamis, as shown in Figure 7. The flat parts in Figure 7 indicate that the water at these wave probe locations is fully displaced, as discussed in the last paragraph. The maximum wave height in scenario 1 of 133.0 m is larger than the 75.4 m that is observed in scenario 2 ( Figure 7). However, the landslide mass is with 6.93 × 10 9 kg smaller in scenario 1 than in scenario 2, in which it is 11.08 × 10 9 kg. This seemingly contradiction can be explained with the larger slide Froude number v/(gh) 1/2 in scenario 1, which has a significantly more dominant effect on the tsunami magnitude than the slide mass for subaerial landslides, according to the impulse product parameter [14].
The maximum wave amplitude in scenario 1 was more than twice as large as the water depth (Figure 2a). This is not unusual in the slide impact zone for large slide impact velocities (see e.g., Figure 7a,b in Huang et al. [60] or Figure 5a in Heller and Hager [14], where the relative maximum wave amplitudes were similarly large). Such large "waves" are the consequences of the impact crater, water splash, and the slide trapped below the primary wave, and does not represent a stable wave. Figure 7 also reveals that the wave decays faster in scenario 1; the water column is elevated by the landslide located below the primary wave at the impact zone, and the sudden drop in wave elevation may be related to the disappearance of this effect once the wave travels faster than the slide.   (Figure 8a), which was mainly due to a second wave caused by reflection from the cliffs surrounding the beach (Figure 9) and shoaling effects. The effect of the reflection can be seen at xb = − 100 m and t = 135 s in Figure 8b.    Figure 8a shows the time series of the water surfaces measured every 200 m along the slide axis for 10 locations in total. The primary wave remains the largest wave up to the beach with a maximum wave amplitude of 28.5 m at x = 850 m. This is around 40% less than the input value that was used in SWASH at x = 750 m, which is consistent with the decay found between x = 500 m and x = 750 m in DualSPHysics (Figure 7a). The minimum primary wave amplitude in Figure 8a (Figure 8a), which was mainly due to a second wave caused by reflection from the cliffs surrounding the beach (Figure 9) and shoaling effects. The effect of the reflection can be seen at xb = − 100 m and t = 135 s in Figure 8b.   Figure 9 shows the wave propagation towards Cala d'Hort with four snapshots at different times. This shows that the tsunami front remains circular and the maximum amplitude occurs on the slide axis. Further, Figure 9d shows that the wave reaches the shore 122 s after slide impact. Figure 9c,d clearly confirm the reflections on the southeast cliffs generating the secondary peak in Figure 8b.

Scenario 1 (Cala d'Hort)
The tsunami run-up is investigated next with Figure 10, and five critical points are examined in more detail. The values shown refer to the wet computational grid points with P1 and P4 showing the wet most inland points. Given that the grid resolution is only 30 m, the maximum run-up height is likely to be underestimated by the specified values. P1 is positioned inland of Cala Carbo, which is a smaller beach north of Cala d'Hort where at a terrain elevation of 6 m a maximum inundation depth of 0.23 m is measured. P2 is positioned on the west side of the island Escull de Cala d'Hort where the tsunami generates a run-up height of 5.7 m with a maximum water depth of 4.7 m. The last three points in Figure 10 are all positioned at Cala d'Hort. Point B is on the slide axis near the beach shore where a maximum inundation depth of 14.43 m on a terrain elevation of 0.42 m is measured. The water level starts to decrease towards P3 at the most inland part of the beach with a run-up of 9.24 m and a maximum inundation depth of 11.5 m over the entire event. Finally, the water reaches the maximum terrain elevation of 19.7 m above SWL at P4 and a maximum inundation depth of 1.8 m.
J. Mar. Sci. Eng. 2018, 9, x FOR PEER REVIEW 13 of 22 Figure 9 shows the wave propagation towards Cala d'Hort with four snapshots at different times. This shows that the tsunami front remains circular and the maximum amplitude occurs on the slide axis. Further, Figure 9d shows that the wave reaches the shore 122 s after slide impact. Figure  9c,d clearly confirm the reflections on the southeast cliffs generating the secondary peak in Figure 8b.
The tsunami run-up is investigated next with Figure 10, and five critical points are examined in more detail. The values shown refer to the wet computational grid points with P1 and P4 showing the wet most inland points. Given that the grid resolution is only 30 m, the maximum run-up height is likely to be underestimated by the specified values. P1 is positioned inland of Cala Carbo, which is a smaller beach north of Cala d'Hort where at a terrain elevation of 6 m a maximum inundation depth of 0.23 m is measured. P2 is positioned on the west side of the island Escull de Cala d'Hort where the tsunami generates a run-up height of 5.7 m with a maximum water depth of 4.7 m. The last three points in Figure 10 are all positioned at Cala d'Hort. Point B is on the slide axis near the beach shore where a maximum inundation depth of 14.43 m on a terrain elevation of 0.42 m is measured. The water level starts to decrease towards P3 at the most inland part of the beach with a run-up of 9.24 m and a maximum inundation depth of 11.5 m over the entire event. Finally, the water reaches the maximum terrain elevation of 19.7 m above SWL at P4 and a maximum inundation depth of 1.8 m.

Scenario 2 (Marina de Formentera)
The tsunami propagation and run-up analysis was performed similarly for scenario 2 as for scenario 1. However, the distance between the slide impact and Marina de Formentera is around 10 times larger. Therefore, Figure 11a shows the water surface elevations on the slide axis for every 2000 m only. The maximum amplitude at x = 1500 m is 31 m. The wave decays fast until x = 13,500 m; afterwards, the amplitude stabilises at approximately 2.5 to 3 m whilst approaching the harbour, as the wave decay may be compensated by shoaling effects.

Scenario 2 (Marina de Formentera)
The tsunami propagation and run-up analysis was performed similarly for scenario 2 as for scenario 1. However, the distance between the slide impact and Marina de Formentera is around 10 times larger. Therefore, Figure 11a shows the water surface elevations on the slide axis for every 2000 m only. The maximum amplitude at x = 1500 m is 31 m. The wave decays fast until x = 13,500 m; afterwards, the amplitude stabilises at approximately 2.5 to 3 m whilst approaching the harbour, as the wave decay may be compensated by shoaling effects.

Scenario 2 (Marina de Formentera)
The tsunami propagation and run-up analysis was performed similarly for scenario 2 as for scenario 1. However, the distance between the slide impact and Marina de Formentera is around 10 times larger. Therefore, Figure 11a shows the water surface elevations on the slide axis for every 2000 m only. The maximum amplitude at x = 1500 m is 31 m. The wave decays fast until x = 13,500 m; afterwards, the amplitude stabilises at approximately 2.5 to 3 m whilst approaching the harbour, as the wave decay may be compensated by shoaling effects.    Figure 12 shows four snapshots of the tsunami propagation towards Formentera, where the primary tsunami impacts the harbour and coast of the island in Figure 12d after 15 min. It should be noted that the range of the legends in Figure 12b-d have been changed in relation to Figure 12a to reflect the strong reduction in wave amplitude. Figure 12b shows the tsunami diffraction around Cap Llentrisca on Ibiza, which subsequently results in reflection from the cliffs. J. Mar. Sci. Eng. 2018, 9, x FOR PEER REVIEW 15 of 22 Figure 11b shows how the tsunami approaches the tip of the harbour with a constant wave amplitude of 2.7 m due to the moderate sea bottom slope. Further, Figure 12 shows four snapshots of the tsunami propagation towards Formentera, where the primary tsunami impacts the harbour and coast of the island in Figure 12d after 15 min. It should be noted that the range of the legends in Figure 12b-d have been changed in relation to Figure 12a to reflect the strong reduction in wave amplitude. Figure 12b shows the tsunami diffraction around Cap Llentrisca on Ibiza, which subsequently results in reflection from the cliffs.

Comparison with La Palma Case
Abadie et al. [12] conducted a closely related numerical study by investigating a potential landslide tsunami from La Palma, Canary Islands, and expressed the tsunami decay as a function of the radial distance r from the slide impact location. The tsunami amplitudes a of the present study in both scenarios are compared in Figure 14 with the wave amplitude decays found in Abadie et al. [12]. The amplitudes are shown for the same positions as in Figures 7, 8a, and 11a. The solid and dashed lines represent the maximum and minimum decay found by Abadie et al. [12] for slide volumes of 80 km 3 and 450 km 3 , respectively. The tsunami amplitude of scenario 1 decays similarly to . , while the tsunami in scenario 2 decays with . , and lays close to the lower boundary found by Abadie et al. [12]. This small difference is likely due to the different bathymetries in the two cases. Scenario 2 shows a larger water depth and a very mildly sloped seabed, whereas in scenario 1, the shallower water depth and rapidly varying seabed interact more with the tsunami, increasing the decay rate. Mainland Ibiza, which is shown on the north side in Figure 12a, may also have a small effect on the lateral energy spread [29,48], and the wave decay may also change with the wave type (landslide-tsunami propagation in 3D typically involves Stokes-like and cnoidal-like waves [29]). Overall, a consistent wave amplitude decay between the present study and Abadie et al. [12] is found.

Comparison with La Palma Case
Abadie et al. [12] conducted a closely related numerical study by investigating a potential landslide tsunami from La Palma, Canary Islands, and expressed the tsunami decay as a function of the radial distance r from the slide impact location. The tsunami amplitudes a of the present study in both scenarios are compared in Figure 14 with the wave amplitude decays found in Abadie et al. [12]. The amplitudes are shown for the same positions as in Figures 7, 8a and 11a. The solid and dashed lines represent the maximum and minimum decay found by Abadie et al. [12] for slide volumes of 80 km 3 and 450 km 3 , respectively. The tsunami amplitude of scenario 1 decays similarly to r −1.19 , while the tsunami in scenario 2 decays with r −0.95 , and lays close to the lower boundary found by Abadie et al. [12]. This small difference is likely due to the different bathymetries in the two cases. Scenario 2 shows a larger water depth and a very mildly sloped seabed, whereas in scenario 1, the shallower water depth and rapidly varying seabed interact more with the tsunami, increasing the decay rate. Mainland Ibiza, which is shown on the north side in Figure 12a, may also have a small effect on the lateral energy spread [29,48], and the wave decay may also change with the wave type (landslide-tsunami propagation in 3D typically involves Stokes-like and cnoidal-like waves [29]). Overall, a consistent wave amplitude decay between the present study and Abadie et al. [12] is found. Figure 14. Comparison of wave amplitude a decay over radial distance r for both scenarios with maximum and minimum decay found by Abadie et al. [12].

Implications and Limitations
The main purpose of this work is to demonstrate a versatile numerical technique that can be applied to other landslide-tsunami events in reservoirs, lakes, fjords, and the sea. This numerical landslide-tsunami hazard assessment technique for site-specific bathymetric and topographic conditions relies on the two complementing numerical codes DualSPHysics and SWASH; DualSPHysics is well-suited for the violent wave generation process, but is computationally expensive, and may result in unphysical wave decay beyond the coupling location. On the other hand, SWASH deals well with long-distance wave propagation and run-up by taking frequency dispersion and bottom friction into account at small computational cost. It was demonstrated that these two codes combined are capable of resolving the entire landslide-tsunami phenomenon from slide impact, wave generation, and wave propagation to run-up, as well as the inundation of coastal areas.
The landslide-tsunami scenarios investigated herein are hypothetical, and no evidence currently suggests that Es Vedrà may become unstable. Further, the predicted maximum wave amplitudes herein would be extreme values. This is because only subaerial landslides were investigated, which are known to generate larger tsunamis than partially submerged or submarine landslides by given slide volume. Further, extreme slide scenarios resulting in the largest waves were investigated, and it was in addition not possible to fully model the slide basal friction at this stage (the slides are rather modelled as hypermobile [59]). This issue has already been pointed out by Heller et al. [8], who reduced the slide impact velocity by approximately 50% to match the experimental wave amplitudes. In the present work, the calibration and validation tests of DualSPHysics v4.0 with the experimental data presented in Heller et al. [8] was once more performed, resulting in similar conclusions as described by Heller et al. [8], namely that the numerical wave amplitudes are overestimated with an unreduced slide impact velocity. More work is required to fully address this slide kinematics issue.
In the meantime, it should be kept in mind that the simulated tsunamis are likely to be larger in the immediate slide impact zone than observed in reality. On the other hand, the wave decay in DualSPHysics was shown to be overpredicted in 3D by Heller et al. [8] (see their Figure 6e,f), such that the too-large slide velocity and wave decay partially compensate each other, and the simulated tsunami amplitudes at the coupling location are expected to be reliable enough for engineering applications. On the other hand, the wave propagation in SWASH is known to work reliably [48].
Finally, even though the wave propagation towards the mainland of Spain was not investigated herein, it is safe to state that a potential landslide tsunami originating from Es Vedrà would be very small at the mainland of Spain due to two main reasons. Firstly, Es Vedrà is less sloped towards the mainland of Spain, such that the slide volume and hence the tsunamis, would be smaller than in the two investigated scenarios. Secondly, the closest point on the mainland of Spain is 85 km away from

Implications and Limitations
The main purpose of this work is to demonstrate a versatile numerical technique that can be applied to other landslide-tsunami events in reservoirs, lakes, fjords, and the sea. This numerical landslide-tsunami hazard assessment technique for site-specific bathymetric and topographic conditions relies on the two complementing numerical codes DualSPHysics and SWASH; DualSPHysics is well-suited for the violent wave generation process, but is computationally expensive, and may result in unphysical wave decay beyond the coupling location. On the other hand, SWASH deals well with long-distance wave propagation and run-up by taking frequency dispersion and bottom friction into account at small computational cost. It was demonstrated that these two codes combined are capable of resolving the entire landslide-tsunami phenomenon from slide impact, wave generation, and wave propagation to run-up, as well as the inundation of coastal areas.
The landslide-tsunami scenarios investigated herein are hypothetical, and no evidence currently suggests that Es Vedrà may become unstable. Further, the predicted maximum wave amplitudes herein would be extreme values. This is because only subaerial landslides were investigated, which are known to generate larger tsunamis than partially submerged or submarine landslides by given slide volume. Further, extreme slide scenarios resulting in the largest waves were investigated, and it was in addition not possible to fully model the slide basal friction at this stage (the slides are rather modelled as hypermobile [59]). This issue has already been pointed out by Heller et al. [8], who reduced the slide impact velocity by approximately 50% to match the experimental wave amplitudes. In the present work, the calibration and validation tests of DualSPHysics v4.0 with the experimental data presented in Heller et al. [8] was once more performed, resulting in similar conclusions as described by Heller et al. [8], namely that the numerical wave amplitudes are overestimated with an unreduced slide impact velocity. More work is required to fully address this slide kinematics issue. In the meantime, it should be kept in mind that the simulated tsunamis are likely to be larger in the immediate slide impact zone than observed in reality. On the other hand, the wave decay in DualSPHysics was shown to be overpredicted in 3D by Heller et al. [8] (see their Figure 6e,f), such that the too-large slide velocity and wave decay partially compensate each other, and the simulated tsunami amplitudes at the coupling location are expected to be reliable enough for engineering applications. On the other hand, the wave propagation in SWASH is known to work reliably [48].
Finally, even though the wave propagation towards the mainland of Spain was not investigated herein, it is safe to state that a potential landslide tsunami originating from Es Vedrà would be very small at the mainland of Spain due to two main reasons. Firstly, Es Vedrà is less sloped towards the mainland of Spain, such that the slide volume and hence the tsunamis, would be smaller than in the two investigated scenarios. Secondly, the closest point on the mainland of Spain is 85 km away from Es Vedrà, and by applying the decay r −1.19 that was found towards Cala d'Hort (which is likely to underpredict the decay of the more freely propagating tsunamis towards the mainland of Spain), the resulting wave would only be 0.16 m based on the largest wave amplitude of 133 m in scenario 1. Further, the affected coast of mainland Spain consists mainly of cliffs, which also reduces the damage potential of a hypothetical landslide tsunami from Es Vedrà.

Conclusions
This work addressed a numerical technique to conduct landslide-tsunami hazard assessments in reservoirs, lakes, fjords, and the sea. The viability of this technique was demonstrated with hypothetical landslide tsunamis originating from Es Vedrà, offshore Ibiza, under consideration of the site-specific bathymetry and topography. The two numerical codes DualSPHysics and SWASH were applied, thereby combining the strengths in modelling the violent wave generation process of the former with the accurate long-distance wave propagation at the small computational cost of the latter. The coupling of the two codes was carried out by importing the wave profiles and wave kinematics from DualSPHysics into SWASH at a distance where the wave profiles were reasonably stable. Two landslide-tsunami scenarios were investigated by focusing on two particularly exposed locations, namely (i) Cala d'Hort (3 km away from Es Vedrà) and (ii) Marina de Formentera (23 km away from Es Vedrà). The main findings of this study are summarised as follows: -Two different slide-wave interaction phases were observed. (i) At the very beginning, the slide moved faster than the waves, such that the slide propagated below the primary wave crest and additionally elevated the water column and free water surface. (ii) The slide then slowed down such that the waves travelled faster and abruptly decayed due to the increased water depth. -In scenario 1 (Cala d'Hort), the maximum wave amplitude was 133.0 m, reducing to a wave amplitude of 14.2 m at 5 m offshore and a maximum run-up height of over 21.5 m. In scenario 2 (Marina de Formentera), the maximum wave amplitude was 75.4 m, reducing to 2.7 m at 5 m offshore Marina de Formentera, such that the inundation depth was 1.2 m in the populated harbor area. This is significantly smaller than at Cala d'Hort, but may still result in significant devastation due to a larger density of buildings and infrastructure. -The proposed numerical technique results likely in an overestimation of the landslide tsunamis because extreme slide scenarios were selected (extreme slide masses, slip orientation, and subaerial slides), and the slide velocity in DualSPHysics is likely to be overpredicted. -The proposed numerical technique also provided new insights into 3D landslide-tsunami propagation by considering site-specific bathymetric and topographic conditions.
The change of the landslide-tsunami features by using the laminar viscosity with sub-particle scale turbulence rather than the artificial viscosity for the dissipative terms, the slide kinematics issue and the definition of an exact criterion for the coupling location remain open for future research.