Boussinesq modeling of surface waves due to underwater landslides

Consideration is given to the influence of an underwater landslide on waves at the surface of a shallow body of fluid. The equations of motion which govern the evolution of the barycenter of the landslide mass include various dissipative effects due to bottom friction, internal energy dissipation, and viscous drag. The surface waves are studied in the Boussinesq scaling, with time-dependent bathymetry. A numerical model for the Boussinesq equations is introduced which is able to handle time-dependent bottom topography, and the equations of motion for the landslide and surface waves are solved simultaneously. The numerical solver for the Boussinesq equations can also be restricted to implement a shallow-water solver, and the shallow-water and Boussinesq configurations are compared. A particular bathymetry is chosen to illustrate the general method, and it is found that the Boussinesq system predicts larger wave run-up than the shallow-water theory in the example treated in this paper. It is also found that the finite fluid domain has a significant impact on the behavior of the wave run-up.


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 [30,45] . 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 [51], and recent years have seen a multitude of works devoted to the study of such underwater landslides and the resulting effect on surface waves [3,12,14,21,29,30,39,40,44,51]. As suggested in [23], 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 drawdown 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 [17,49].
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 [30] 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 [30] relies on integrating the full water-wave equations using an irrotational boundary-element code, and using an open boundary with transmission conditions [28,27]. While most works have considered a given dynamics for the landslide, the bottom motion in [30] 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 [30] is considered moving on a straight inclined bottom with constant slope.
More recently, Khakimzyanov and Shokina [35], and Chubarov et. al. [12] 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 time-dependent bathymetry is then used to drive a numerical solver of the shallow-water equations. An advantage of this approach when compared to [30] 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 introduce the effect of dispersion into the description of the wave motion at the surface of the fluid while keeping the simplicity of the shallow-water approach. To this end, we utilise the so-called Peregrine system which is a particular case of a general class of model systems which arise in the Boussinesq scaling [10]. 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 time-dependent bottom topography. Such a system was derived by Wu [56], and can be used in connection with the dynamic bathymetry.
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 [44] can lead to significantly different behavior when compared to tsunamis on an open ocean.
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 [32], 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 [13,42,54]. A more general approach was recently pioneered by Khakimzyanov and Shokina in [35], 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 [35].
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: (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 reasonable shape.
Since the landslide motion is translational, its shape at time t is given by the function z = ζ(x, t) = ζ 0 (x − x 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 To simplify the subsequent presentation, we introduce the classical arc-length parameterisation, where the parameter s = s(x) is given by the formula 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 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 where c w is the added mass coefficient. As explained in [6], 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 The last integral can be computed exactly for the particular choice (2.1) of the landslide shape to give 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 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 dt is needed to dissipate the landslide kinetic energy independently of its direction of motion. Thus the force F r takes the form 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: 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:  Here κ(x) is the signed curvature of the bottom which can be computed using the formula 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 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  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.
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(t) in the initial Cartesian coordinate system. This yields the bottom motion that drives the fluid solver.

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 [10].
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 β = 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 [43], and falls into a general class of Boussinesq systems, as shown in the systematic studies [9,38]. As opposed to the shallow-water 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 [55].
The derivation of (3.1) given in [43] 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 [56]. Given a bottom topography described by z = −h(x, t), the system takes the form 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 [7], but the amplitudes occurring in the present problem are far from the breaking limit. One may argue that the dispersion in (3.2) is too strong in comparison with dispersion in realistic water waves. However, the linear dispersion relation of (3.2) with fixed even bottom is closer to the dispersion relation of the original water-wave problem than most other standard Boussinesq equations [7].

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 and substituting this representation into the governing equations (3.1), there appears 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: The momentum balance equation can now be integrated to yield 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 where the linear operator L and the nonlinear operator N , are defined respectively by 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.
The Petviashvili method was analyzed in [41], and can be very efficiently implemented using the Fast Fourier Transform [22]. The iteration can be started for instance with the third-order asymptotic solution of Grimshaw [31]. The iterative procedure is continued until the L ∞ norm between two successive iteration is on the order of machine precision. Figure 4 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 [31]. 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 heightz given byz = −1 + 1 3. Figure 5 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 [20].

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  The system (5.1),(5.2) can be formally rewritten in the form where the density V and the advective flux F(V) are defined by .
The source term is defined by , and the dispersive term is defined by 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 and it is clear that A(V) has the two distinct eigenvalues The corresponding right and left eigenvectors are the columns of the matrices

We consider a partition of the real line R into cells (or finite volumes)
. 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., A simple integration of (5.3) over the cell C i leads to the exact relation 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 denotes the reconstructions of the conservative variablesV 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 In order to discretise the advective flux F(V), we follow the method of [25,26] and use the following FVCF scheme 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 where s ± ≡ sign(λ ± ). After some simple algebraic computations, one can find 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: Recall that the bathymetry is composed of the static part and of the landslide subject to a translational motion: The derivative h xtt can be readily obtained from the formula

5.1.
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 [36,52,53], UNO [34], ENO [33], WENO [57] and many others. In recent studies on unidirectional wave models [19] and on Boussinesq-type equations [18], 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 piecewiseparabolic interpolant Q(x) to a piecewise smooth function V(x) (see [34] for more details). On each segment containing the face is locally a quadratic polynomial and wherever v(x) is smooth we have 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 and where minmod(x, y) is the usual minmod function defined as 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 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 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,5]. A brief summary of the UNO2 reconstruction can be also found in [18,19].

5.2.
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).
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 where F

Time stepping.
We assume that the linear system of equations is already inverted and we have the following system of ODEs: 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 [8]. 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 [46]. The one step of the Bogacki-Shampine method is given by:  Figure 4 is compared to a translated profile. It appears that the second-order convergence is achieved.
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 [47,48], the proportionality factor ρ n is given by where ǫ n is a local error estimation at time step t n , and the constants β 1 , β 2 and α are defined by The parameter p gives the order of the scheme, and p = 3 in our case.
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: One possible choice was suggested in [48]: In our computations the parameter κ is set to 1.

5.4.
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  6, 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.

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 where the function p(x) is defined as 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 7. 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 8 shows wave records at six virtual wave gauges for both the dispersive system (3.2) and the shallow-water 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 is shown in Figure 9. Figure 10 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 [50]. Energy issues connected to water wave models of Boussinesq type have also been studied before [1,2,16]. 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 [32]. In order to compute the wave energy in the fluid, we use the integral which arises from the shallow-water theory. The kinetic energy of the landslide is given by with the generalised mass m given by (2.3), and v = ds dt as defined in Section 2. Figure  10 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 F r = 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. [11] (a similar formula was also derived in [15]). 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 travel time of a wave from the edge of 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 11. This figure shows 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.
computational domain to the shore is computed as Then the formula for the wave run-up R at the shore reads 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 11 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 shallow-water 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 10. It can be seen there that the wave energy in the shallow-water model starts to diverge from the Boussinesq model at non-dimensional 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 11 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 12, 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 [30], 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 [30], 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 have been studied. The key features of the study have been that the motion of the underwater landslide have been determined by integrating a second-order ordinary differential equation derived from first principles of Newtonian mechanics, and that the wave motion has been 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 [35]. The Boussinesq model which has been utilised here allows for a dynamic bathymetry, and was derived in [56]. The numerical method used in this paper is an extension of the method put forward in [4,5]. 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. 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 [30]) 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.