Nonlinear deformation and run-up of single tsunami waves of positive polarity: numerical simulations and analytical predictions

The estimate of individual wave run-up is especially important for tsunami warning and risk assessment as it allows to evaluate the inundation area. Here as a model of tsunami we use the long single wave of positive polarity. The period of such wave is rather long which makes it different from the famous Korteweg-de Vries soliton. This wave nonlinearly deforms during its propagation in the ocean, what results in a steep wave front formation. Situations, when waves approach the coast with a steep front are often observed during large tsunamis, e.g. 2004 Indian Ocean and 2011 15 Tohoku tsunamis. Here we study the nonlinear deformation and run-up of long single waves of positive polarity in the conjoined water basin, which consists of the constant depth section and a plane beach. The work is performed numerically and analytically in the framework of the nonlinear shallow water theory. Analytically, wave propagation along the constant depth section and its run-up on a beach are considered independently without taking into account wave interaction with the toe of the bottom slope. The propagation along the bottom of constant depth is described by Riemann wave, while the wave 20 run-up on a plane beach is calculated using rigorous analytical solutions of the nonlinear shallow water theory following the Carrier-Greenspan approach. Numerically, we use the finite volume method with the second order UNO2 reconstruction in space and the third order Runge-Kutta scheme with locally adaptive time steps. During wave propagation along the constant depth section, the wave becomes asymmetric with a steep wave front. Shown, that the maximum run-up height depends on the front steepness of the incoming wave approaching the toe of the bottom slope. The corresponding formula for maximum 25 run-up height, which takes into account the wave front steepness, is proposed.


Introduction
Evaluation of wave run-up characteristics is one of the most important tasks in coastal oceanography especially when estimating tsunami hazard. This knowledge is required as for planning coastal structures and protection works, as for shortterm tsunami forecast and tsunami warning. Its importance is also confirmed by a number of scientific papers, e.g. see recent 30 works (Tang et al. 2017;Touhami and Khellaf 2017;Zainali et al. 2017;Raz et al. 2018;Yao et al. 2018).
Many of these analytical formulas have been validated experimentally in laboratory tanks (Synolakis 1987, Li andRaichlen 2002;Lin et al. 2009;Didenkulova et al. 2013). For most of them, the solitary waves have been used. The soliton is rather easy to generate in the flume, therefore, laboratory studies of run-up of solitons are the most popular. However, (Madsen et 40 al. 2008) pointed out that the solitons are inappropriate to describe the real tsunami and proposed to use waves of longer duration than solitons and downscaled records of real tsunami. Schimmels et al. (2016) and Sriram et al. (2016)  It should be mentioned that the shape of tsunami varies a lot depending on its origin and the propagation path. One of the 45 best examples of tsunami wave shape variability is given in Shuto (1985) for the 1983 Japan Sea tsunami, where the same tsunami event resulted in very different tsunami approaches in different locations along Japanese coast. These wave shapes included: single positive pulses, undergoing both surging and spilling breaking scenarios, breaking bores, periodic wave trains, surging and breaking as well, a sequence of two or three waves and undular bores. This is why there is no such term as "typical tsunami wave shape", and therefore in the papers on wave run-up cited above many different wave shapes, such 50 as single pulses, N-waves, periodic symmetric and asymmetric wave trains, are considered. In this paper, we focus on the nonlinear deformation and run-up of long single pulses of positive polarity on a plane beach.
A similar study was performed for periodic sine waves (Didenkulova et al. 2007;Didenkulova 2009). It was shown that the run-up height increases with an increase in the wave asymmetry (wave front steepness) which is a result of nonlinear wave deformation during its propagation in a basin of constant depth. It was found analytically that the run-up height of this 55 nonlinearly deformed sine wave is proportional to the square root of the wave front steepness. Later on, this result was also confirmed experimentally (Didenkulova et al. 2013).
It should be noted that these analytical finding also match tsunami observations. Steep tsunami waves are often witnessed and reported during large tsunami events, such as 2004 Indian Ocean and 2011 Tohoku tsunamis. Sometimes the wave, which approaches the coast, represents a "wall of water" or a bore, which is demonstrated by numerous photos and videos of 60 these events.
The nonlinear steepening of the long single waves of positive polarity has also been observed experimentally in ), but its effect on wave run-up has not been studied yet. In this paper, we study this effect both analytically and numerically. Analytically, we apply the methodology developed in (Didenkulova 2009;Didenkulova et al. 2014), where we consider the processes of wave propagation in the basin of constant depth and the following wave run-up on a plane beach 3 independently, not taking into account the point of merging of these two bathymetries. Numerically, we solve the nonlinear shallow water equations.
The paper is organized as follows. In Section 2, we give the main formulas and briefly describe the analytical solution. The numerical model is described and validated in Section 3. The nonlinear deformation and run-up of the long single wave of positive polarity is described in Section 4. The main results are summarized in Section 5. 70

