Momentum Conservative Scheme for Simulating Wave Runup and Underwater Landslide

This paper focuses on the numerical modelling and simulation of tsunami waves triggered by an underwater landslide. The equation of motion for water waves is represented by the Nonlinear Shallow Water Equations (NSWE). Meanwhile, the motion of underwater landslide is modeled by incorporating a term for bottom motion into the NSWE. The model is solved numerically by using a finite volume method with a momentum conservative staggered grid scheme that is proposed by Stelling & Duinmeijer 2003 [12]. Here, we modify the scheme for the implementation of bottom motion. The accuracy of the implementation for representing wave runup and rundown is shown by performing the runup of harmonic wave as proposed by Carrier & Greenspan 1958 [2], and also solitary wave runup of Synolakis, 1986 [14], for both breaking and non-breaking cases. For the underwater landslide, result of the simulation is compared with simulation using the Boundary Integral Equation Model (BIEM) that is performed by Lynett and Liu, 2002 [9].


I. INTRODUCTION
n general, tsunami waves are commonly caused by an earthquake or an underwater volcano explosion.But in some tsunami cases, it can also be caused by an underwater landslide.Previously, topics on the tsunami that caused by underwater landslide get less attention, since the tsunami waves are commonly generated from seismic activities [3].In the last few decades, it has been known that the underwater landslide may trigger the generation of large scale of tsunami wave.In October 16, 1979, a tsunami has been generated by an underwater landslide in Nice, France.Although it is still debatable, the tsunami was not generated by seismic activity.Penilovsky et al [12] suspect that the cause of the tsunami is due to instability of sediment near the location.
An underwater landslide itself can be triggered by many factors.Instability and less solid sediment near shoreline are potential to trigger the generation of an underwater landslide [7].When wave approaches the sediment, an amount of wave energy is transferred to the sediment.The other way around, when the wave goes backward, the water draws the sediment seaward.As the event continues, the density of the sediment may decrease and prone to trigger generation of underwater landslide.In Indonesia, there are some tsunami events triggered by underwater landslide.An example is a tsunami is Flores in 1992 where an earthquake with 7.7 Richter scale generated tsunami wave up to 25 m and costing 2200 lives [11].The tsunami was suspected caused by an underwater landslide that is generated by an earthquake [15].
One method to understand the tsunami generation and propagation is to simulate the phenomenon numerically using an appropriate, accurate and efficient wave model and numerical scheme.Especially for modeling of tsunami generated by underwater landslide, effects of bottom motion should be incorporated in the fluid motion; hence an accurate numerical implementation is needed to achieve a good result.Moreover, the numerical scheme should able to simulate wave runup accurately.Research related with tsunami generated by underwater landslide, have been done in the last few decades.Hammack, in 1973 [6], performs physical experiments as well as theoretical modeling regarding wave generation by bottom motion.Lynett and Liu [9] simulate tsunami wave generated by underwater landslide by using a Boussinesq model that is implemented using high order finite difference scheme.Later, Fuhrman & Madsen [3] reproduce the simulation using their higher order Boussinesq model.This paper focuses on simulation of tsunami generation of underwater landslide.Rather than using a higher order Boussinesq models as in [9,4], here we use a rather simple wave model; the non-dispersive Nonlinear Shallow Water Equations (NSWE).To be able to simulate accurately wave runup and motion of bottom, the model is implemented by using finite volume method a momentum conservative staggered grid scheme that is first proposed by Stelling & Duinmeijer [13].A special treatment in the numerical implementation is made for describing the motion of bottom.The obtained numerical implementation is validated by comparing solitary wave runup experiment performed by Synolakis, 1986 [14] and analytical solution of harmonic wave runup obtained by Carrier and Greenspan [2].Moreover, we also reproduce simulation of tsunami wave generation by underwater landslide performed by [9,4].The organization of this paper is as follows.Section 2 describes the SWE model, is then followed by description of the implementation of staggered grid scheme in Section 3. In Section 4, the model is tested for simulation of solitary wave runup of Synolakis [14] and Carrier & Greenspan [2] and for simulation of tsunami wave generated by underwater landslide of Lynett & Liu [9].Finally we conclude the paper in the last section.

