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

The estimate of an individual wave run-up is especially important for tsunami warning and risk assessment, as it allows for evaluating the inundation area. Here, as a model of tsunamis, we use the long single wave of positive polarity. The period of such a 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, which results in a steep wave front formation. Situations in which waves approach the coast with a steep front are often observed during large tsunamis, e.g. the 2004 Indian Ocean and 2011 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 the Riemann wave, while the wave 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. It is 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 run-up height, which takes into account the wave front steepness, is proposed.

1 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 for planning coastal structures and protection works as well as for short-term tsunami forecasts and tsunami warning. Its importance is also confirmed by a number of scientific papers (see recent works, e.g. Tang et al., 2017; Touhami and Khellaf, 2017; Zainali et al., 2017; Raz et al., 2018; Yao et al., 2018).

The general solution of the nonlinear shallow-water equations on a plane beach was found by Carrier and Greenspan (1958) using the hodograph transformation. Later on, many other authors found specific solutions for different types of waves climbing the beach (see, for instance, Pedersen and Gjevik, 1983; Synolakis, 1987; Synolakis et al., 1988; Mazova et al., 1991; Pelinovsky and Mazova, 1992; Tadepalli and Synolakis, 1994; Brocchini and Gentile, 2001; Carrier et al., 2003; Kânoğlu, 2004; Tinti and Tonini, 2005; Kânoğlu and Synolakis 2006; Madsen and Fuhrman, 2008; Didenkulova et al., 2007; Didenkulova, 2009; Madsen and Schäffer, 2010).

Many of these analytical formulas have been validated experimentally in laboratory tanks (Synolakis, 1987; Li and Raichlen, 2002; Lin et al., 1999; 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 al., 2008) pointed out that the solitons are inappropriate for describing 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) generated such long waves in the Large Wave Flume of Hanover (GWK FZK) using the piston type of wave maker, while McGovern et al. (2018) did it using the pneumatic wave generator.

It should be mentioned that the shape of tsunami varies a lot depending on its origin and the propagation path. One of the best examples of tsunami wave shape variability is given in Shuto (1985) for the 1983 Sea of Japan tsunami, where the same tsunami event resulted in very different tsunami approaches in different locations along the Japanese coast. These wave shapes included the following: single positive pulses, undergoing both surging and spilling breaking scenarios; breaking bores; periodic wave trains, surging as well as breaking; and a sequence of two or three waves and undular bores. This is why there is no “typical tsunami wave shape”, and therefore in the papers on wave run-up cited above, many different wave shapes, such as single pulses, N waves, and 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 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 findings 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 these events.

The nonlinear steepening of the long single waves of positive polarity has also been observed experimentally in Sriram et al. (2016), 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) and 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 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 Sect. 2, we give the main formulas and briefly describe the analytical solution. The numerical model is described and validated in Sect. 3. The nonlinear deformation and run-up of the long single wave of positive polarity are described in Sect. 4. The main results are summarized in Sect. 5.

2 Analytical solution