Analytical solution
We solve the nonlinear shallow water equations for the bathymetry shown in Fig. 1: (2) 75 Here (x, t) is the vertical displacement of the water surface with respect to the still water level, u(x,t)depth-averaged water flow, h(x)unperturbed water depth, g is the gravitational acceleration, x is the coordinate directed onshore, and t is time. The system of Eqs.(1),(2) is solved independently for the two bathymetries shown in Fig. 1: a basin of constant depth h 0 and length X 0 and a plane beach, where the water depth h(x) = -x tan.
Eqs. (1),(2) can be solved exactly for a few specific cases. In the case of constant depth, the solution is described by the 80 Riemann wave (Stoker 1957). Its adaptation for the boundary problem can be found in Zahibo et al. (2008). In the case of a plane beach, the corresponding solution was found by Carrier and Greenspan (1958). Both solutions are well-known and widely used and we do not reproduce them here, but just provide some key formulas.
As already mentioned, during its propagation along the basin of constant depth h 0 , the wave transforms as a Riemann wave (Zahibo et al. 2008): 85 where 00 ( is the water displacement at the left boundary. After the propagation over the section of constant depth h 0 , the incident wave has the following shape: Following the methodology developed in , we let this nonlinearly deformed wave described by Eq. (5) run-up on a plane beach, characterized by the water depth h(x) =x tan. This approach does not take into account the merging point of the two bathymetries and, therefore, does not account for reflection from the toe of the slope and wave interaction with the reflected wave. To do this we represent the input wave 0 X  as a Fourier integral: Its complex spectrum B(ω) can be found in an explicit form in terms of the inverse Fourier transform: Eq. (7) can be re-written in terms of the water displacement, produced by the wave maker at the left boundary (Zahibo et al.

2008): 105
In this study we consider long single pulses of positive polarity: 5 where A is the input wave height and T is the effective wave period at the location with the water depth h 0 . The wave described by Eq. (9) has an arbitrary height and period and, therefore, does not satisfy properties of the soliton, but just has a 110 sech 2 shape. Substituting Eq. (9) into Eq. (8), we can calculate the complex spectrum B(ω).
Wave run-up oscillations at the coast r(t) and the velocity of the moving shoreline u(t) can be found from ): This solution we also compare with the run-up of a single wave of positive polarity described by Eq. (9) (without nonlinear deformation). The maximum run-up height R max of such wave (9) can be found from Sriram et al. 120 2016): If the initial wave is soliton, Eq. (14) coincides with the famous Synolakis formula (Synolakis, 1987).

Numerical model
Numerically, we solve the nonlinear shallow water equations Eqs. (1),(2) written in a conservative form for a total water 125 depth. We include the effect of the varying bathymetry (in space) and neglect all friction effects. However, the resulting numerical model will take into account for some dissipation thanks to the numerical scheme dissipation, which is necessary for the stability of the scheme and should not influence much run-up characteristics. Namely, we employ the natural numerical method, which was developed especially for conservation laws -the finite volume schemes.
The numerical scheme is based on the second order in space UNO2 reconstruction, which is briefly described in (Dutykh et 130 al. 2011b). In time we employ the third order Runge-Kutta scheme with locally adaptive time steps in order to satisfy the CFL stability condition. The numerical technique to simulate the wave run-up was described previously in (Dutykh et al. 6 2011a). The bathymetry source term is discretized using the hydrostatic reconstruction technique, which implies the wellbalanced property of the numerical scheme (Gosse, 2013).  Fig. 2) as also on the beach (Fig. 3). The comparison of run-up height calculated numerically and analytically using the approach described in Section 2 145 with the experimental record in shown in Fig. 3. It can be seen that the experimentally recorded wave is slightly smaller which may be caused by the bottom friction and especially on the slope. Both numerical and analytical models describe the first wave of positive polarity rather well. The numerical prediction of run-up height is slightly higher than the analytical one. This additional increase in the run-up height in numerical model may be explained by the nonlinear interaction with the reflected wave, which is not taken into account in the analytical model. The wave of negative polarity is much more sensitive 150 to all the effects mentioned above than the wave of positive polarity and, therefore, looks different for all three lines in Fig. 3. By introducing additional dissipation in numerical model one can easily reach perfect agreement between the numerical simulations and experimental data. However, we do not do so, since below we are focusing on the analysis of analytical results and for clarity would like to avoid additional parameters in the numerical model. Also, we focus on the maximum run-up height and, therefore, expect small differences between the results of analytical and numerical models. 155