II. NONLINEAR SHALLOW WATER EQUATIONS
The Nonlinear Shallow Water Equations (NSWE) can describe wave dynamic especially for long wave such as tsunami.In the SWE, the vertical variation in fluid is assumed to be small; therefore it is valid only for long waves.Compared to the Boussinesq type of model (BTM) that is dispersive, the SWE is a non- dispersive model, and it has a simpler model compared to the BTM, results in a less computational cost than the BTM, hence it is preferred for simulating tsunami wave propagation and runup.
We restrict the problem to be 1 dimensional (1D) problem.Let x, z, and t be the horizontal, vertical coordinates and the time, respectively.The wave elevation and the horizontal velocity are denoted by !(!, !) and !!, ! as illustrated in Fig. 1.We define !(!, !) as the water depth measured from still water level (z=0), ℎ !, != !!, !+ !(!, !) as the total depth, and g as the gravity accelleration with != 9.81 !/! ! .The NSWE is given by the following equations Here, !!denotes the coefficient for the bottom roughness.The equation ( 1) and ( 2) are known as the continuity and the momentum equation, respectively.Note that when there is no bottom motion or !!!= 0, then eq. ( 1) can be simplified as

III. NUMERICAL IMPLEMENTATION
In this paper, the NSWE in equations ( 1) and ( 2) are implemented numerically by using the Finite Volume Method (FVM) with momentum conservative Staggered Grid (SG) scheme.Firstly, we consider a discretization in horizontal space with left and right boundary are denoted by !! and !!, respectively.Let !!∈ !, with != 1, 2, . .! !, be the grid location that is named as full-grid, and !!!!/! ∈ !, with != 1, 2, . .! !, be the grid location that is named as half-grid, as illustrated in Fig. 2. horizontal velocity is placed in the half-grid.By using first order forward time discretization and second order spatial discretization, the discretization for eq.( 3) is given as follows * is not defined in full grid as well as in half grid.Here, ℎ * is approximated by using upwind method as follows , jika 0 , jika 0 The condition (5) states that when the wave is moving to the right (or !!!!/! ≥ 0), ℎ * is assigned with ℎ ! or its value in the left side.The other way, when the wave is moving to the left, ℎ * is assigned with ℎ !!! or its value in the right side.This condition guarantees non-dry condition for the total depth, see [12,13].For simulating effects of bottom motion, instead of discretizing directly the eq.( 1) in the form of ℎ, we use the surface elevation !dan depth !, where ℎ = !+ !, such that the discretization of eq. ( 1) can be rewritted as follows ( ) The discretization in the momentum equation ( 2) is given as follows The last term in the left hand side of ( 7) is the nonlinear friction term that is approximated as follows Here, ℎ in eq. ( 8) is calculated as an averaged value in the half grid, that is defined as In eq. ( 7), a special treatment is given for calculating the advection term !! !!.As suggested in [12,13], the term is crucial and must treated carefully.Here, the term is calculated by using horizonal momentum != ℎ!, such that !! !!can be described as ( ) Since ℎ is always positive, is guaranteed by the upwind condition in (5), the expression (10) at time !!can be approximated by the following expression where ( ) ( ) Note that the horizontal velocity ! is discretized in the half grid, i.e. !!!!/! .Similar with the approximation for ℎ * , the value for !* is in the full grid, that is approximated using upwind condition as follows , jika 0 , jika 0 In summary, for an initial value problem as well as for bottom motion problem, the SWE in ( 1) and ( 2) can be solved as follows.
2. By using output from step (1), i.e. !!!!! and ℎ !!!! , and initial horizontal velocity !!, we calculate the momentum eq. ( 2) by using (7).Here, the advection term is calculated by using (11) and expressions (12,13,14).From this second step, we obtain For testing the capability of the numerical implementation for simulating physical tsunami phenomena such as wave runup and bottom motion, we test the numerical implementation with 3 test cases for wave runup and a test case of underwater landslide in the next section.

IV. RESULTS OF NUMERICAL IMPLEMENTATION
In a tsunami event generated by underwater landslide, not only the wave generation process that is important, but also the capability of the numerical model for simulating wave runup phenomenon is also crucial.In this section, to test the capability of the numerical implementation, we simulate 3 runup cases.The first one is the classical harmonic wave runup of Carrier and Greenspan 1958 [2], the others are breaking and non-breaking solitary wave runup performed by Synolakis [14].Finally, to test the performance of wave generated by bottom motion, we simulate an underwater landslide proposed by Lynett & Liu [9].