We solve the nonlinear shallow-water equations for the bathymetry shown in Fig. 1:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial u}{\partial t}}+u{\displaystyle \frac{\partial u}{\partial x}}+g{\displaystyle \frac{\partial \mathit{\eta }}{\partial x}}=\mathrm{0},\text{(2)}& {\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left[\left(h\left(x\right)+\mathit{\eta }\right)u\right]=\mathrm{0}.\end{array}$$

Here η(x,t) is the vertical displacement of the water surface with respect to the still water level, u(x,t) is the depth-averaged water flow, h(x) is the unperturbed water depth, g is the gravitational acceleration, x is the coordinate directed onshore and t is time. The system of Eqs. (1) and (2) is solved independently for the two bathymetries shown in Fig. 1: a basin of constant depth h₀ and length X₀ and a plane beach, where the water depth $h\left(x\right)=-x$ tan α.

Equations (1) and (2) can be solved exactly for a few specific cases. In the case of constant depth, the solution is described by the 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₀, the wave transforms as a Riemann wave (Zahibo et al., 2008):

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathit{\eta }\left(x,t\right)={\mathit{\eta }}_{\mathrm{0}}\left[t-{\displaystyle \frac{x+X_{\mathrm{0}}+L}{V\left(x,t\right)}}\right],\text{(4)}& {\displaystyle }V\left(x,t\right)=\mathrm{3}\sqrt{g\left[h_{\mathrm{0}}+\mathit{\eta }\left(x,t\right)\right]}-\mathrm{2}\sqrt{g{h}_{\mathrm{0}}},\end{array}$$

where ${\mathit{\eta }}_{\mathrm{0}}(x=-L-X_{\mathrm{0}},t)$ is the water displacement at the left boundary. After the propagation over the section of constant depth h₀, the incident wave has the following shape:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathit{\eta }}_{X\mathrm{0}}\left(t\right)& ={\mathit{\eta }}_{\mathrm{0}}\left[t-{\displaystyle \frac{X_{\mathrm{0}}}{V\left(x,t\right)}}\right],V_{X\mathrm{0}}\left(t\right)\\ & =\mathrm{3}\sqrt{g\left[h_{\mathrm{0}}+{\mathit{\eta }}_{X\mathrm{0}}\left(t\right)\right]}-\mathrm{2}\sqrt{g{h}_{\mathrm{0}}}.\end{array}\end{array}$$

Following the methodology developed in Didenkulova (2008), 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.

Figure 1 Bathymetry sketch. The wavy curve at the toe of the slope regards analytical solution, which does not take into account merging between the constant depth and sloping beach sections.

To do this, we represent the input wave η_X0 as a Fourier integral:

$$\begin{array}{}\text{(6)}& {\mathit{\eta }}_{X\mathrm{0}}=\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}B\left(\mathit{\omega }\right)\exp \left(i\mathit{\omega }t\right)\mathrm{d}\mathit{\omega }.\end{array}$$

Its complex spectrum B(ω) can be found in an explicit form in terms of the inverse Fourier transform:

$$\begin{array}{}\text{(7)}& B\left(\mathit{\omega }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{\mathit{\eta }}_{X\mathrm{0}}\left(t\right)\exp \left(-i\mathit{\omega }t\right)\mathrm{d}t.\end{array}$$

Equation (7) can be rewritten in terms of the water displacement, produced by the wave maker at the left boundary (Zahibo et al., 2008):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}B\left(\mathit{\omega }\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }i\mathit{\omega }}}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_{\mathrm{0}}}{\mathrm{d}z}}\exp \left(-i\mathit{\omega }\left[z+{\displaystyle \frac{x+X_{\mathrm{0}}+L}{V\left({\mathit{\eta }}_{\mathrm{0}}\right)}}\right]\right),\\ & \mathrm{d}z,z=t-{\displaystyle \frac{x+X_{\mathrm{0}}+L}{V\left({\mathit{\eta }}_{\mathrm{0}}\right)}}.\end{array}\end{array}$$

In this study we consider long single pulses of positive polarity:

$$\begin{array}{}\text{(9)}& {\mathit{\eta }}_{\mathrm{0}}\left(t\right)=A{\mathrm{sech}}^{\mathrm{2}}\left({\displaystyle \frac{t}{T}}\right),\end{array}$$

where A is the input wave height and T is the effective wave period at the location with the water depth h₀. 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 sech² 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 Didenkulova et al. (2008):