Results of numerical and analytical calculations 160
It is reported in (Didenkulova et al 2007;Didenkulova 2009) for a periodic sine wave, that the extreme run-up height increases proportionally to the square root of the wave front steepness. In this section, we study the nonlinear deformation and steepening of waves described by Eq. (9) and its effect on the extreme wave run-up height. The corresponding bathymetry used in analytical and numerical calculations is normalized on the water depth in the section of constant depth h 0 , and is shown in Fig. 1. The input wave parameters such as wave amplitude, A/h 0 , and effective wave length, λ/X 0 , where 165 0 T gh   , are changed. The beach slope is taken tan α = 1:20 for all simulations.
We underline that in order to have analytical solution, the criterion of no wave breaking should be satisfied. Therefore, all analytical and numerical calculations below are chosen for non-breaking waves. Fig. 4 shows the dimensionless maximum run-up height, R max /A as a function of the initial wave amplitude, A/h 0 . The incident wave propagates over different distances to the bottom slope, X 0 /λ = 1.7, 3.4, 5.1, 6.8; kh 0 = 0.38. Analytical 170 solution described in Section 2 is shown with lines, and numerical solution described in Section 3 is shown with symbols (diamonds, triangles, squares and circles). It can be seen that in most cases and especially for small values of X 0 /λ = 1.7 and 3.4, numerical simulations give larger run-up heights than analytical predictions. These differences can be explained by the effects of wave interaction with the toe of the underwater beach slope, which are not taken into account in the analytical solution. For larger distances X 0 /λ = 6.8, both analytical and numerical solutions give similar results, supported by the 175 numerical scheme dissipation described in Section 3, which can be considered a "numerical error". It should be mentioned that we use zero physical dissipation rate for these simulations, however, small dissipation for stability of the numerical 8 scheme is still needed and this may become noticeable at large distances. For the sech 2 -shaped wave (A/h 0 = 0.03, λ/X 0 = 0.12) propagation, the reduction of initial wave amplitude constitutes ~2 %.
It is worth mentioning, that for small initial wave amplitudes all run-up heights are close to each other and are close to the 180 thick black line, which corresponds to Eq. (14) for wave run-up on a beach without constant depth section. This means that the effects we are talking about are important only for nonlinear waves and irrelevant for weakly nonlinear or almost linear waves.
The same effects can be seen in Fig. 5, which shows the maximum run-up height, R max /A as a function of distance to the slope, X 0 /λ, for different amplitudes of the initial wave, A/h 0 . The distance X 0 /λ changes from 0.8 to 9.4; kh 0 = 0.38. The 185 analytical solution is shown with lines while the numerical solution is shown with symbols (triangles, squares and circles). It can be seen in Fig. 5, for smaller values of X 0 /λ < 6 numerical predictions provide relatively larger run-up values, as compared with analytical predictions, while for higher values of X 0 /λ > 6 the differences are significantly reduced. A relevant change of this behaviour is given for A/h 0 = 0.06. We can observe that numerical predictions for this amplitude become smaller than analytical predictions for X 0 /λ > 8. As stated above, we believe that this can be a result of interplay of 190 two effects: interaction with the underwater bottom slope, which is not taken into account in the analytical prediction and the numerical scheme dissipation ("numerical error"), which affects the numerical results.  The dependence of maximum run-up height, R max /A on kh 0 is shown in Fig. 6 for A/h 0 = 0.03. It can be seen that the difference between numerical and analytical results decreases with an increase in kh 0 . We relate this effect with the wave 205 interaction with the slope, which is not properly accounted in our analytical approach. As one can see in Fig. 7, this difference for a milder beach slope tan α = 1:50 is reduced. The next Fig. 8 supports all the conclusions drawn above. It also shows that difference between analytical and numerical results increases with an increase in wave period. As pointed out before for small wave periods the numerical solution may 215 coincide with the analytical one or even become smaller as it happens for kh 0 = 0.38 for X 0 /λ > 8.
In order to calculate the incident wave front steepness in the beginning of the bottom slope from results of numerical 230 simulations we should separate the incident wave and the wave reflected from the bottom slope. At the same time, the wave steepening along the basin of constant depth is very well described analytically as demonstrated in Fig. 9.

11
It can be seen that the wave transformation described by the analytical model is in a good agreement with numerical simulations. Therefore, below we reference to the analytically defined wave front steepness having in mind that it well coincides with the numerical one. Having said this, we approach the main result of this paper, which is shown in Fig. 10. The 235 red solid line gives the analytical prediction. It is universal for single waves of positive polarity for different amplitudes A/h 0 and kh 0 and can be well approximated by the power fit (coefficient of determination R-squared = 0.99): where R max /A is the maximum run-up height in the conjoined basin (with a section of constant depth), R 0 /A is the corresponding maximum run-up height on a plane beach (without a section of constant depth). 240 Figure 8: Maximum run-up height, R max /A as a function of the distance to the slope, X 0 /λ for different kh 0 . Analytical solution described in Section 2 is shown by lines and numerical solution described in Section 3 is shown by symbols (triangles, squares and circles) with matching colours; A/h 0 = 0.03.