A. Harmonic Wave Runup
The first test case for wave runup is the classical problem introduced by Carrier and Greenspan 1958 [2].A harmonic wave with amplitude of 0.003m and period of 10s is propagating into a linearly sloping topography with slope of 1:25.An analytical solution is available, provided by [2].The initial depth is != 0.5m.For the numerical setting, we use a spatial discretization Δ! = 0.04m and with a temporal discretization Δ! = 0.05s.The simulation is run for 80s.The harmonic wave is influxed into the domain of computation from the left boundary at ! = −30.5musing an embedded influxing that is proposed by Liam e.a.2014 [8].In Fig. 3, we compare results of the simulation at many snapshot time with envelopes of the analytical solution of Carrier and Greenspan [2].It can be seen from the figure that the analytical solution agree quite well with the results of the simulation.Moreover, in Fig. 4, we compare horizontal shore line movement between the analytical solution with the result of the simulation.The figure shows a good agreement between the results of simulation with the analytical solution.

B. Solitary Wave Runup
The second test case is the runup of solitary waves, for non-breaking as well as breaking case.The test cases are based on the study of Synolakis 1986 [14], where experimental data from hydrodynamic laboratory are available for comparison.In his study, many physical experiment of solitary waves are performed.In this paper, we only use two cases, i.e. the non-breaking case; where the ratio between waveheight ! and the depth ! is !/! =0.0185 and the breaking case where !/! =0.To show quantitative comparison between results of simulation with the experiment data of Synolakis [14], we calculate the Root Mean Square Error (RMSE) and the Correlation coefficient (CorrCoef) for the nonbreaking and breaking wave runup.These values are summarized in Table I and II for the non-breaking case and the breaking case, respectively.From the Table I, the RMSE values show a very good agreement between results of simulation of non-breaking solitary runup with experiment data of Synolakis, while for the Correlation Coefficient values, the values also show a good agreement, i.e.CorrCoef > 0.9 for all snapshot times.
For the breaking solitary wave runup, the RMSE values show higher errors, especially during the runup process at ! = 3.19, 4.78, 6.38s.For the last two time frames, i.e. != 7.98s and 9.57!, the RMSE values are better.The correlation coefficient values also show similar trend.Here, the initial waveheight is much higher than the non-breaking case, results in a steeper wave form when the solitary wave reaches the sloping bottom.As the dispersive property is missing in the NSWE, i.e. all wave frequencies is moving with the same velocity and the nonlinear effects are not compensated by dispersive effects, results in steep wave crest especially at ! = 3.19 and 4.78!.At ! = 6.38! the result of simulation shows that the wave is already breaking, while the experimental data shows higher waveheight and nearly breaking.After the breaking process, at ! = 7.98!, the simulation can follow again the experiment, just as well at ! = 9.57!.

C. Underwater Landslide
To show the capability of the numerical scheme for simulation bottom motion phenomenon, we test the numerical implementation for simulating the wave generated by underwater landslide that is proposed by Lynett & Liu. 2002 [9] and later by Fuhrman & Madsen [4].Both in Lynett & Liu [9] and Fuhrman & Madsen [4], nonlinear and dispersive models are used for simulating the underwater landslide.Both models are compared with results of simulation by using a Boundary Integral Equation Model (BIEM) which solve potential flow.The BIEM is introduced in Grilli e.a.1989 [5].In this paper, we use simpler model, i.e. the non-dispersive NSWE for simulating the underwater landslide phenomenon.Let ! and !denote the horizontal coordinate and time, respectively.As described in [9], the dynamic of the bottom motion is described by ℎ(!, !) as follows The initial position a sliding mass is denoted by !!, where !(!) corresponds to motion of the sliding mass parallel with the slope that is defined as An initial acceleration for the landslide and terminal landslide velocity are denoted by !! and !!, respectively.They are defined as For the numerical setting, zero initial conditions for surface elevation and velocity are used.We use spatial discretization Δ! = 0.05! and temporal discretization Δ! = 0.005!.The simulation is performed for 7s.Snapshot of the surface elevation generated by underwater landslide is shown in the upper plot of Fig. 9, whereas the corresponding position of the bottom is shown in the lower plot.As in [9,4], we also compare the NSWE simulation with results of simulation of BIEM at several snapshot times, i.e. != 1.51, 3, 4.51, and 5.86!.These comparisons are shown in Fig. 10.From the figure, it shows that results of NSWE simulation match quite well with the results of BIEM as in Lynett & Liu [9].We also calculate quantitative comparison for this case, i.e. the RMSE and CorrCoef values as in the runup cases.Results of comparison are shown in Table III.In the table, even though the RMSE values tend to increase over the time, but it is considered to be relatively very low, as also can seen qualitatively in Fig. 10.Similarly with the CorrCoef values, the CorrCoef values tend to decrease over the time, but it shows relatively very high correlation values.As visually noticed, the RMSE values are very low and Correlation Coefficients are very high (CorrCoef > 0.98), which show that the results of the NSWE simulation with staggered grid agree very well with the results of BIEM in [9].