$$\begin{array}{}\text{(10)}& {\displaystyle }r\left(t\right)=R\left(t+{\displaystyle \frac{u}{g\tan \mathit{\alpha }}}\right)-{\displaystyle \frac{u^{\mathrm{2}}}{\mathrm{2}g}},\text{(11)}& {\displaystyle }u\left(t\right)=U\left(t+{\displaystyle \frac{u\left(t\right)}{g\tan \mathit{\alpha }}}\right),\text{(12)}& {\displaystyle }\begin{array}{rl}R\left(t\right)& =\sqrt{\mathrm{2}\mathit{\pi }\mathit{\tau }\left(L\right)}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\sqrt{\mathrm{|}\mathit{\omega }\mathrm{|}}H\left(\mathit{\omega }\right)\\ & \exp \left\{i\left(\mathit{\omega }\left(t-\mathit{\tau }\left(L\right)\right)+{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\mathrm{sign}\left(\mathit{\omega }\right)\right)\right\}\mathrm{d}\mathit{\omega },\end{array}\text{(13)}& {\displaystyle }U\left(t\right)={\displaystyle \frac{\mathrm{1}}{\tan \mathit{\alpha }}}{\displaystyle \frac{\mathrm{d}R}{\mathrm{d}t}},\end{array}$$

where $\mathit{\tau }=\mathrm{2}L/\sqrt{g{h}_{\mathrm{0}}}$ is the travel time to the coast.

We also compare this solution 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 a wave (Eq. 9) can be found from Didenkulova et al. (2008) and Sriram et al. (2016):

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{R_{max}}{A}}=\mathrm{2.8312}\sqrt{\cot \mathit{\alpha }}{\left({\displaystyle \frac{\mathrm{1}}{g{h}_{\mathrm{0}}}}{\left({\displaystyle \frac{\mathrm{2}{h}_{\mathrm{0}}}{\sqrt{\mathrm{3}T}}}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{4}}.\end{array}$$

If the initial wave is a soliton, Eq. (14) coincides with the famous Synolakis formula (Synolakis, 1987).

3 Numerical model

Numerically, we solve the nonlinear shallow-water equations Eqs. (1) and (2), written in a conservative form for a total water depth. We include the effect of the varying bathymetry (in space) and neglect all friction effects. However, the resulting numerical model will be taken into account for some dissipation thanks to the numerical scheme dissipation, which is necessary for the stability of the scheme and should not influence many 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 al. (2011b). In time we employ the third-order Runge–Kutta scheme with locally adaptive time steps in order to satisfy the Courant–Friedrichs–Lewy stability condition along with the local error estimator to bound the error term to the prescribed tolerance parameter. The numerical technique to simulate the wave run-up was described previously in Dutykh et al. (2011a). The bathymetry source term is discretized using the hydrostatic reconstruction technique, which implies the well-balanced property of the numerical scheme (Gosse, 2013).

Figure 2 Water elevations along the 251 m long constant depth section of the Large Wave Flume (GWK), where h₀=3.5 m, A=0.1 m, T=20 s and tan $\mathit{\alpha }=\mathrm{1}:\mathrm{6}$. Results of numerical simulations are shown by the red line, and experimental data are shown by the blue line.

The numerical scheme is validated against experimental data of wave propagation and run-up in the Large Wave Flume (GWK) in Hanover, Germany. The experiments were set with a flat bottom, with a constant depth of h₀=3.5 m, length of [a, b] =251 m and a plane beach with a slope of tan $\mathit{\alpha }=\mathrm{1}:\mathrm{6}$ (see Fig. 1). The flume had 16 wave gauges along the constant depth section and a run-up gauge on the slope. The incident wave had an amplitude of A=0.1 m and period of T=20 s. The detailed description of the experiments can be found in Didenkulova et al. (2013). The results of numerical simulations are in good agreement with the laboratory experiments along the constant depth section (see Fig. 2) as well as on the beach (Fig. 3). The comparison of the run-up height is calculated numerically and analytically using the approach described in Sect. 2 and with the experimental record shown in Fig. 3. It can be seen that the experimentally recorded wave is slightly smaller, which may be caused by the bottom friction, 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 the 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 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 the numerical model, one can easily reach perfect agreement between the numerical simulations and experimental data. However, we do not do so, since below we focus 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. The data used for all figures of this paper are available at https://doi.org/10.13140/rg.2.2.27658.41922 (Abdalazeez et al., 2019).

Figure 3 Run-up height of the long single wave with A=0.1 m and T=20 s on a beach slope, where tan $\mathit{\alpha }=\mathrm{1}:\mathrm{6}$. The numerical solution is shown by the red dotted line, the analytical solution is shown by the blue dashed line and the experimental record is shown by the black solid line.

4 Results of numerical and analytical calculations

It is reported in Didenkulova et al. (2007) and Didenkulova (2009), for a periodic sine wave, that the extreme run-up height increases proportionally with 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 their 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₀ and is shown in Fig. 1. The input wave parameters such as wave amplitude, A∕h₀, and effective wave length, λ∕X₀, where $\mathit{\lambda }=T\sqrt{g{h}_{\mathrm{0}}}$, are changed. The beach slope is taken as tan $\mathit{\alpha }=\mathrm{1}:\mathrm{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.

Figure 4 Maximum run-up height, R_max∕A, as a function of initial wave amplitude, A∕h₀, for different distances to the slope, X₀∕λ. Analytical solution described in Sect. 2 is shown by lines, and numerical solution described in Sect. 3 is shown by symbols (diamonds, triangles, squares and circles) with matching colours. The thick black line corresponds to Eq. (14) for wave run-up on a beach without constant depth section, where kh₀=0.38.

Figure 4 shows the dimensionless maximum run-up height, R_max∕A, as a function of the initial wave amplitude, A∕h₀. The incident wave propagates over different distances to the bottom slope, $X_{\mathrm{0}}/\mathit{\lambda }=\mathrm{1.7}$, 3.4, 5.1 and 6.8, where kh₀=0.38. The analytical solution described in Sect. 2 is shown with lines, and the numerical solution described in Sect. 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_{\mathrm{0}}/\mathit{\lambda }=\mathrm{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_{\mathrm{0}}/\mathit{\lambda }=\mathrm{6.8}$, both analytical and numerical solutions give similar results, supported by the numerical scheme dissipation described in Sect. 3, which can be considered a “numerical error”. It should be mentioned that we use a physical dissipation rate of zero for these simulations; however, a small dissipation for stability of the numerical scheme is still needed, and this may become noticeable at large distances. For the sech²-shaped wave ($A/h_{\mathrm{0}}=\mathrm{0.03}$, $\mathit{\lambda }/X_{\mathrm{0}}=\mathrm{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 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.

Figure 5 Maximum run-up height, R_max∕A, as a function of distance to the slope, X₀∕λ, for different amplitudes of the initial wave, A∕h₀. Analytical solution described in Sect. 2 is shown by lines, and numerical solution described in Sect. 3 is shown by symbols (triangles, squares and circles) with matching colours, where kh₀=0.38.

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₀∕λ, for different amplitudes of the initial wave, A∕h₀. The distance X₀∕λ changes from 0.8 to 9.4, where kh₀=0.38. The 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 that for smaller values of X₀∕λ < 6, numerical predictions provide relatively larger run-up values compared with analytical predictions, while for higher values of X₀∕λ > 6, the differences are significantly reduced. A relevant change of this behaviour is given for $A/h_{\mathrm{0}}=\mathrm{0.06}$. We can observe that numerical predictions for this amplitude become smaller than analytical predictions for X₀∕λ > 8. As stated above, we believe that this can be a result of interplay of 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₀ is shown in Fig. 6 for $A/h_{\mathrm{0}}=\mathrm{0.03}$. It can be seen that the difference between numerical and analytical results decreases with an increase in kh₀. We relate this effect with the wave 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 $\mathit{\alpha }=\mathrm{1}:\mathrm{50}$ is reduced.

Figure 6 Maximum run-up height, R_max∕A, as a function of kh₀ for different distances to the slope, X₀∕λ. Analytical solution described in Sect. 2 is shown by lines, and numerical solution described in Sect. 3 is shown by symbols (diamonds, triangles, squares and circles) with matching colours. The thick black line corresponds to Eq. (14) for wave run-up on a beach without constant depth section ($A/h_{\mathrm{0}}=\mathrm{0.03}$).

The next figure, Fig. 8, supports all the conclusions drawn above. It also shows that the difference between analytical and numerical results increases with an increase in the wave period. As pointed out before for small wave periods, the numerical solution may coincide with the analytical one or even become smaller as in kh₀=0.38 for X₀∕λ > 8.

Figure 7 Maximum run-up height, R_max∕A, as a function of initial effective wave length, λ∕X₀ (blue axes) and kh₀ (black axes). Analytical solutions for tan $\mathit{\alpha }=\mathrm{1}:\mathrm{20}$ and tan $\mathit{\alpha }=\mathrm{1}:\mathrm{50}$ are shown by dotted and dashed lines, respectively, while numerical simulations for tan $\mathit{\alpha }=\mathrm{1}:\mathrm{20}$ and tan $\mathit{\alpha }=\mathrm{1}:\mathrm{50}$ are shown by circles and crosses, respectively ($A/h_{\mathrm{0}}=\mathrm{0.03}$).

It is important that both analytical and numerical results in Figs. 5 and 8 demonstrate an increase in maximum run-up height with an increase in the distance X₀∕λ. This result is in agreement with the conclusions of Didenkulova et al. (2007) and Didenkulova (2009) for sinusoidal waves. In order to be consistent with the results of Didenkulova et al. (2007) and Didenkulova (2009), we connect the distance X₀∕λ with the incident wave front steepness in the beginning of the bottom slope. The wave front steepness s is defined as the maximum of the time derivative of water displacement, $\mathrm{d}(\mathit{\eta }/A)/\mathrm{d}(t/T)$, and is studied in relation with the initial wave front steepness, s₀, where

$$\begin{array}{}\text{(15)}& s\left(x\right)={\displaystyle \frac{max\left(\mathrm{d}\mathit{\eta }(x,t)/\mathrm{d}t\right)}{A/T}},s_{\mathrm{0}}={\displaystyle \frac{max\left(\mathrm{d}\mathit{\eta }(x=a,t)/\mathrm{d}t\right)}{A/T}}.\end{array}$$

In order to calculate the incident wave front steepness in the beginning of the bottom slope from results of numerical 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.

It can be seen that the wave transformation described by the analytical model is in a good agreement with numerical simulations. Therefore, below we make reference to the analytically defined wave front steepness, keeping in mind that it coincides well with the numerical one. Having said this, we approach the main result of this paper, which is shown in Fig. 10. The red solid line gives the analytical prediction. It is universal for single waves of positive polarity for different amplitudes A∕h₀ and kh₀ and can be approximated well by the power fit (coefficient of determination R² =0.99):

$$\begin{array}{}\text{(16)}& R_{max}/R_{\mathrm{0}}={\left(s/s_{\mathrm{0}}\right)}^{\mathrm{0.42}},\end{array}$$

where R_max∕A is the maximum run-up height in the conjoined basin (with a section of constant depth); R₀∕A is the corresponding maximum run-up height on a plane beach (without a section of constant depth).

Figure 8 Maximum run-up height, R_max∕A, as a function of the distance to the slope, X₀∕λ , for different values of kh₀. Analytical solution described in Sect. 2 is shown by lines, and numerical solution described in Sect. 3 is shown by symbols (triangles, squares and circles) with matching colours ($A/h_{\mathrm{0}}=\mathrm{0.03}$).

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 the sinusoidal wave has a sign-variable form and, therefore, excites a higher run-up. For possible mechanisms, see the discussion on N waves in Tadepalli and Synolakis (1994).

Figure 9 Wave evolution at different locations, $x/\mathit{\lambda }=\mathrm{0}$, 0.85, 1.71, 2.56, 3.41, 4.27 and 5.12, along the section of constant depth for a basin with $X_{\mathrm{0}}/\mathit{\lambda }=\mathrm{5.12}$ and tan $\mathit{\alpha }=\mathrm{1}:\mathrm{20}$. Numerical results are shown by solid lines, while the analytical predictions are given by the black dotted lines. The parameters of the wave are $A/h_{\mathrm{0}}=\mathrm{0.03}$ and kh₀=0.19.

Figure 10 The ratio of maximum run-up height in the conjoined basin, R_max∕A, and the maximum run-up height on a plane beach, R₀∕A, versus the wave front steepness, s∕s₀, for $A/h_{\mathrm{0}}=\mathrm{0.057}$, kh₀=0.38 (brown points); $A/h_{\mathrm{0}}=\mathrm{0.086}$, kh₀=0.38 (red plus signs); $A/h_{\mathrm{0}}=\mathrm{0.057}$, kh₀=0.29 (blue points); $A/h_{\mathrm{0}}=\mathrm{0.086}$, kh₀=0.29 (turquoise plus signs); $A/h_{\mathrm{0}}=\mathrm{0.057}$, kh₀=0.22 (violet points); $A/h_{\mathrm{0}}=\mathrm{0.086}$, kh₀=0.22 (pink plus signs); $A/h_{\mathrm{0}}=\mathrm{0.057}$, kh₀=0.19 (dark-green points); and $A/h_{\mathrm{0}}=\mathrm{0.086}$, kh₀=0.19 (light-green plus signs). All markers correspond to the results of numerical simulations, while the asymptotic analytical predictions are given by the red solid line. Black dashed line corresponds to the power fit of the analytical results of Eq. (16).

The results of numerical simulations are shown in Fig. 10 with different markers. It can be seen that numerical data for the same period but different amplitudes follow the same curve. The run-up is higher for waves with smaller kh₀. In our opinion, this dependence on kh₀ is a result of merging a plane beach with a flat bottom. This effect can be parameterized with the factor $(L/\mathit{\lambda }{)}^{\mathrm{1}/\mathrm{4}}$. The result of this parameterization is shown in Fig. 11. Here we can see that for smaller face front wave steepness, s∕s₀ < 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₀ > 1.5, the dependence on s∕s₀ is weaker. This dependence for all numerical run-up height data, presented in Fig. 11, can be approximated by the power fit (coefficient of determination R² =0.85):

$$\begin{array}{}\text{(17)}& R_{max}/R_{\mathrm{0}}=\mathrm{1.17}{\left(\mathit{\lambda }/L\right)}^{\mathrm{1}/\mathrm{4}}{\left(s/s_{\mathrm{0}}\right)}^{\mathrm{1}/\mathrm{4}}.\end{array}$$

Figure 11 The normalized maximum run-up height, R_max∕R₀ (L∕λ)^1∕4, calculated numerically versus the wave front steepness, s∕s₀, for the same values of A∕h₀ and kh₀ as in Fig. 10. Red solid line is proportional to the “analytically estimated” Eq. (16), while black solid line corresponds to Eq. (17).

5 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 wave front, as was observed during large tsunami events, e.g. the 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 a Riemann wave, while the wave run-up on a plane beach is calculated 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 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.

• It is found analytically that the 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: $R_{max}/R_{\mathrm{0}}={\left(s/s_{\mathrm{0}}\right)}^{\mathrm{0.42}}$.

• 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 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₀ < 1.5), numerical curves of the maximum tsunami run-up height are parallel to the analytical ones, while for larger face front wave steepness (s∕s₀ > 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, which gives the best fit with the data of numerical simulations in general, is $R_{max}/R_{\mathrm{0}}=\mathrm{1.17}{\left(\mathit{\lambda }/L\right)}^{\mathrm{1}/\mathrm{4}}{\left(s/s_{\mathrm{0}}\right)}^{\mathrm{1}/\mathrm{4}}$.

• These results can also be used in tsunami forecasts. Sometimes, in order to save time for tsunami forecasts, 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 tsunami wave can be estimated from the data of the virtual (computed) or real tide-gauge stations and then be used to estimate the 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 to neglect dispersive effects. In this context, it is important to mention the recent work of Larsen and Fuhrman (2019). They used Reynolds-averaged Navier–Stokes (RANS) equations and k–ω model for turbulence closure to simulate the propagation and run-up of positive 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 Madsen et al. (2008), namely that these short waves would have little effect on the overall run-up and inundation of tsunamis (though they found that they could significantly increase the maximum flow velocities).