245
The fit is shown in Fig. 10 by the black dashed line. For comparison, the dependence of the maximum run-up height on the wave front steepness obtained using the same method for a sine wave is stronger than for a single wave of positive polarity (Didenkulova et al. 2007) and is proportional to the square root of the wave front steepness. This is logical as sinusoidal wave has a sign-variable form and, therefore, excites a higher run-up. For possible mechanisms, see discussion on N-waves in (Tadepalli and Synolakis 1994). 250 The results of numerical simulations are shown in Fig. 10 with different markers. It can be seen that numerical data for the 265 same period, but different amplitudes follow the same curve. The run-up is higher for waves with smaller kh 0 . In our opinion, 13 this dependence on kh 0 is a result of merging plane beach with a flat bottom. This effect can be parameterized with the factor (L/λ) 1/4 . The result of this parameterization is shown in Fig. 11. Here we can see that for smaller face front wave steepness, s/s 0 < 1.5, the run-up height is proportional to the analytically estimated curve shown by Eq. (16), while for larger face front wave steepness, s/s 0 > 1.5, the dependence on s/s 0 is weaker. This dependence for all numerical run-up height data, presented 270 in Fig. 11, can be approximated by the power fit (coefficient of determination R-squared = 0.85):

Conclusions and Discussion
In this paper, we study the nonlinear deformation and run-up of tsunami waves, represented by single waves of positive polarity. We consider the conjoined water basin, which consists of a section of constant depth and a plane beach. While propagating in such basin, the wave shape changes forming a steep front. Tsunamis often approach the coast with a steep 280 wave front, as it was observed during large tsunami events, e.g. 2004 Indian Ocean Tsunami and 2011 Tohoku tsunami.
The study is performed both analytically and numerically in the framework of the nonlinear shallow water theory. The analytical solution considers nonlinear wave steepening in the constant depth section and wave run-up on a plane beach independently and, therefore, does not take into account wave interaction with the toe of the bottom slope. The propagation along the bottom of constant depth is described by Riemann wave, while the wave run-up on a plane beach is calculated 285 using rigorous analytical solutions of the nonlinear shallow water theory following the Carrier-Greenspan approach. The numerical scheme does not have this limitation. It employs the finite volume method and is based on the second order UNO2 14 reconstruction in space and the third order Runge-Kutta scheme with locally adaptive time steps. The model is validated against experimental data.
The main conclusions of the paper are the following. 290  Found analytically, that maximum tsunami run-up height on a beach depends on the wave front steepness at the toe of the bottom slope. This dependence is general for single waves of different amplitudes and periods and can be approximated by the power fit:  This dependence is slightly weaker than the corresponding dependence for a sine wave, proportional to the square root of the wave front steepness (Didenkulova et al. 2007). The stronger dependence of a sine wave run-up on the 295 wave front steepness is consistent with the philosophy of N-waves (Tadepalli and Synolakis 1994).
 Numerical simulations in general support this analytical finding. For smaller face front wave steepness (s/s 0 < 1.5) numerical curves of maximum tsunami run-up height are parallel to the analytical one, while for larger face front wave steepness (s/s 0 > 1.5), this dependence is milder. The latter may be a result of numerical dissipation (error), which is larger for a longer wave propagation and, consequently, larger wave steepness. The suggested formula, 300 which gives the best fit with the data of numerical simulations in general is      These results can also be used in tsunami forecast. Sometimes, in order to save time for tsunami forecast, especially for long distance wave propagation, the tsunami run-up height is not simulated directly, but estimated using analytical or empirical formulas (Glimsdal et al. 2019;Løvholt et al. 2012). In these cases we recommend using formulas, which take into account the face front wave steepness. The face front steepness of the approaching 305 tsunami wave can be estimated from the data of the virtual (computed) or real tide-gauge stations and then be used to estimate tsunami maximum run-up height on a beach.
The nonlinear shallow water equations which are used in this study and commonly utilized for tsunami modelling, are also known as to neglect dispersive effects. In this context, it is important to mention the recent work of Larsen and Fuhrman (2019). They used RANS equations and k-ω model for turbulence closure to simulate propagation and run-up of positive 310 single waves, including full resolution of dispersive short waves (and their breaking) that can develop near a positive tsunami front. They similarly showed that this effect depends on the propagation distance prior to the slope, if a simple toe with a slope type of bathymetry is utilized. This work shows that these short waves have little effect on the overall run-up, and hence give additional credence to the use of shallow water equations. These results largely confirm what was previously hypothesized by , that these short waves would have little effect on the overall run-up and inundation of 315 tsunamis (though they found that they could significantly increase the maximum flow velocities).