V. CONCLUSION
In this paper we have modified the staggered grid scheme proposed by Stelling & Duimeijer [12] for simulating bottom motion phenomenon, especially for case of wave generated by underwater landslide.Three runup cases have been performed to show that the numerical implementation is able to simulate non-breaking as well as breaking wave runup.In the test case of breaking solitary wave of Synolakis [14], with !/! = 0.3, the wave breaking process is captured automatically without introducing any ad hoc treatment in the numerical implementation.In Lynett & Liu [9] and in Fuhrman & Madsen [4], the wave simulations of underwater landslide are performed by using nonlinear and dispersive wave models, i.e.Boussinesq type of model, in which much more complicated than the Nonlinear Shallow Water Equations (NSWE).By using the numerical implementation proposed in this paper, the rather simple NSWE can reconstruct the wave generation by underwater landslide as in [9,4] with very small error, i.e.RMSE of 0.0002 and high correlation value, i.e. 0.99.

Fig. 1 .
Fig. 1.Illustration of variables for the shallow water equations.

Fig. 2 .
Fig. 2. Illustration of staggered grid scheme.The time discretization is denoted by integer n, hence !!!≈ !(! !, ! ! ) is an approximation value for ! at full grid !! and at time !! .In the staggered grid scheme, the surface elevation is placed in the full-grid, while the
and !!!!/! !!! , values in the next time step are obtained by repeating step 1 and 2 above.

Fig. 3 .
Fig. 3. Comparison between envelopes of the analytical solution (solid black lines) with results of simulation at many time frames (solid grey lines).

Fig. 4 .
Fig. 4. Comparison between horizontal shoreline movements of the analytical solution (solid black line) with the result of simulation (dashed red line).

3 .
Here, we consider a domain of computation !∈ [−100!, 10!].The center of the solitary wave is placed at ! != −40!.In the right side of the domain, a sloping bottom with slope 1:19.85 is placed, starting a !!= −19.85!.The illustration of the domain of computation is shown in Fig. 5.For the non-breaking case, i.e. !/! =0.0185, snapshots of the results of simulation at time != 12.77, 15.96, 19.15, 22.34!, are shown in Fig. 6.As shown in these plots, results of the simulation can follow the laboratory experiment data relatively very well during the runup and the rundown process.The runup process for the breaking case !/! =0.3 is shown in Fig. 7.

Fig. 5 .
Fig. 5. Illustration of domain of simulation for solitary wave runup of Synolakis.

Fig. 6 .
Fig. 6.Snapshot of non-breaking solitary wave runup of Synolakis for case H/d = 0.0185 at various time.The experiment data are denoted by blue dots, while the simulation is denoted by solid red line.

Fig. 7 .
Fig. 7. Snapshot of breaking solitary wave runup of Synolakis for case H/d = 0.3 at various time.The experiment data are denoted by blue dots, while the simulation is denoted by solid red line.
where != tan ! is the slope of the beach where the underwater landslide is occured, ℎ′ is defined as proposed in Fuhrman & Madsen 2009 [4], we use the following coefficients γ = 2, C d = C m = 1, b = 1 m, dan Δh = 0,05 m, dan != 6 ! .An illustration of the underwater landslide and its parameters is shown in Fig. 8.

Fig. 9 .
Fig. 9. Snapshot of wave simulation at t=5.9s for the underwater landslide (upper plot) and the bottom profile (lower plot).

TABLE I RMSE
AND CORRELATION COEFFICIENT FOR NON-BREAKING CASE OF SYNOLAKIS

TABLE III RMSE
AND CORRELATION COEFFICIENT FOR UNDERWATER LANDSLIDE.NON-BREAKING CASE OF SYNOLAKIS Ind. Journal on Computing Vol. 4, Issue. 1, March 2019