Denys Dutykh
email: denys.dutykh@ucd.ie

Henrik Kalisch
email: henrik.kalisch@math.uib.no

Boussinesq modeling of surface waves due to underwater landslides

Keywords: Surface waves, Boussinesq model, submarine landslides, wave run-up, tsunami PACS: 45.20.D, 47.11.Df, 47.35.Bb, 47.35.Fg, 47.85.Dh

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Surface waves originating from sudden perturbations of the bottom topography are often termed tsunamis. Two distinct generation mechanisms of a tsunami are underwater earthquakes, and submarine mass failures. Among the broad class of submarine mass failures, landslides can be characterised as translational failures which travel considerable distances along the bottom profile [36,[START_REF] Prior | Submarine landslides: geometry and nomenclature[END_REF]. In the past, the role of landslides and rock falls in the excitation of tsunamis may have been underestimated, as most known occurrences of tsunamis were accredited to seismic activity. However, it is now more accepted that submarine mass failures also contribute to a large portion of tsunamis [START_REF] Tinti | Tsunami Excitation by Submarine Slides in Shallow-water Approximation[END_REF], and recent years have seen a multitude of works devoted to the study of such underwater landslides and the resulting effect on surface waves [START_REF] Bardet | Landslide Tsunamis: Recent Findings and Research Directions[END_REF][START_REF] Chubarov | Numerical modelling of surface water waves arising due to movement of underwater landslide on irregular bottom slope[END_REF][START_REF] Didenkulova | Tsunami waves generated by submarine landslides of variable volume: analytical solutions for a basin of variable depth[END_REF][START_REF] Fernández-Nieto | A new Savage-Hutter type model for submarine avalanches and generated tsunami[END_REF]35,36,[START_REF] Okal | Normal mode energetics for far-field tsunamis generated by dislocations and landslides[END_REF]47,[START_REF] Poncet | A study of the tsunami effects of two landslides in the St. Lawrence estuary[END_REF][START_REF] Tinti | Tsunami Excitation by Submarine Slides in Shallow-water Approximation[END_REF]. As suggested in [START_REF] Fritz | Extreme runup from the 17 July 2006 Java tsunami[END_REF], it is possible for underwater landslides and earthquakes to act in tandem, and produce very large surface waves A natural question to ask is whether the effect of underwater landslides on surface waves can be such that they may pose a danger for civil engineering structures located near the shore. Consequently, one important issue is the wave action and in particular the run-up and draw-down at beaches in the vicinity of the landslide. While the draw-down itself may not pose a threat, one consequence of a large draw-down can be the amplification of the run-up of the following positive wave crest [START_REF] Dutykh | Dispersive wave runup on non-uniform shores[END_REF][START_REF] Tadepalli | Model for the leading waves of tsunamis[END_REF].

There have been many numerical and a few experimental studies devoted to this subject, but it is generally difficult to include many of the complex parameters and dependencies of a realistic landslide into a physical model. Therefore, most workers attempt to distill the problem to a model setup where many effects such as turbulence and sedimentation are disregarded. For example, Grilli and Watts [36] study tsunami sensitivity to several landslide parameters in the case of a landslide in a coastal area of an open ocean. In particular, dependence on the landslide shape and the initial depth of the landslide location are studied, and it is found that the landslide with the smallest length produced the largest waveheight and run-up, and that the wave run-up at an adjacent beach is inversely proportional to the initial depth. The work in [36] relies on integrating the full water-wave equations using an irrotational boundary-element code, and using an open boundary with transmission conditions [START_REF] Grilli | A fully nonlinear model for three-dimensional overturning waves over an arbitrary bottom[END_REF][START_REF] Grilli | Progress in fully nonlinear potential flow modeling of 3D extreme ocean waves[END_REF]. While most works have considered a given dynamics for the landslide, the bottom motion in [36] is described by an ordinary differential equation similar to the one used here. Thus the motion of the landslide is computed using a differential equation derived from first principles using Newtonian mechanics. However to expedite comparison with experiments, the landslide in [36] is considered moving on a straight inclined bottom with constant slope.

More recently, Khakimzyanov and Shokina [START_REF] Khakimzyanov | Numerical modelling of surface water waves arising due to a movement of the underwater landslide on an irregular bottom[END_REF], and Chubarov et. al. [START_REF] Chubarov | Numerical modelling of surface water waves arising due to movement of underwater landslide on irregular bottom slope[END_REF] have also used a differential equation to find the bottom motion. One major novelty in their work is that the landslide motion is computed on a bottom with an arbitrary shape. The timedependent bathymetry is then used to drive a numerical solver of the shallow-water equations. An advantage of this approach when compared to [36] is the reduced computation time. On the other hand, the description of the wave motion in the shallow-water theory is only approximate, and in particular, one important effect of surface waves, namely the influence of frequency dispersion is neglected.

The main aim of the current work is to study the dispersive wave generation in a closed basin [7] using a more realistic landslide model [START_REF] Chubarov | Simulation of surface waves by an underwater landslide[END_REF] while keeping the simplicity of the shallow-water approach. To this end, we use the so-called Peregrine system which is a particular case of a general class of model systems which arise in the Boussinesq scaling [START_REF] Boussinesq | Théorie de l'intumescence liquide appelée onde solitaire ou de translation se propageant dans un canal rectangulaire[END_REF]. A common feature of all Boussinesq-type systems is that they allow a simplified study of surface waves in which both nonlinear and dispersive effects are taken into account. In the present case, we need to use a Boussinesq system which can handle complex and timedependent bottom topography. Such a system was derived by Wu [START_REF] Wu | Generation of upstream advancing solitons by moving disturbances[END_REF], and can be used in connection with the dynamic bathymetry. An example of the type of situation considered here is shown in Figure 1, which shows how the bathymetry is given by the combination of a solid bottom, and a landslide profile sliding along the fixed bottom.

We conduct two main experiments. First, a comparison with the shallow-water theory is carried out. Second, the dependence of the tsunami characteristics on the initial depth of the landslide is investigated. The main findings of the present work are that the predictions of the shallow-water and Boussinesq theory are divergent for the cases treated in this paper, and that the effect of a finite fluid domain, such as a river, lake or fjord [START_REF] Poncet | A study of the tsunami effects of two landslides in the St. Lawrence estuary[END_REF] can lead to significantly different behavior when compared to tsunamis on an open ocean (see also [7]). The Boussinesq model in this paper is based on the assumption of an inviscid fluid, and irrotational flow. These are standard assumptions in the study of surface waves, and generally give good results, unless there are strong background currents in the fluid. Another effect which is not taken account of here is the wave resistance on the landslide due to waves created by the motion of the landslide. However, as observed in [START_REF] Harbitz | Mechanisms of tsunami generation by submarine landslides -a short review[END_REF], this effect is negligible for most realistic cases of underwater landslides. Viscosity is included in the dynamic model for the landslide as will be shown in the next section. In order to capture the effect of slide deformation during the evolution, a damping term in the equation of motion is included to model the internal friction in the landslide mass.

The paper is organised in the following way. In Section 2, the equation of motion for the landslide is developed. Then in Section 3, the Boussinesq model is recalled. In Section 4, solitary-wave solutions of the Peregrine system are found numerically. In Section 5, the numerical scheme for the Boussinesq system is explained and the numerical method is tested using the exact solutions of Section 4. Section 6 contains results of numerical runs for a few specific cases of bottom bathymetry, a parameter study of wave run-up in relation to the initial depth of the landslide, and a comparison with the shallow-water theory.

The landslide model

In this section we briefly present a mathematical model of underwater landslide motion. This process has to be addressed carefully since it determines the subsequent formation of water waves at the free surface. In the present study, we will assume the movable mass to be a solid body with a prescribed shape and known physical properties. Since the landslide mass and volume is preserved during the evolution, it is sufficient to determine the position of the barycenter x = x c (t) as a proxy for the motion of the whole body. As observed in the introduction, most studies of wave generation due to underwater landslides are based on prescribed bottom motion, or on solving the equation of motion on a uniform slope while taking account of different types of friction and viscous terms. Examples of such works are [START_REF] Di Risio | Three-dimensional experiments on landslide generated waves at a sloping coast[END_REF][START_REF] Pelinovsky | Simplified model of tsunami generation by submarine landslides[END_REF][START_REF] Watts | Comparing model simulations of three benchmark tsunami generation cases[END_REF]. A more general approach was recently pioneered by Khakimzyanov and Shokina in [START_REF] Khakimzyanov | Numerical modelling of surface water waves arising due to a movement of the underwater landslide on an irregular bottom[END_REF], where curvature effects of the bottom topography were taken into account. Since this model is applicable to a wider range of cases, we follow the approach of [START_REF] Khakimzyanov | Numerical modelling of surface water waves arising due to a movement of the underwater landslide on an irregular bottom[END_REF]. However, in addition to the effects included by Khakimzyanov and Shokina, our model also incorporates the effect of internal friction in the slide material (given by the dissipative force F i ) and the action of bottom friction, given by F b .

The static bathymetry is prescribed by a sufficiently smooth single-valued function z = -h 0 (x), and the landslide shape is initially prescribed by a localised function z = ζ 0 (x). To be specific, in this study we choose the following shape function for the landslide mass:

ζ 0 (x) = A 1 2 1 + cos( 2π(x-x 0 ) ℓ ) , |x -x 0 | ≤ ℓ 2 0, |x -x 0 | > ℓ 2 .
(2.1)

In this formula, A is the maximum height, ℓ is the length of the slide and x 0 is the initial position of its barycenter. It is clear that the model description given below and the method of numerical integration used in the present work is applicable to any other smooth profile, as long as it is sufficiently localized and fully submerged.

Since the landslide motion is translational, its shape at time t is given by the function z = ζ(x, t) = ζ 0 (xx c (t)). Recall that the landslide center is located at a point with abscissa x = x c (t). Then, the impermeable bottom for the water wave problem can be easily determined at any time by simply superposing the static and dynamic components. Thus the bottom boundary conditions for the fluid are to be imposed at

z = -h(x, t) = -h 0 (x) + ζ(x, t).
To simplify the subsequent presentation, we introduce the classical arc-length parameterisation, where the parameter s = s(x) is given by the formula

s = L(x) = x x 0 1 + (h ′ 0 (ξ)) 2 dξ. (2.
2)

The function L(x) is monotone and can be efficiently inverted to yield the original Cartesian abscissa x = L -1 (s). Within the parameterisation (2.2), the center of the landslide is initially located at a point with the curvilinear coordinate s = 0. The local tangential direction is denoted by τ and the normal direction by n.

A straightforward application of Newton's second law reveals that the landslide motion is governed by the differential equation

m d 2 s dt 2 = F τ (t)
, where m is the landslide mass and F τ (t) is the tangential component of the sum of forces acting on the moving submerged body. In order to project the forces onto the axes of the local coordinate system, the angle θ(x) between τ and Ox is needed. This angle is determined by θ(x) =arctan h ′ 0 (x) . Let us denote by ρ w and ρ ℓ the densities of the water and landslide material correspondingly. If V is the volume of the slide, then the total mass m is given by the expression

m := ρ ℓ + c w ρ w V, (2.3)
where c w is the added mass coefficient. As explained in [START_REF] Batchelor | An introduction to fluid dynamics, volume 61 of Cambridge mathematical library[END_REF], a portion of the water mass has to be added to the mass of the landslide since it is entrained by the underwater body motion. For a cylinder, the coefficient c w is equal exactly to one, but in the present case, the coefficient has to be estimated. The volume of the sliding material is given by V = W • S, where W is the landslide width in the transverse direction, and S can be computed by

S = R ζ 0 (x) dx.
The last integral can be computed exactly for the particular choice (2.1) of the landslide shape to give

V = 1 2 ℓAW.
The total projected force acting on the landslide can be conventionally represented as a sum of the force F g representing the joint action of gravity and buoyancy, and the total contribution of various dissipative forces.

The gravity and buoyancy forces act in opposite directions and their horizontal projection F g can be easily computed by

F g (t) = (ρ ℓ -ρ w )W g R ζ(x, t) sin θ(x) dx.
Now, let us specify the dissipative forces. The water resistance to the motion of the landslide F r due to viscous dissipation is proportional to the maximal transverse section of the moving body and to the square of its velocity. In addition, the coefficient sign ds is needed to dissipate the landslide kinetic energy independently of its direction of motion. Thus the force F r takes the form

F r = -sign ds dt 1 2 c d ρ w AW ds dt 2 ,
where c d is the resistance coefficient of the water. The friction force F f is proportional to the normal force exerted on the body due to the weight:

F f = -c f sign ds dt N(x, t).
The normal force N(x, t) is composed of the normal components of gravity and buoyancy forces, but also of the centripetal force due to the variation of the bottom slope:

N(x, t) = (ρ ℓ -ρ w )gW R ζ(x, t) cos θ(x) dx + ρ ℓ W R ζ(x, t)κ(x) ds dt 2 dx.
Here κ(x) is the signed curvature of the bottom which can be computed using the formula

κ(x) = h ′′ 0 (x) 1 + (h ′ 0 (x)) 2 3 2
.

We note that the last term vanishes for a plane bottom since κ(x) ≡ 0 in this particular case. Energy loss inside the sliding material due to internal friction is modeled by where c v is an internal friction coefficient. Finally, dissipation in the boundary layer between the landslide and the solid bottom is taken account of by the term

F i = -c v ρ ℓ W S ds dt ,
F b = -c b ρ w W ℓ ds dt ds dt ,
where c b is the Chézy coefficient. Finally, if we sum up the contributions of all the forces described above, we obtain the second order differential equation

(γ + c w )S d 2 s dt 2 = (γ -1)g I 1 (t) -c f σ(t)I 2 (t) -σ(t) c f γI 3 (t) + 1 2 c d A ds dt 2 -c v γS ds dt -c b ℓ ds dt ds dt , (2.4)
where γ := ρ ℓ ρw > 1 is the ratio of densities, σ(t) := sign ds dt and the integrals I 1,2,3 (t) are defined by In order to obtain a well-posed initial value problem, equation (2.4) has to be supplemented with initial conditions for s(0) and s ′ (0). In the remainder we always take homogeneous initial conditions, and consider the motion driven only by the gravitational acceleration of the landslide. However, different boundary conditions might also be reasonable from a modeling point of view.

I 1 (t) = R ζ(x, t) sin θ(x) dx, I 2 (t) = R ζ(x, t) cos θ(x) dx, I 3 (t) = R ζ(x, t)κ(x) dx.
In order to approximate solutions of equation (2.4), we employ the Bogacki-Shampine third-order Runge-Kutta scheme. The integrals I 1,2,3 (t) are computed using the trapezoidal rule, and once the landslide trajectory s = s(t) is found, we use equation (2.2) to find its motion x = x c (t) in the initial Cartesian coordinate system. This yields the bottom motion that drives the fluid solver.

For illustrative purposes we show a few examples of landslide trajectories over the bottom profile depicted in Figure 1. The other parameters used in the simulations are given in Section 6 and also in Table 2. We performed a series of simulations in order to study the effect of various dissipative terms on the landslide trajectory. The dependence on the friction coefficient c f is shown in Figure 2 where the landslide barycenter position x c (t) and its velocity v c (t) are shown as functions of time for c f = tan(1 • ), tan(2 • ) and tan(3 • ). n the case of the weak friction c f = tan(1 • ), the landslide reaches a sufficient speed to escape from the basin depicted on Figure 1. For the latter case (c f = tan(3 • )) we show also simultaneously the landslide speed v c (t) := dxx dt and its acceleration a c (t

) := dvc dt = d 2 xc dt 2
in Figure 3. In particular, one can see that the acceleration is a discontinuous function whose jumps correspond exactly to moments of time where the speed v c changes its sign, in accordance with the employed model (2.4). However, in our model there are also two new dissipative terms F i and F b whose importance has to be studied also. We fix the value of c f = tan(3 • ) for all subsequent experiments and we will vary only the two other coefficients c v and c b for fixed other parameters given in Table 2. These numerical results are presented in Figures 4 and5. One can see that the influence of these parameters on the landslide trajectory is weaker. However, we choose to keep them in the model in order to have more latitude for fine tuning the slide trajectory if need be.

The Boussinesq model

Once the motion of the landslide is determined, and therefore the time-dependent bathymetry h(x, t) = h 0 (x)ζ(x, t) is given, the next step is to consider the coupling between the bathymetry variations and the evolution of surface waves. The main assumptions on the fluid are that it is inviscid and incompressible, and that the flow is irrotational. Under these assumptions, the potential-flow free surface problem governs the motion of the fluid. However, in the present case, the fluid is shallow, and the waves at the surface are of small amplitude when compared to the depth of the fluid. In that case, the potential-flow problem may be simplified, and the model used in this paper is a variant of the so-called classical Boussinesq system derived by Boussinesq [START_REF] Boussinesq | Théorie de l'intumescence liquide appelée onde solitaire ou de translation se propageant dans un canal rectangulaire[END_REF].

Let us first consider the case of an even bottom, and a constant fluid depth d 0 . Denote a typical wave amplitude by a, and a typical wavelength by λ. The parameter α = a d 0 then describes the relative amplitude of the waves, and the parameter β = d 2 0 λ 2 measures the 'shallowness' of the fluid in comparison to the wavelength. In the case when both α and β are small and approximately of the same order of magnitude, the system

η t + d 0 u x + (ηu) x = 0, u t + gη x + uu x - d 2 0 3 u xxt = 0 (3.1)
may be used as an approximate model for the description of the evolution of the surface waves and the fluid flow. In (3.1), η denotes the deflection of the free surface from its rest position, and u denotes the horizontal fluid velocity at a height z = d 0 (-1 + 1/3) in the fluid column if z is measured from the rest position of the free surface. The same equation appears if the velocity is taken to be the average of the horizontal velocity over the flow depth.

The system (3.1) was first derived by Peregrine in [50], and falls into a general class of Boussinesq systems, as shown in the systematic studies [START_REF] Bona | Boussinesq equations and other systems for small-amplitude long waves in nonlinear dispersive media. I: Derivation and linear theory[END_REF][START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF]. As opposed to the shallowwater approximation, the pressure is not assumed to be hydrostatic, and the horizontal velocity varies with depth. In fact, the horizontal velocity profile is a quadratic function of z [START_REF] Whitham | Linear and nonlinear waves[END_REF]. Non-hysrostatic effects lead the appearance of linear dispersive terms in the governing equations. The problem of landslide generated waves has been addressed in the fully nonlinear shallow water framework [65, [START_REF] Chubarov | Simulation of surface waves by an underwater landslide[END_REF][START_REF] Yu | Modelling surface waves generated by a moving landslide with allowance for vertical flow structure[END_REF][START_REF] Beisel | Surface wave modeling generated by an underwater landslide moving along a nonuniform slope[END_REF]. Nevertheless, several authors obtained recently interesting results even in the linear [53,[START_REF] Seo | Edge waves generated by the landslide on a sloping beach[END_REF] or nonlinear [START_REF] Fernández-Nieto | A new Savage-Hutter type model for submarine avalanches and generated tsunami[END_REF][START_REF] Didenkulova | Tsunami waves generated by submarine landslides of variable volume: analytical solutions for a basin of variable depth[END_REF]7] hydrostatic models.

The derivation of (3.1) given in [50] also featured and extension to non-constant but time-independent bathymetry. However, the present case of a dynamic bottom profile calls for a system which allows for time-dependent bathymetry, and such a system was derived in [START_REF] Wu | Generation of upstream advancing solitons by moving disturbances[END_REF]. Given a bottom topography described by z = -h(x, t), the system takes the form

η t + (h + η)u x + h t = 0, u t + gη x + uu x = 1 2 h h t + (hu) x xt - h 2 6 u xxt . (3.2)
In order for this system to be asymptotically valid, we need α ∼ β as before. Moreover, concerning the unsteady bottom profile, we make the assumptions that h x ≤ O(αβ 1/2 ), and

h t ≤ O(αβ 1/2 ).
In comparison to the shallow-water equations with a time-dependent bottom topography, the system (3.2) has additional terms on the right-hand side of the second equation. The effect of these terms is to incorporate frequency dispersion into the model. One practical aspect of this modification is that wave breaking can be completely avoided as long as the amplitude of the waves is small enough. Wave breaking is also possible in evolution systems of Boussinesq type [START_REF] Bjørkavåg | Wave breaking in Boussinesq models for undular bores[END_REF][START_REF] Briganti | Boussinesq modelling of breaking waves: description of turbulence[END_REF], but the amplitudes occurring in the present problem are far from the breaking limit. The phase speed of a small-amplitude linear wave of wavelength 2π/k in the equation (3.2) with a stationary even bottom has the form

c 2 = gd 0 1 + d 2 0 3 k 2 ,
while the phase speed is given by c 2 = gd 0 tanh(kd 0 ) kd 0 in the linearized full water wave problem. Thus one might argue that the dispersion in (3.2) is too strong in comparison with dispersion in realistic water waves. However, as discussed in [START_REF] Bjørkavåg | Wave breaking in Boussinesq models for undular bores[END_REF], the linear dispersion relation of (3.2) is still closer to the dispersion relation of the original water-wave problem than most other standard Boussinesq equations which feature even faster decay of the phase speed with increasing k.

Solitary waves

Before the numerical method for approximating solutions of (3.2) is presented, we digress for a moment, and explain how to find numerically exact solutions of the system (3.1). These solutions will later be used to test the implementation of the numerical procedure. Assuming the special form

η(x, t) = η(ξ), u(x, t) = u(ξ), ξ := x -c s t,
and substituting this representation into the governing equations (3.1), there appears

-c s η ′ + (d + η)u ′ = 0, -c s u ′ + 1 2 (u 2 ) ′ + gη ′ + c s d 2 3 u ′′′ = 0.
Assuming decay of both η and u to zero as |x| → ∞, the integration of the mass conservation equation from -∞ to ξ gives the following relation between η and u:

u = c s η d + η , η = d • u c s -u . (4.1)
The momentum balance equation can now be integrated to yield

-c s u - d 2 3 u ′′ + 1 2 u 2 + gη = 0. (4.2)
Finally, in order to obtain a closed form equation in terms of the velocity u, we substitute the expression (4.1) for η into (4.2). The resulting differential equation can be written in operator notation as

Lu = N (u),
where the linear operator L and the nonlinear operator N , are defined respectively by

Lu := c s u - d 2 3 u ′′ , and N (u) := 1 2 u 2 + gdu c s -u .
While nothing formal appears to be known about existence of localised solutions of (4.1), (4.2), it is straightforward to compute approximations of solitary waves numerically. In particular, one may use the well known Petviashvili iteration method which takes the form The exponent q is usually defined as a function of the degree p of the nonlinearity, with the rule of thumb that the expression q := p p-1 generally works well. In our case, the nonlinearities are quadratic, so that we choose p = 2, and hence q = 1.

u n+1 = L -1 • N (u n ) • (u n , N (u n )) (u n , Lu n ) -q . ( 4.3)
The Petviashvili method was analyzed in [48], and can be very efficiently implemented using the Fast Fourier Transform [START_REF] Frigo | The Design and Implementation of FFTW3[END_REF]. The iteration can be started for instance with the third-order asymptotic solution of Grimshaw [START_REF] Grimshaw | The solitary wave in water of variable depth. Part 2[END_REF]. The iterative procedure is continued until the L ∞ norm between two successive iteration is on the order of machine precision. Figure 6 shows approximate solitary-wave solutions of (3.1) with various wave speeds, and compares them to the third-order asymptotic approximation of solitary-wave solutions of the full water-wave problem obtained by Grimshaw [START_REF] Grimshaw | The solitary wave in water of variable depth. Part 2[END_REF]. The left panel shows comparisons of the free-surface excursion, while the right panel shows a comparison of the horizontal component of the velocity field, evaluated at the non-dimensional height z given by z = -1 + 1/3. Figure 7 shows a comparison of the wavespeed-amplitude relation between the solitary-wave approximation of (3.1) and the ninth-order asymptotic approximation to the full water-wave problem obtained by Fenton [START_REF] Fenton | A ninth-order solution for the solitary wave[END_REF].

The numerical scheme

For the numerical discretisation, a finite-volume discretisation procedure similar to the one used in [4, 5] is employed. Let us take as a unit of length the undisturbed depth d 0 of the fluid above the barycenter of the landslide, and as a unit of time the ratio d 0 g . Then the Peregrine system (3.2) is rewritten in terms of the total water depth H as

H t + [Hu] x = 0, (5.1)
u t + 1 2 u 2 + (H -h) x = 1 2 hh xtt + 1 2 h(hu) xxt - 1 6 h 2 u xxt , (5.2)
The system (5.1), (5.2) can be formally rewritten in the form

V t + [ F(V) ] x = S b + M(V), (5.3)
where the density V and the advective flux F(V) are defined by

V ≡ H u , F(V) ≡ H u 1 2 u 2 + (H -h) .
The source term is defined by

S b ≡ 0 1 2 hh xtt ,
and the dispersive term is defined by

M(V) ≡ 0 1 2 h(hu) xxt -1 6 h 2 u xxt .
We begin our presentation by a discretisation of the hyperbolic part of (5.1), (5.2), which is the classical nonlinear shallow-water system, and then discuss the treatment of dispersive terms. The Jacobian of the advective flux F(V) is easily computed to be

A(V) = ∂ F(V) ∂V = u H 1 u ,
and it is clear that A(V) has the two distinct eigenvalues

λ ± = u ± c s , c s ≡ √ H.
The corresponding right and left eigenvectors are the columns of the matrices

R = H -H c s c s , L = R -1 = 1 2
H -1 c -1 s -H -1 c -1 s .
We consider a partition of the real line R into cells (or finite volumes)

C i = [x i-1 2 , x i+ 1 2 ] with cell centers x i = 1 2 (x i-1 2 + x i+ 1 2
) (i ∈ Z). Let ∆x i denote the length of the cell C i . In the sequel we will consider only uniform partitions with ∆x i = ∆x, ∀i ∈ Z. We would like to approximate the solution V(x, t) by discrete values. In order to do so, we introduce the cell average of V on the cell C i (denoted with an overbar), i.e.,

Vi (t) ≡ Hi (t) , ūi (t) = 1 ∆x C i V(x, t) dx.
A simple integration of (5.3) over the cell C i leads to the exact relation

d V dt + 1 ∆x F(V(x i+ 1 2 , t)) -F(V(x i-1 2 , t)) = 1 ∆x C i S b (V) dx ≡ Si .
Since the discrete solution is discontinuous at cell interfaces x i+ 1 2 (i ∈ Z), we replace the flux at the cell faces by the so-called numerical flux function

F(V(x i± 1 2 , t)) ≈ F i± 1 2 ( VL i± 1 2 , VR i± 1 2 ),
where VL,R i± 1 2 denotes the reconstructions of the conservative variables V from left and right sides of each cell interface (the reconstruction procedure employed in the present study will be described below). Consequently, the semi-discrete scheme takes the form

d Vi dt + 1 ∆x F i+ 1 2 -F i-1 2 = Si . (5.4)
In order to discretise the advective flux F(V), we follow the method of [START_REF] Ghidaglia | Une méthode volumes-finis à flux caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de conservation[END_REF][START_REF] Ghidaglia | On the numerical solution to two fluid models via cell centered finite volume method[END_REF] and use the following FVCF scheme

F(V, W) = F(V) + F(W) 2 -U(V, W) • F(W) -F(V) 2 .
The first part of the numerical flux is centered, the second part is the upwinding introduced through the Jacobian sign-matrix U(V, W) defined by

U(V, W) = sign A( 1 2 (V + W)) , sign(A) = R • diag(s + , s -) • L,
where s ± ≡ sign(λ ± ). After some simple algebraic computations, one can find

U = 1 2 s + + s - (H/c s ) (s + -s -) (c s /H) (s + -s -) s + + s - ,
the sign-matrix U being evaluated at the average state of left and right values.

Finally the source term S b (x, t) = (0, 1 2 hh xtt ), which is due to the moving bottom, is discretised by evaluating the bathymetry function and its derivatives at cell centers:

1 ∆x C i S b (x, t) dx ≈ 0, 1 2 h(x i , t) h xtt (x i , t) .
Recall that the bathymetry is composed of the static part and of the landslide subject to a translational motion:

h(x, t) = h 0 (x) -ζ(x, t) = h 0 (x) -ζ 0 x -x c (t) .
The derivative h xtt can be readily obtained from the formula

h xtt (x, t) = d 2 x c dt 2 d 2 ζ 0 dx 2 (x -x c (t)) - dx c dt 2 d 3 ζ 0 dx 3 (x -x c (t)).

High-order reconstruction

In order to obtain a higher-order scheme in space, we need to replace the piecewise constant data by a piecewise polynomial representation. This goal is achieved by various so-called reconstruction procedures such as MUSCL TVD [START_REF] Kolgan | Finite-difference schemes for computation of three dimensional solutions of gas dynamics and calculation of a flow over a body under an angle of attack[END_REF]62,[START_REF] Van Leer | Upwind and High-Resolution Methods for Compressible Flow: From Donor Cell to Residual-Distribution Schemes[END_REF], UNO [41], ENO [START_REF] Harten | ENO schemes with subcell resolution[END_REF], WENO [START_REF] Xing | High order finite difference WENO schemes with the exact conservation property for the shallow water equations[END_REF] and many others. In recent studies on unidirectional wave models [24] and on Boussinesq-type equations [START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF], the UNO2 scheme showed a good performance with small dissipation in realistic propagation and run-up simulations. Consequently, we retain this scheme for the discretisation of the advective flux of the Peregrine system (5.1), (5.2).

The main idea of the UNO2 scheme is to construct a non-oscillatory piecewise-parabolic interpolant Q(x) to a piecewise smooth function V(x) (see [41] for more details). On each segment containing the face x i+ 1 2 ∈ [x i , x i+1 ], the function Q(x) = q i+ 1 2 (x) is locally a quadratic polynomial and wherever v(x) is smooth we have

Q(x) -V(x) = 0 + O(∆x 3 ), dQ dx (x ± 0) - dV dx = 0 + O(∆x 2 ).
Also, Q(x) should be non-oscillatory in the sense that the number of its local extrema does not exceed that of V(x). Since q i+ 1 2 (x i ) = Vi and q i+ 1 2 (x i+1 ) = Vi+1 , it can be written in the form

q i+ 1 2 (x) = Vi + d i+ 1 2 {V} × x -x i ∆x + 1 2 D i+ 1 2 {V} × (x -x i )(x -x i+1 ) ∆x 2 ,
where

d i+ 1 2 {V} ≡ Vi+1 -Vi and D i+ 1 2
V is closely related to the second derivative of the interpolant since D i+ 1 2 {V} = ∆x 2 q ′′ i+ 1 2 (x). The polynomial q i+ 1 2 (x) is chosen to be the least oscillatory between two candidates interpolating V(x) at (x i-1 , x i , x i+1 ) and (x i , x i+1 , x i+2 ). This requirement leads to the following choice of

D i+ 1 2 {V} ≡ minmod D i {V}, D i+1 {V} with D i {V} = Vi+1 -2 Vi + Vi-1 , D i+1 {V} = Vi+2 -2 Vi+1 + Vi ,
and where minmod(x, y) is the usual minmod function defined as minmod(x, y) ≡ 1 2 [ sign(x) + sign(y) ] × min(|x|, |y|). To achieve the second order O(∆x 2 ) accuracy, it is sufficient to consider piecewise linear reconstructions in each cell. Let L(x) denote this approximately reconstructed function which can be written in this form

L(x) = Vi + S i • x -x i ∆x , x ∈ [x i-1 2 , x i+ 1 2 ].
In order to L(x) be a non-oscillatory approximation, we use the parabolic interpolation Q(x) constructed below to estimate the slopes S i within each cell

S i = ∆x × minmod dQ dx (x i -0), dQ dx (x i + 0) .
In other words, the solution is reconstructed on the cells while the solution gradient is estimated on the dual mesh as it is often performed in more modern schemes [4,[START_REF] Barth | Fundamentals, chapter Finite[END_REF]. A brief summary of the UNO2 reconstruction can be also found in [START_REF] Dutykh | Finite volume schemes for dispersive wave propagation and runup[END_REF]24].

Treatment of the dispersive terms

In this section, we explain how we treat numerically the dispersive terms of the Peregrine system (5.1), (5.2) which are present only in the momentum conservation equation (5.2). We propose the following approximation for the second component of M( V) of M( V):

M i ( V) = 1 2 hi hi+1 (ū t ) i+1 -2 hi (ū t ) i + hi-1 (ū t ) i-1 ∆x 2 - 1 6 h2 i (ū t ) i+1 -2(ū t ) i + (ū t ) i-1 ∆x 2 = hi 2∆x 2 hi-1 - 1 3 hi (ū t ) i-1 - 2 3∆x 2 h2 i (ū t ) i + hi 2∆x 2 hi+1 - 1 3 hi (ū t ) i+1 .
Note that this spatial discretisation is of the second order O(∆x 2 ) so as to be consistent with the UNO2 advective flux discretisation presented above. If we denote by I the identity matrix, we can now rewrite the semi-discrete scheme in the form

d H dt + 1 ∆x F (1) + ( V) -F (1)
-( V) = 0, (I -M) • dū dt + 1 ∆x F (2) + ( V) -F (2)
-( V) = S (2) b ,
where

F (1,2)
± ( V) are the two components of the advective numerical flux vector F at the right (+) and left (-) faces correspondingly, and S In order to advance the numerical solution forward in time one has to invert the matrix (I -M) at every time step. This is no problem in practice, since the matrix appears to be well conditioned in all cases we have considered. In fact, the invertibility of the matrix (I -M) can be rigorously shown to hold for small enough ∆x since the matrix is then diagonally dominant. The criterion for diagonal dominance in the present case is seen to be

1 + 2 3∆x 2 h2 i > - 1 6 h2 i ∆x 2 + hi 2∆x 2 hi-1 + - 1 6 h2 i ∆x 2 + hi 2∆x 2 hi+1 .
Using a Taylor expansion to express the terms hi-1 and hi+1 as hi-1 = hi -∆x h′ (x i ) + O(∆x 2 ) and hi+1 = hi + ∆x h′ (x i ) + O(∆x 2 ) , respectively, the criterion reduces to

1 + 1 3∆x 2 h2 i > hi h′ i ∆x + O(1),
and this is guaranteed to hold for small enough ∆x.

Time stepping

We assume that the linear system of equations is already inverted and we have the following system of ODEs:

V t = N (V, t), V(0) = V 0 .
In order to solve numerically the last system of equations, we apply the Bogacki-Shampine method proposed by Przemyslaw Bogacki and Lawrence F. Shampine in 1989 [START_REF] Bogacki | A 3(2) pair of Runge-Kutta formulas[END_REF]. It is a Runge-Kutta scheme of the third order with four stages. It has an embedded second order method which is used to estimate the local error and thus, to adapt the time step size. Moreover, the Bogacki-Shampine method enjoys the First Same As Last (FSAL) property so that it needs approximately three function evaluations per step. This method is also implemented in the ode23 function in Matlab [START_REF] Shampine | The MATLAB ODE Suite[END_REF]. The one step of the Bogacki-Shampine method is given by:

k 1 = N (V (n) , t n ), k 2 = N (V (n) + 1 2 ∆t n k 1 , t n + 1 2 ∆t), k 3 = N (V (n) ) + 3 4 ∆t n k 2 , t n + 3 4 ∆t), V (n+1) = V (n) + ∆t n 2 9 k 1 + 1 3 k 2 + 4 9 k 3 , k 4 = N (V (n+1) , t n + ∆t n ), V (n+1) 2 = V (n) + ∆t n 4 24 k 1 + 1 4 k 2 + 1 3 k 3 + 1 8 k 4 .
Here V (n) ≈ V(t n ), ∆t is the time step and

V (n+1) 2
is a second order approximation to the solution V(t n+1 ), so the difference between V (n+1) and V (n+1) 2

gives an estimation of the local error. The FSAL property consists in the fact that k 4 is equal to k 1 in the next time step, thus saving one function evaluation.

If the new time step ∆t n+1 is given by ∆t n+1 = ρ n ∆t n , then according to the H211b digital filter approach [57,[START_REF] Söderlind | Adaptive time-stepping and computational stability[END_REF], the proportionality factor ρ n is given by

ρ n = δ ǫ n β 1 δ ǫ n-1 β 2 ρ -α n-1 , (5.5)
where ǫ n is a local error estimation at time step t n , and the constants β 1 , β 2 and α are defined by

α = 1 4 , β 1 = 1 4p , β 2 = 1 4p .
The parameter p gives the order of the scheme, and p = 3 in our case. It appears that the second-order convergence is achieved.

Remark 1. The adaptive strategy (5.5) can be further improved if we smooth the factor ρ n before computing the next time step ∆t n+1 :

∆t n+1 = ρn ∆t n , ρn = ω(ρ n ).
The function ω(ρ) is called the time step limiter and should be smooth, monotonically increasing and should satisfy the following conditions:

ω(0) < 1, ω(+∞) > 1, ω(1) = 1, ω ′ (1) = 1.
One possible choice was suggested in [START_REF] Söderlind | Adaptive time-stepping and computational stability[END_REF]:

ω(ρ) = 1 + κ arctan ρ -1 κ .
In our computations the parameter κ is set to 1.

Validation

The scheme described in this section is implemented in MATLAB, and runs on a workstation. To check whether the implementation is correct, we use the approximate solitary waves of (3.1), computed in the last section. These are used as initial data in the fully discrete scheme, and integrated forward in time. The computed solutions are then compared to the same solitary waves, but shifted forward in space by ct 0 , where c is the wave speed, and t 0 is the final time. This procedure is repeated a number of times with different spatial gridsizes. As a result, it is possible to find the spatial convergence rate of the scheme. As is visible in Figure 8, the convergence achieved by the practical implementation of the discretisation is very close to the theoretical convergence rate. Since the temporal discretisation is adaptive, we do not present a convergence study in terms of the timestep ∆t.

Wave generation by moving bottom

We have just shown the convergence of our scheme under the mesh refinement. Even if the solution we used in validation is fully nonlinear, it only exists on the flat bottom. Since in the present study we are mainly interested in the wave generation by bottom motion, the next validation test will be entirely devoted to this question. Namely, we are going to use an analytical solution to the linearized full Euler equations also known as the Cauchy-Poisson problem. The use of this solution in tsunami generation problems was first proposed by J. [START_REF] Hammack | A note on tsunamis: their generation and propagation in an ocean of uniform depth[END_REF] [START_REF] Hammack | A note on tsunamis: their generation and propagation in an ocean of uniform depth[END_REF].

We consider the linearized water wave problem for a fluid layer of uniform depth z = -d 0 = const. However, a portion of the bottom can move vertically and its deformation is given by a smooth function ζ(x, t) such that ζ(x, 0) ≡ 0 at time t the bottom profile is given by z = -d 0 + ζ(x, t). Moreover, we will make a special assumption about the structure of the bottom deformation:

ζ(x, t) = T (t)ζ 0 (x), T (t)
:= 1 -e -αt , α > 0, t ≥ 0.
Obviously, we have to assume that ||ζ 0 || ≪ 1 so that the linear approximation be valid. Then, the free surface elevation at any time is given by the following formula [START_REF] Hammack | A note on tsunamis: their generation and propagation in an ocean of uniform depth[END_REF][START_REF] Dutykh | Linear theory of wave generation by a moving bottom[END_REF]:

η(x, t) = - α 2 2π R ζ0 (k) cosh(kd 0 ) • e -αt -cos(ωt) -ω α sin(ωt) α 2 + ω 2 • e -ikx dk, (5.6)
where ζ0 (k) is the Fourier transform of ζ 0 (x) and ω := gk tanh(kd 0 ) is the wave frequency corresponding to the wavenumber k. The above integral can be easily computed using the FFT algorithm. To fix the ideas for numerical computations, we will take the following localized oscillatory bottom deformation:

ζ 0 (x) = a cos(k 0 x)e -λ 0 x 2 , λ 0 > 0.
The values of all parameters used in numerical simulation are given in Table 1. The nonlinearity parameter a/d 0 is chosen to be 0.05, which is far above the nonlinearity of the earthquake generated tsunamis. However, we think that this value corresponds better to the scope of the present study. In order to simulate this set-up using the Peregrine system, we consider a symmetric 1D computational domain [-220, 220] discretized into N = 2000 equal control volumes. The time stepping tolerance parameter was set far below the spatial discretization error (∼ O(∆x 2 )). First, we will take a moderately fast bottom uplift corresponding to the parameter α = 1.0. Computational results are presented in Figures 9(a-e). One can see that the overall agreement is fairly good even if some small differences can be noticed on Figures 9(c-d). However, the resulting wave form predicted by the Peregrine system follows closely the linearized full Euler solution (5.6). Now, we will double the bottom uplift speed (α = 2.0). This result is presented in Figures 10(a-e).

One can see more substantial differences during the generation phase (see panels b-d). However, here again the resulting wave is surprisingly well represented by the Boussinesqtype equations. The observed discrepancies during the generation phase are essentially due to the simplified structure of the vertical speed in Boussinesq-type equations [START_REF] Dutykh | On the use of the finite fault solution for tsunami generation problems[END_REF].

Numerical results and discussion

Let us consider a one-dimensional computational domain I = [a, b] = [0, 220] composed of two regions: the generation region and a sloping beach on the right. More specifically, the static bathymetry function h 0 (x) is given by a smoothed out profile generated from the expression

h 0 (x) = d 0 + tan δ • (x -a) + p(x), a ≤ x ≤ m, d 0 + tan δ • (m -a) -tan δ • (x -m), x > m,
where the function p(x) is defined as

p(x) = A 1 sech (k 1 (x -x 1 )) + A 2 sech (k 2 (x -x 2 )),
In essence, this function represents a perturbation of the sloping bottom by two underwater bumps. We made this nontrivial choice in order to illustrate the advantages of our landslide model, which was designed to handle general non-flat bathymetries. The parameters can be chosen in order to fit a given bathymetry, but the particular values used here are A 1 = 4.75, A 2 = 8.85, k 1 = 0.06, k 2 = 0.13, x 1 = 45, x 2 = 80, and m = 120. The bottom profile for these parameters is depicted in Figure 11. Of course, in general, if the bottom topography is known, then a numerical bathymetry map could also be used.

We now present some results of the solution of the surface wave problem using the model in Section 3, integrated numerically with the method of Section 5. A landslide is introduced on the left side of the bathymetry, and using the method of Section 2, its path along the bottom is determined by following the barycenter. Simultaneously, the system (3.2) is solved with the time-dependent bottom topography given from the solution of the landslide problem. The problem is integrated up to a final time T . Figure 12 shows wave records at six virtual wave gauges for both the dispersive system (3.2) and the shallowwater system. It appears from this figure that the shallow-water system underpredicts the development of free-surface oscillations. In particular, the wave gauges located at x = 40 and x = 60 show similar waveheights for both the shallow-water, and the dispersive system, but a qualitative divergence, as small oscillations are already developing which are not captured by the shallow-water system. Once the waves have propagated to the wave gauges located at x = 80, the dispersive oscillations have amplified, so that the waveheight is larger by a factor of 2 to 3 than the waveheight predicted by the shallow-water system. Going further to the wave gauges located at x = 100 and x = 120, the now rising bottom starts to have a damping effect on the waves.

The maximum and minimum free surface elevation over the whole domain is shown in Figure 13. On the lower panel of the same Figure 13 we show the maximal unsigned horizontal velocity. One can see that for short times the hydrostatic and dispersive models give very close extreme values. Later the differences start to appear due to the accumulation of dispersive effects.

Figure 14 shows the development of the kinetic energy of the landslide mass and simultaneously the total (kinetic plus potential) energy contained in the body of the fluid and the surface waves. Energy development is an important question in the study of tsunamis, and there have been studies exclusively devoted to this question [60]. Energy issues connected to water wave models of Boussinesq type have also been studied before [START_REF] Ali | Energy balance for undular bores[END_REF][START_REF] Ali | Mechanical balance laws for Boussinesq models of surface water waves[END_REF][START_REF] Dutykh | Energy of tsunami waves generated by bottom motion[END_REF]. While these models contained a source of energy, in the case at hand, the work done by friction as the landslide slides down the bottom acts as a drain of energy, and after the landslide has come to rest, all energy has been transferred to the fluid. However, not all energy can be considered as residing in the wave motion, because a significant amount of energy is needed to lift the water from the final position of the landslide to the initial position of the landslide. This results in a large increase in potential energy of the fluid, and only a fraction of the potential energy of the landslide is transferred to the wave motion. This fact has also been explained in previous works [START_REF] Harbitz | Mechanisms of tsunami generation by submarine landslides -a short review[END_REF].

In order to compute the wave energy in the fluid, we use the integral

E w = b a g 2 η 2 + 1 2 (h 0 + η)u 2 dx, (6.1)
which arises from the shallow-water theory. The kinetic energy of the landslide is given by

E sl = 1 2 mv 2 , (6.2)
with the generalised mass m given by (2.3), and v = ds dt as defined in Section 2. Figure 14 shows the development of the wave energy and kinetic energy of the landslide. The upper panel shows the energy according to the shallow-water and dispersive model. The lower panel shows the kinetic energy of the landslide.

We have also computed the Froude number Fr =

v √ gh(xc)
during the evolution. Here v is the x-component of the velocity of the barycenter of the landslide, x c is the position of the barycenter, and h(x c ) is the corresponding local water depth. This number was always found to be much less than 1 in all numerical experiments. The maximum value was generally about 0.5.

To compute the wave run-up and draw-down, we use exact representations given by Choi et al. [START_REF] Choi | Rapid forecasting of tsunami runup heights from 2-D numerical simulations[END_REF] (a similar formula was also derived in [START_REF] Didenkulova | Run-up of long waves on a beach: the influence of the incident wave form[END_REF]). On the right beach, the undisturbed water depth at the edge of the computational domain is h = 3, and the distance from the computational domain to the shore line is L = 30. Using the shallow-water wave speed, the Run-up on the left and right beach using (6.4), computed with the dispersive system (solid curve) and the nonlinear shallow-water system (dashed curve) as a function of (non-dimensional) time.

travel time of a wave from the edge of the computational domain to the shore is computed as

T = 2L √ gh = 2 L gα . (6.3)
Then the formula for the wave run-up R at the shore reads

R = t-T 0 t -τ (t -τ ) 2 -T 2 • dη dτ (x, τ ) dτ (6.4)
with x = 220. At the left beach, the undisturbed water depth is h = 1.642, and the distance to the beach is L = 11.2814. A similar formula can be then be computed for x = 0. Figure 15 shows the run-up on the left and right beaches both in the Boussinesq scaling and in the shallow-water theory. While the agreement is fair on the left beach, it appears immediately that the Boussinesq theory predicts a wave run-up on the right beach which is much larger (roughly by a factor of two) than the wave run-up according to the shallowwater theory. A possible explanation for this divergence is the nature of the numerical solver when applied to the shallow-water system. In this case, there is continuous numerical dissipation through the handling of hyperbolic wave breaking. Since the waves do not break in the Boussinesq scaling, the dissipation is not present, or at least much smaller. The difference can also be read off from the comparison of the wave energy in the Boussinesq and shallow-water system provided in Figure 14. It can be seen there that the wave energy in the shallow-water model starts to diverge from the Boussinesq model at nondimensional time t = 50. The difference between the two increases continuously, until at the final time, the Boussinesq energy is about 50% larger than the shallow-water energy. Note that significant run-up in Figure 15 does not happen until non-dimensional time t = 75, at which time the energy in the Boussinesq system is already much larger than in the shallow-water system.

In Figure 16, we have plotted the maximum wave amplitude, the minimum wave amplitude, and the maximum wave run-up on the left and right beaches. In comparison to previous studies, such as [36], where an open domain was used, it appears that in our case, the maximal amplitude, as well as the run-up have a minimum at d 0 between 1 and 1.5. In [36], it was found that maximum wave amplitude and run-up (on the left beach) were strictly decreasing functions of d 0 . The phenomenon of rising amplitude and run-up may be accredited to resonant effects which are absent on an open domain (such as an ocean beach), but cannot be neglected for tsunamis generated by landslides in rivers and lakes.

Conclusion

The influence of an underwater landslide on surface waves in a closed basin has been studied. The key features of the study have been that the motion of the underwater landslide is determined by integrating a second-order ordinary differential equation derived from first principles of Newtonian mechanics, and that the wave motion is studied in the Boussinesq scaling which allows for both nonlinear and dispersive effects. The dynamics of the motion of the bottom have been developed following recent work in [START_REF] Khakimzyanov | Numerical modelling of surface water waves arising due to a movement of the underwater landslide on an irregular bottom[END_REF]. The Boussinesq model which has been utilised here allows for a dynamic bathymetry, and was derived in [START_REF] Wu | Generation of upstream advancing solitons by moving disturbances[END_REF]. The numerical method used in this paper is an extension of the method put forward in [4,[START_REF] Barth | Fundamentals, chapter Finite[END_REF].

The results presented in Section 6 clearly show that dispersion may have a strong effect on the run-up and draw-down at the beaches, but it is not clear which of the two models (the shallow-water, or the Boussinesq model) paints a more realistic picture of the actual wave conditions. We have no way of quantifying the energy dissipation in the shallow-water solver other than comparing the total wave energy with that of the Boussinesq model. As the difference is rather large, we expect a significant amount of numerical dissipation in the shallow-water simulation. While the shallow-water model simulation might be closer to physical reality where actual damping occurs because of molecular viscosity and fluidstructure interactions, it is likely that the Boussinesq model exhibits a closer resemblance to the Euler equations which are taken as the basic governing equations in this work.

Of course, this difference could be more or less pronounced depending on the particular case under study. For example, the divergence between the shallow-water theory and the dispersive model is stronger at the right beach than at the left beach. The results also show that a finite domain exhibits different behavior than a half-open domain (such as used in [36]) with respect to the dependence of the wave run-up on the initial depth of the landslide. While the run-up is a strictly decreasing function of the initial depth in an open domain, a closed domain appears to exhibit resonant effects, which make the dependence more complex.

Figure 1 .

 1 Figure 1. The fixed bathymetry z = h 0 (x) and the position of the landslide after 50 s. The position of the barycenter is indicated by a black dot.

Figure 2 .

 2 Figure 2. Position and velocity of the barycenter of the landslide as functions of dimensional time for three different values of the friction coefficient c f .

Figure 3 .

 3 Figure 3. Velocity and acceleration of the barycenter of the landslide as a function of dimensional time. The friction coefficient is c f = tan(3 • ). The discontinuities in the acceleration are due to the coefficient sign ds dt

2 Figure 4 .

 24 Figure 4. Position and velocity of the barycenter of the landslide as functions of dimensional time for three different values of the friction coefficient c v .

Figure 5 .

 5 Figure 5. Position and velocity of the barycenter of the landslide as functions of dimensional time for three different values of the friction coefficient c b .

Figure 6 .

 6 Figure 6. Comparison of the numerical approximation of solitary wave solutions of (3.1) to Grimshaw's third-order asymptotic approximation of solitary waves using the Euler equations for the full water wave problem. The left panel shows the surface elevation, and the right panel shows the horizontal velocity at z = d 0 (-1 + 1/3).

Figure 7 .

 7 Figure 7. Amplitude-speed relation of solitary wave solutions of (3.1) and of Fenton's ninth-order asymptotic approximation of solitary waves using the Euler equations for the full water wave problem.

  of the second component of S b .

Figure 8 .

 8 Figure 8. Convergence rate of the finite-volume scheme in the L ∞ -norm (left panel) and the L 2 -norm (right panel). The numerical integration of a solitary wave as shown in Figure6is compared to a translated profile. It appears that the second-order convergence is achieved.

Figure 9 .Figure 10 .

 910 Figure 9. Free surface waves generated by a moderately fast bottom motion. The blue dashed line corresponds to the analytical Cauchy-Poisson solution, while the solid black line is our numerical solution to the Peregrine system.

Figure 11 .

 11 Figure 11. The physical setup of the problem. The river bed is indicated in dark grey. The computational fluid domain is shaded light grey, and the landslide is visible in black. Note the difference in horizontal and vertical scales in the left panel. The right panel shows a closeup of the left beach and part of the landslide in a one-to-one aspect ratio.

Figure 12 .Figure 13 .

 1213 Figure 12. Time series of the surface elevation at wave gauges located at x = 40, x = 60, x = 80, x = 100 and x = 120. The solid (blue) curve depicts the wave elevation computed with the dispersive system (3.2), and the dashed curve represents results obtained from the shallow-water system. All variables are non-dimensional.

Figure 14 .

 14 Figure 14. Development of the wave energy, and the kinetic energy of the landslide as a function of (non-dimensional) time. Note that the kinetic energy of the landslide starts from 0 (all energy is potential) and also ends at 0 (all energy has been dissipated or transferred to the fluid).

  Figure 15. Run-up on the left and right beach using (6.4), computed with the

  Maximal run-up on the right

Figure 16 .

 16 Figure 16. Maximal and minimal wave amplitude, and the maximum run-up on the left and right beaches as a function of the initial depth of the center of the landslide d 0 .

Table 1 .

 1 Values of various parameters used to simulate the wave generation by moving bottom.

	Gravity acceleration: g	1.0
	Undisturbed water depth: d 0	1.0
	Bottom displacement amplitude: a	0.05
	Bottom oscillation inverse length: k 0 Bottom localization parameter: λ 0 Vertical uplift speed: α	π 40 0.7 × 10 -3 1.0 and 2.0

Table 2 .

 2 Values of various parameters used in numerical computations.

Acknowledgements D. Dutykh would like to thank the University of Bergen for support and hospitality during the preparation of this manuscript. Support from the Agence Nationale de la Recherche under project ANR-08-BLAN-0301-01 (MathOcéan) and from ERC under project ERC-2011-AdG 290562-MULTIWAVE is also gratefully acknowledged. H. Kalisch acknowledges support of the Research Council of Norway through grant no. NFR 213474/F20.