An Efficient Two-Layer Non-Hydrostatic Model for Investigating Wave Run-Up Phenomena

Abstract

In this paper, we study the maximum run-up of solitary waves on a sloping beach and over a reef through a non-hydrostatic model. We do a modification on the non-hydrostatic model derived by Stelling and Zijlema. The model is approximated by resolving the vertical fluid depth into two-layer system. In contrast to the two-layer model proposed by Stelling, here, we have a block of a tridiagonal matrix for the hydrodynamic pressure. The equations are then solved by applying a staggered finite volume method with predictor-corrector step. For validation, several test cases are presented. The first test is simulating the propagation of solitary waves over a flat bottom. Good results in amplitude and shape preservation are obtained. Furthermore, run-up simulations are conducted for solitary waves climbing up a sloping beach, following the experimental set-up by Synolakis. In this case, two simulations are performed with solitary waves of small and large amplitude. Again, good agreements are obtained, especially for the prediction of run-up height. Moreover, we validate our numerical scheme for wave run-up simulation over a reef, and the result confirms the experimental data.

1. Introduction

Wave run-up is the maximum vertical limit of wave upsurge which occurs on a beach or structure above the still water level. In the tsunami mitigation plan, a good estimation of wave run-up is needed to determine the evacuation area to protect people from the tsunami. Many researchers developed empirical formulas to predict the run-up for regular waves [1,2,3]. Analytical formulation of wave run-up was also derived for different types of incident waves based on Carrier and Greenspan’s solution [4,5,6]. Both empirical and analytical formulas are useful for practical application but they are generally limited to a simplified sloping structure. For more complex topography and general type of waves, the numerical approach is necessary to predict the wave run-up.

Several wave models can be used to approximate the wave run-up. Shallow water equations are often used in wave run-up models, especially in tsunami cases. The numerical scheme based on shallow water equations is good enough to simulate wave run-up for long waves [7,8]. However, for dispersive waves, i.e., solitary waves, this model failed to simulate the propagation with the correct speed. Thus, another model is needed in dispersive wave modelling.

There are two approaches to model the dispersive waves, i.e., Boussinesq-type and non-hydrostatic approaches. The Boussinesq-type model was constructed for irrotational and incompressible fluid flow without friction [9,10,11]. Nevertheless, the difficulty in dealing with the Boussinesq model is due to the complexity of non-linearity and dispersion terms. An alternative model that involves non-linearity and dispersion is the depth integrated non-hydrostatic model which was firstly developed by [12] to solve Euler equations with total pressure. Then, this total pressure was separated by [13] and [14] into hydrostatic and hydrodynamic pressure to improve the efficiency of the computational time.

Stelling and Zijlema (2003) [15] in SWASH and Magdalena [16] proposed an alternate depth-integrated formulation, which introduces non-hydrostatic pressure and vertical velocity terms in the non-linear shallow water equations (NSWE) to describe weakly dispersive waves. The paper presents some validations of the scheme by making a comparison between the numerical results and analytical solution, also with the experimental data. To obtain a good agreement with both analytical and experimental results, it is necessary to use two or more layers. Based on the scheme introduced by Stelling and Zijlema (2003), we propose a modification of two-layer non-hydrostatic approach.

Here, we study the performance of the modified two-layer scheme. It can handle wave frequency correctly including short waves up to kd ≈ 6. For the application of two-layer fluid system, this scheme can simulate soliton propagation with undisturbed shape. Therefore, we can conclude that the modified two-layer scheme can balance the non-linearity and dispersion effects. This is the most important test case for non-linear and dispersive model. In this research, first, we construct a momentum conservative numerical scheme without considering the dispersion effect using a staggered finite volume method. This method, which was initially investigated by [17] to solve differential equations, was further extended and implemented by many researchers [15,18,19]. Working with a staggered finite volume gives us more advantages than working in a collocated grid because we do not end up with the Riemann problem that is relatively difficult to solve, as seen in the NHWAVE proposed by [20]. Besides, this scheme could give accurate result for subcritical flow [21] and relatively simple even for three-dimensional extent [14,22,23]. Secondly, to derive the model with dispersion, we propose a non-hydrostatic approach that incorporates hydrodynamic pressure. By taking the hydrodynamic pressure into account, we consider inhomogeneity in vertical direction. This two-layer approach has a tridiagonal matrix for pressure, whereas the standard two-layer approach in SWASH by Stelling and Zijlema (2003) has a pentadiagonal matrix. Dealing with a tridiagonal matrix allows us to implement the Thomas algorithm, which is faster than a pentadiagonal matrix that needs to be solved using an iterative method. Thus, the computational cost of this two-layer approach is comparable with that of the one-layer approach [24], but the performance is comparable with the two-layer approach as conducted in Stelling and Zijlema (2003). We can argue that this two-layer approach is efficient.

In this paper, we will see how wide the range of problems that can be covered by the numerical scheme is. A series of test cases will be conducted to validate and verify the two-layer non-hydrostatic numerical scheme. First, we will simulate the propagation of solitary waves over a flat bottom topography. Then, the topography will be extended to a sloping beach with and without a reef. The results will be compared to the data measured in laboratory experiments. Finally, the numerical scheme is implemented to simulate the 2004 Indian Ocean tsunami run-up on the west coast of Aceh, Indonesia.

2. The Two-Layer Non-Hydrostatic Model

In this section, we explain briefly the derivation of the non-hydrostatic model using the two-layer approach to study wave run-up over a sloping beach. First, we consider our observation domain as a two-layer fluid system, upper and lower layer domain, that are indexed by subscript 1 and 2, respectively. The sketch of this problem is shown in Figure 1. Here, $d\left(x\right)$ denotes the bottom topography, while ${\eta }_k$ denotes the free surface elevation. The interface between the two layers is located at $z_{\alpha }$. The water thickness for each layer is $h_k$ for $k=1,2$. The horizontal and vertical velocity is denoted by $u_k$ and $w_k$, respectively.

Our derivation based on the Euler equations model in terms of the velocity $(u,w)$ read as

$$\begin{array}{c}\hfill \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0,\end{array}$$

(1)

$$\begin{array}{c}\hfill \frac{\partial u}{\partial t}+\frac{\partial u^2}{\partial x}+\frac{\partial uw}{\partial z}+g\frac{\partial \eta }{\partial x}+\frac{\partial q}{\partial x}=0,\end{array}$$

(2)

$$\begin{array}{c}\hfill \frac{\partial w}{\partial t}+\frac{\partial uw}{\partial x}+\frac{\partial w^2}{\partial z}+\frac{\partial q}{\partial x}=0,\end{array}$$

(3)

where q is the essential term in this model that is our hydrodynamic pressure. Respectively, the kinematic free surface and bottom boundary conditions for the flow read as:

$$\begin{array}{c}\hfill w_{\eta }=\frac{\partial \eta }{\partial t}+u_{\eta }\frac{\partial \eta }{\partial x},\end{array}$$

(4)

$$\begin{array}{c}\hfill w_d=-u_d\frac{\partial d}{\partial x}.\end{array}$$

(5)

The subscripts $\eta$ and d denote the position where quantity evaluated. Based on the Euler equations above, we integrate Equations (1) and (2) from the bottom to the interface. By implementing the kinematic boundary condition (5) and Leibniz’s rule, we obtain

$$\begin{array}{c}\hfill \frac{\partial h_2{u}_2}{\partial x}-u_{z_{\alpha }}\frac{\partial z_{\alpha }}{\partial x}+w_{z_{\alpha }}=0,\end{array}$$

(6)

$$\begin{array}{c}\hfill \frac{\partial u_2}{\partial t}+u_2\frac{\partial u_2}{\partial x}+g\frac{\partial \eta }{\partial x}+\frac{1}{2}\frac{\partial }{\partial x}(q_d+q_{z_{\alpha }})-\frac{1}{h_2}{q}_d\frac{\partial d}{\partial x}-\frac{1}{h_2}{q}_{z_{\alpha }}\frac{\partial z_{\alpha }}{\partial x}.\end{array}$$

(7)

Similarly, for the upper layer, incorporating the kinematic boundary condition (4) and Equation (6) into the integration of Equations (1) and (2) from the interface to the surface, we obtain

$$\frac{\partial u_1}{\partial t}+u_1\frac{\partial u_1}{\partial x}+g\frac{\partial \eta }{\partial x}+\frac{1}{h_1}=\frac{1}{2}(q_{z=-d}+q_{z_{\alpha }})-q_{z=-d}\frac{\partial d}{\partial x}-q_{z_{\alpha }}\frac{\partial z_{\alpha }}{\partial x}.$$

(8)

Under the assumption that $q_{z=\eta }=0$, then the momentum Equation (8) becomes

$$\frac{\partial u_1}{\partial t}+u_1\frac{\partial u_1}{\partial x}+g\frac{\partial \eta }{\partial x}+\frac{1}{2}\frac{\partial q}{\partial x}+\frac{1}{2h_1}\frac{\partial (\eta -z_{\alpha })}{\partial x}=0.$$

(9)

In this non-hydrostatic approach, vertical homogeneities are adopted. The relations between vertical velocity and hydrodynamic pressure q for the two layers come from the linearized momentum equation with the Keller-box scheme read as

$$\begin{array}{c}\hfill \frac{1}{2}\left(\frac{\partial w_1}{\partial t}+\frac{\partial w_2}{\partial t}\right)+\frac{q_{\eta }-q_{z_{\alpha }}}{h_2}=0,\end{array}$$

(10)

$$\begin{array}{c}\hfill \frac{1}{2}\left(\frac{\partial w_0}{\partial t}+\frac{\partial w_1}{\partial t}\right)+\frac{q_{z_{\alpha }}-q_d}{h_1}=0.\end{array}$$

(11)

The closure relations needed to calculate $w_1$ and $w_2$ come from the mass conservation for both layers as follows.

$$\begin{array}{c}\hfill (w_1-w_2)+h_1{\left(u_1\right)}_x=0,\end{array}$$

(12)

$$\begin{array}{c}\hfill (w_2-0)+h_2{\left(u_2\right)}_x=0.\end{array}$$

(13)

To eliminate the vertical velocity and the pressure at the bottom from the set of equations, first, Equation (12) is multiplied by 2, then added to Equation (13):

$$(w_1+w_2)+h_1{\left(u_1\right)}_x+2h_2{u}_2=0.$$

(14)

From now on, we treat $w_{12}=w_1+w_2$ as a new variable.

Here, we assume the fluid density profile increases with depth only in the upper part of a fluid column. In the lower part, fluid density value is nearly constant, as commonly observed in nature. Therefore, the hydrodynamic pressure q is assumed to be a constant function throughout the lower layer. In this two-layer formulation, we denote $q_d=q_{z_{\alpha }}\equiv q$, so that Equation (7) is simplified to

$$\frac{\partial u_2}{\partial t}+u_2\frac{\partial u_2}{\partial x}+g\frac{\partial \eta }{\partial x}+\frac{\partial q}{\partial x}-\frac{1}{h_2}q\frac{\partial d}{\partial x}-\frac{1}{h_2}q\frac{\partial z_{\alpha }}{\partial x}.$$

(15)

Finally, the full set of non-hydrostatic model using two-layer approach is

$$\begin{array}{c}\hfill \frac{\partial \eta }{\partial t}+\frac{\partial h_2{u}_2}{\partial x}+\frac{\partial h_1{u}_1}{\partial x}=0,\end{array}$$

(16)

$$\begin{array}{c}\hfill \frac{\partial u_1}{\partial t}+u_1\frac{\partial u_1}{\partial x}+g\frac{\partial \eta }{\partial x}+\frac{1}{2}\frac{\partial q}{\partial x}+\frac{1}{2h_1}\frac{\partial (\eta -z_{\alpha })}{\partial x}=0,\end{array}$$

(17)

$$\begin{array}{c}\hfill \frac{\partial u_2}{\partial t}+u_2\frac{\partial u_2}{\partial x}+g\frac{\partial \eta }{\partial x}+\frac{\partial q}{\partial x}-\frac{1}{h_2}q\frac{\partial d}{\partial x}-\frac{1}{h_2}q\frac{\partial z_{\alpha }}{\partial x},\end{array}$$

(18)

$$\begin{array}{c}\hfill \frac{1}{2}\frac{\partial w_{12}}{\partial t}+\frac{q_{\eta }-q_{z_{\alpha }}}{h_2}=0,\end{array}$$

(19)

$$\begin{array}{c}\hfill w_{12}+h_1{\left(u_1\right)}_x+2h_2{u}_2=0.\end{array}$$

(20)

The system consists of five equations with five unknown variables $\eta ,u_1,u_2,w_{12}$, and $q$. This gives us an advantage in the computation where we do not need to compute Equation (11) anymore.

3. Dispersion Relation

In this section, we derive the dispersion relation of our linearized two-layer non-hydrostatic model over a constant depth $d=d_1+d_2$ where $d_1=\alpha d$ and $d_2=(1-\alpha )d$. To derive the dispersion relation of the proposed model, first, we assume that the wave is a periodic wave with a small amplitude in the form of

$$\begin{array}{c}\hfill \eta (x,t)=\overline{\eta }{e}^{i(kx-\omega t)},\end{array}$$

(21)

$$\begin{array}{c}\hfill u_1(x,t)=\overline{u_1}{e}^{i(kx-\omega t)},\end{array}$$

(22)

$$\begin{array}{c}\hfill u_2(x,t)=\overline{u_2}{e}^{i(kx-\omega t)},\end{array}$$

(23)

$$\begin{array}{c}\hfill w_{12}(x,t)=\overline{w_{12}}{e}^{i(kx-\omega t)},\end{array}$$

(24)

$$\begin{array}{c}\hfill q(x,t)=\overline{q}{e}^{i(kx-\omega t)}.\end{array}$$

(25)

Notation $\omega$ denotes wave celerity, and k is the wavenumber. Substituting Equations (21)–(25) into the linearized governing equations written as

$$\begin{array}{c}\hfill \frac{\partial \eta }{\partial x}+d_2\frac{\partial u_2}{\partial x}+d_1\frac{\partial u_1}{\partial x}=0,\end{array}$$

(26)

$$\begin{array}{c}\hfill \frac{\partial u_1}{\partial t}+g\frac{\partial \eta }{\partial x}+\frac{1}{2}\frac{\partial q}{\partial x}=0,\end{array}$$

(27)

$$\begin{array}{c}\hfill \frac{\partial u_2}{\partial t}+g\frac{\partial \eta }{\partial x}+\frac{\partial q}{\partial x}=0,\end{array}$$

(28)

$$\begin{array}{c}\hfill w_{12}+2d_2\frac{\partial u_2}{\partial x}+d_1\frac{\partial u_1}{\partial x}=0,\end{array}$$

(29)

$$\begin{array}{c}\hfill \frac{1}{2}\frac{\partial w_{12}}{\partial t}-\frac{q}{d_2}=0,\end{array}$$

(30)

we obtain the following dispersion relation,

$$\frac{{\omega }^2}{g{k}^2}=\frac{d_1{d}_2^2{k}^2+4d_1+4d_2}{d_1{d}_2{k}^2+4d_2^2{k}^2+4}.$$

(31)

To check the accuracy of our dispersion relation, we compare it with the dispersion relation proposed by [10] and [25]. Nwogu [10] obtained the dispersion relation of Boussinesq-type model written as

$$c^2=gd\frac{1-({\alpha }_N+1/3){\left(kd\right)}^2}{1-{\alpha }_N{\left(kd\right)}^2},$$

(32)

where ${\alpha }_N$ is a parameter that was adjusted to be equal to −1/3 in order to get the result closer to the exact dispersion relation. The linear dispersion relation of the hybrid system obtained by Bai and Cheung wrote as

$$c^2=gd\frac{1+1/16{\left(kd\right)}^2}{1+(3/16+1/4\beta ){\left(kd\right)}^2},$$

(33)

where $\beta$ is a free parameter.

For the comparison between the exact dispersion relation with the ones obtained using our non-hydrostatic model, Boussinesq-type model in [10] and hybrid system in [25], we use $d_1=0.5d$ and $d_2=0.5d$. The comparison is shown in Figure 2. The figure shows that our numerical dispersion relation is in a good agreement with the exact solution. Moreover, the accuracy of our result is up to $kd\simeq 10$ with 1% error.

4. Numerical Formulation

In this section, we implement a staggered finite volume method with predictor-corrector step to solve our governing equations numerically. The predictor step is used to calculate the governing equations with neglecting the hydrodynamic pressure term. In the next step, we correct the velocity by including hydrodynamic pressure and considering inhomogeneity in the vertical direction. We use all values in the hydrostatic part as the predictor in the correction step.

4.1. Predictor Part

Following the staggered momentum conservative scheme for the NSWE [18,19], here, we implement the analogous scheme for our non-hydrostatic model. Consider an observation domain $[0,L]$. Applying a staggered partition, then our grid points are $x_{\frac{1}{2}}=0,x_1,x_{3/2},x_2,\dots ,x_{Nx}$, $x_{Nx+\frac{1}{2}}=L$. The rules of implementing staggered arrangement, the values of ${\eta }_1,{\eta }_2$ are computed using mass conservation (26) at every full grid points $x_i=i\Delta x$, with $i=1,2,\dots ,Nx$. Whereas velocities $u_1,u_2$ are computed using momentum equation ((27), (28)) at every staggered grid points $x_{i+\frac{1}{2}}=(i+1/2)\Delta x$, with $i=0,1,2,\dots ,Nx$, see Figure 3. Since $d_1,d_2$, and $d\left(x\right)$ are functions of x, once ${\eta }_1$, and ${\eta }_2$ are known, $h_1$, and $h_2$ can be obtained.

Then, the approximated Equations of (26)–(28) are:

$$\begin{array}{c}\hfill \frac{{\eta }_i^{n+1}-{\eta }_i^n}{\Delta t}+\frac{{}^{*}h{u}_{i+1/2,1}^n-{}^{*}h{u}_{i-1/2,1}^n}{\Delta x}=0,\end{array}$$

(34)

$$\begin{array}{c}\hfill \frac{u_{i+1/2,1}^{n+1}-u_{i+1/2,1}^n}{\Delta t}+g\frac{{\eta }_{i+1,1}^{n+1}-{\eta }_{i,1}^{n+1}}{\Delta x}+{\left(u{u}_x\right)}_{i+1/2,1}=0,\end{array}$$

(35)

$$\begin{array}{c}\hfill \frac{u_{i+1/2,2}^{n+1}-u_{i+1/2,2}^n}{\Delta t}+g\frac{{\eta }_{i+1,1}^{n+1}-{\eta }_{i,1}^{n+1}}{\Delta x}+{\left(u{u}_x\right)}_{i+1/2,2}=0,\end{array}$$

(36)

where the asterisk denotes the unknown values that we need to approximate. To estimate the values of water depth $h_{i+1/2,k}$, for $k=1,2$ at half grid points, we use first order upwind method read as:

$${}^{*}{h}_{i+\frac{1}{2},k}=\left\{\begin{array}{cc}{h}_{i,k},\hfill & \mathrm{if}{u}_{i+1/2,k}\ge 0,\hfill \\ h_{i+1,k},\hfill & \mathrm{if}{u}_{i+1/2,k}<0.\hfill \end{array}\right.$$

(37)

The discretization of the advection term written as

$${\left(u{u}_x\right)}_{i+\frac{1}{2},k}=\frac{1}{{\overline{h}}_{i+\frac{1}{2},k}}\left(\frac{{\overline{q}}_{i+1,k}{}^{*}{u}_{i+1,k}-{\overline{q}}_{i,k}{}^{*}{u}_{i,k}}{\Delta x}-u_{i+\frac{1}{2},k}\frac{{\overline{q}}_{i+1,k}-{\overline{q}}_{i,k}}{\Delta x}\right),$$

(38)

whereas

$${\overline{h}}_{i+\frac{1}{2},k}=\frac{1}{2}(h_{i,k}+h_{i+1,k}),{\overline{q}}_{i,k}=\frac{1}{2}(q_{i+\frac{1}{2},k}+q_{i-\frac{1}{2},k}),q_{i+\frac{1}{2},k}={}^{*}{h}_{i+\frac{1}{2},k}{u}_{i+\frac{1}{2},k},$$

with the following upwind approximation of $u_{i,k}$ for both layers $k=1,2$,

$${}^{*}{u}_{i,k}=\left\{\begin{array}{cc}{u}_{i-\frac{1}{2},k},\hfill & \mathrm{if}{\overline{q}}_{i,k}\ge 0,\hfill \\ u_{i+\frac{1}{2},k},\hfill & \mathrm{if}{\overline{q}}_{i,k}<0.\hfill \end{array}\right.$$

(39)

4.2. Corrector Part

In this step, we correct the values of $\eta ,u_1$, and $u_2$ by including the hydrodynamic pressure term in the model. The staggered grid arrangement of the new calculated variables is depicted in Figure 4. The pressure variable q and vertical velocities $w_1$ and $w_2$ are calculated at the full grid points. Here, we consider the values of $\eta ,u_1$, and $u_2$ that we have obtained from the hydrostatic model as the predictor values. Those are denoted by superscript *. Taking the assumption that hydrodynamic pressure is zero at the surface, and it is constant throughout the lower layer, we only need one vector array to store q values. It is less complicated compared to the two-layer approach proposed by [26], where two vector arrays for q are needed.

In this non-hydrostatic step, we begin with Equation (20) and its discrete form which is

$$w_{i,12}^{n+1}+2h_{i,2}^{n+1}\frac{u_{i+\frac{1}{2},2}^{n+1}-u_{i-\frac{1}{2},2}^{n+1}}{\Delta x}+h_{i,1}^{n+1}\frac{u_{i+\frac{1}{2},1}^{n+1}-u_{i-\frac{1}{2},1}^{n+1}}{\Delta x}=0.$$

(40)

We substitute the values of $u_1,u_2$ and $w_{12}$ in the correction step that are given by:

$$\begin{array}{ccc}\hfill u_{i+1/2,1}^{n+1}& =& u_{i+1/2,1}^{*}-\Delta t\frac{q_{i+1}^{n+1}-q_i^{n+1}}{2\Delta x}\hfill \\ & & -\frac{1}{2h_{i,1}}\frac{({\eta }_{i+1}^{n+1}-d_{i+1,1}+{\eta }_{i+1,2}^{n+1})-({\eta }_{i,1}^{n+1}-d_{i,1}+{\eta }_{i,2}^{n+1})}{\Delta x}{p}_i,\hfill \end{array}$$

(41)

$$\begin{array}{ccc}\hfill u_{i+1/2,2}^{n+1}& =& u_{i+1/2,2}^{*}-\Delta t\frac{q_{i+1}^{n+1}-q_i^{n+1}}{\Delta x},\hfill \end{array}$$

(42)

$$\begin{array}{ccc}\hfill w_{i,12}^{n+1}& =& w_{i,12}^{*}+\Delta t\frac{q_i^{n+1}}{h_{i,2}^{n+1}},\hfill \end{array}$$

(43)

to get a tridiagonal system of equations for $q_i$. Then, we solve the tridiagonal system by means of Thomas Algorithm.

Related to [15,22,25,26]’s work, the performance of our scheme is comparable to the two-layer approach conducted by [15] with the computational cost is more efficient and comparable with one-layer approximation of [15]. In the standard two-layer approach, we have a pentadiagonal matrix, and here we have a tridiagonal matrix. In comparison, solving a standard 6 x 6 matrix with pentadiagonal linear system in MATLAB with a processor Intel(R) Core(TM) i7-8565U CPU @1.80GHz 1.99GHz will take 0.002561 seconds, and for the tridiagonal system, we only need 0.000835 seconds. Therefore, we can conclude that our numerical scheme is more efficient than the two-layer approach proposed by [15].

5. Numerical Simulations

From the previous section, we achieve the dispersion relation of the two-layer non-hydrostatic model and it gives us a good sign that our model is able to simulate non-linear and dispersive wave phenomena. Here, we shall present some benchmark tests to demonstrate the capability of our model in describing non-breaking and breaking waves evolution. First, we start with the numerical experiments of solitary wave propagation in a flat bottom channel. This test examines whether our numerical model has a balance between dispersion and nonlinearity or not. The second test case is a comparison between wave run-up calculation and experimental data for solitary waves run up on a sloping beach. Next, we show the capability of the proposed model with its wet-dry procedure in describing wave propagation over a dry area and holding a breaking wave by comparing the numerical results to the experimental data. This test case is also used to investigate the appearance of undular bore and its propagation, as well as wave run-up over a fringing reef. In the end, we implemented the numerical model to simulated the 2004 Indian Ocean tsunami run-up using simplified Aceh bathymetry.

5.1. Solitary Waves Propagate over a Flat Bottom

In this subsection, we present a benchmark test for solitary wave propagation over a flat bottom. This is a benchmark test for the dispersive and nonlinear waves. For computation, a channel with $L=1000$ m long and $d=10$ m water depth is considered. The left and right absorbing boundaries are applied. The initial wave is located at $x=X_1=200$ and has the following surface profile

$$\begin{array}{ccc}\hfill \eta (x,0)& =& Asec{h}^2(xi/2(x-200)),\hfill \end{array}$$

(44)

$$\begin{array}{ccc}\hfill u_1(x,0)& =& \sqrt{g/d}\eta (x,0),\hfill \end{array}$$

(45)

$$\begin{array}{ccc}\hfill u_2(x,0)& =& \sqrt{g/d}\eta (x,0),\hfill \end{array}$$

(46)

with $A=1,xi=\sqrt{\frac{3A/d^2}{d+A}}$, and $g=9.81$. For numerical computation, we choose $\Delta t=0.01,\Delta x=1$ and $d_1=d_2=d/2$. We compare the numerical solution to the exact solution written as:

$${\eta }_{exact}(x,t)=Asec{h}^2\left(\frac{xi}{2}(x-200-\sqrt{g(d+A)}t)\right).$$

Owing to the balance between nonlinearity and dispersion, in Figure 5, we can see that the solitary waves travel with the correct speed, constant amplitude, and conserved shape. We plot the numerical result together with the analytical solution and SWASH result and they show a good comparison.

5.2. Solitary Waves Run up on a Sloping Beach

This section presents two benchmark tests to examine the profile of solitary waves propagate to a sloping area for the non-breaking and breaking case. The breaking criteria are based on the ratio of amplitude and water depth $(A/d>0.055)$, as stated in [6], and the wave parameter for the breaking simulations is based on experimental data in that paper. The initial waves have the same profile as used in Synolakis’ experiment [6]. This experiment was conducted at 31.73 m length, 60. 96 cm depth, and 39.97 cm width wave tank. Since the length is much longer than the width, we could use this experimental data to validate a one-dimensional problem. We use a simulation set-up where the bathymetry is a flat bottom connected to a sloping beach with the beach slope is 1:19.85. The initial condition is defined in the flat bottom area as

$$\begin{array}{ccc}\hfill \eta (x,0)& =& Asec{h}^2\left(\gamma (x-x_0)\right),\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill u_1(x,0)& =& \sqrt{g/d_0}\eta (x,0),\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill u_2(x,0)& =& \sqrt{g/d_0}\eta (x,0),\hfill \end{array}$$

(49)

where A is wave height, $d_0$ is undisturbed water depth, g is gravity acceleration, $x_0$ is wave crest position, and $\gamma =\sqrt{\frac{3A}{4d_0}}$. This initial condition causes solitary waves to propagate to the beach on the right side. For this simulation, we set $g=1$ and $d_0=1$. We perform two simulations which are the non-breaking case $(A=0.0185,x_0=38.34)$ and the breaking case $(A/d_0=0.3)$. We choose $\Delta x=0.3$ and $\Delta t=0.003$ to handle the moving boundary at the shoreline. In this simulation, the wet and dry procedure is applied.

We will start with a small normalized wave height $A/d_0=0.0185$. The result is summarized in Figure 6 which shows the comparison between the experimental measurements (dots), the hydrostatic model results [8] (dash line), and the two-layer non-hydrostatic model results (solid line). The comparison is noticed at four different times, at $t=40$, $t=50$, $t=60$, and $t=70$. Both of the numerical results are similar and both of them are in a good agreement with the experimental data.

The next simulation is performed for larger normalized wave height, i.e, $A/d_0=0.3$. Here, the surface profiles are captured at $t=15$, $t=20$, $t=25$, and $t=30$. The comparison between laboratory experiment, the hydrostatic model, and the two-layer non-hydrostatic simulation is shown in Figure 7.

As shown in Figure 7, the hydrostatic model produces wave speeds that are faster than the experiment results. Implementing the non-hydrostatic model will produce better propagation speed than the hydrostatic model. Error comparison between hydrostatic and non-hydrostatic model at $T=15s$ and $T=20s$ is shown in

Table 1. Error for breaking solitary wave ($A/d_0=0.3$).

Time	RMSE
	Hydrostatic Model	Non-Hydrostatic Model
15 s	5.11%	2.63%
20 s	5.19%	3.46%

.

Next, we also compare our numerical model with BoussClaw numerical model as shown in Figure 8. BoussClaw is an extension of GeoClaw (part of Clawpack developed by LeVeque [27]) which solves the Boussinesq-type model. Both numerical models are in a good agreement with the experimental data, yet, the non-hydrostatic model produces a better approximation.

5.3. Run-up Height of Solitary Waves

In this section, the two-layer non-hydrostatic scheme will be implemented to measure the maximum run-up of solitary waves over a sloping beach. We conduct several simulations with different depth and wave height, using the same bathymetry and initial condition as shown in the previous simulation. The value $x_0$ in initial condition changes with different wave height, because solitary waves with different height have different “wave length”. The value of $x_0$ is chosen so that the solitary hump is on the flat bottom area. In this simulation we choose $x_0=x_s+L$ as suggested by [6], where $L=\frac{1}{\gamma }{cosh}^{-1}\left(\sqrt{20}\right)$ is the distance of wave crest from the toe of the beach. It is to assure that all the solitary waves propagate with the same relative distance to the beach. The experimental data in [6] were used to be compared with our model. The water depth ranged from 6.25 cm to 38.32 cm and the wave height ratio $A/d_0$ ranged from 0.005 to 0.633. As stated in [6], the wave breaks during run-up when $A/d_0>0.055$. Run-up height for different wave height is presented in Figure 9. Compared to the experimental data, the run-up height produced by our scheme is good enough (the average error is less than 12%) for solitary waves with the non-breaking case. For solitary waves with the breaking case, run-up height that is produced by our scheme is higher than experimental data. Nevertheless, from our numerical results, we can conclude that the relative run-up height increases as the relative wave height increases, which is confirmed by the experimental data.

5.4. Solitary Waves Propagation over a Sloping Beach with Reefs

Here, we implement our numerical scheme to simulate wave run-up over a sloping beach with a fringing reef. This simulation offers a challenge for researchers who are working on numerical mathematics field because the existence of reef over a sloping beach makes the surf-zone processes more complex. The complexities are in describing the moving waterlines in the dry area and the breaking waves phenomenon. Therefore, a complete dispersive and nonlinear numerical scheme is not enough to simulate this case, so we need to include a correct wet-dry procedure.

In O.H. Hinsdale Wave Research Laboratory of Oregon State University, [28] has investigated the transformation of solitary waves over a fringing reef under wet and dry flat conditions. He has done several experiments with different set-ups. Here, we choose two set-ups with different reef configurations and scales to demonstrate and validate our numerical model. For the first test case, the fringing reef is in a 40 m long flume with 1.0 m water depth over a sloping slope 1/5. For the computation, we choose the spatial grid and time step as $\Delta x=0.2$ and $\Delta t=0.01s$, respectively. The initial condition is calm water. The left boundary is a solitary wave with $A/h=0.3$, and the right boundary is a vertical wall as in the experiment. A good agreement given by the comparison of our numerical results and the measured surface elevation along the flume is shown in Figure 10.

The second test involves a fringing reef with a 1/12 slope and a 0.2 m high reef crest in a 100 m long and 2.56 m deep flume. The grid size and time step are $\Delta x=0.2$ and $\Delta t=0.01$. The comparison between the computed and measured surface elevations for $A/h=0.3$ is presented in Figure 11. When the solitary waves propagate over the slope starting at $x=25.9$, the waves begin to shoal and then we see the appearance of a bore on the reef crest as the bore propagates and attenuates over time. This test case demonstrates the balance between dispersion and nonlinearity with the hydraulic jump near the shock.

5.5. Run-up Simulation Using Simplified Aceh Bathymetry

Finally, this numerical model is implemented to simulate wave run-up on simplified Aceh bathymetry. Aceh is one of the worst affected areas by the 2004 Indian Ocean tsunami. We use the same bathymetry data as mentioned in [29] which located on the west coast of Aceh. Its location is about 100 km from the epicenter of the earthquake. We simplify the bathymetry using linear piecewise function as follows.

$$b\left(x\right)=\left\{\begin{array}{cc}-d_0=-1041\hfill & x<-80000\hfill \\ \frac{989}{50000}x+541.4\hfill & -80000\le x<-30000\hfill \\ -52\hfill & -30\le x<-18000\hfill \\ \frac{63}{22000}x\hfill & x>-18000.\hfill \end{array}\right.$$

(50)

Based on the report in [30], we use the leading depression N-waves with the wave amplitude of 0.6 m as initial waves, defined as follows.

$$\begin{array}{ccc}\hfill \eta (x,0)& =& \frac{0.6}{N}exp\left(-\frac{{(x+100)}^2}{{35000}^2}\right)\left(-\frac{2(x+100)}{{35000}^2}\right),\hfill \end{array}$$

(51)

$$\begin{array}{ccc}\hfill u_1(x,0)& =& \sqrt{g/d_0}\eta (x,0),\hfill \end{array}$$

(52)

$$\begin{array}{ccc}\hfill u_2(x,0)& =& \sqrt{g/d_0}\eta (x,0),\hfill \end{array}$$

(53)

$$\begin{array}{ccc}\hfill N& =& max\left(exp\left(-\frac{{(x+100)}^2}{{35000}^2}\right)\left(-\frac{2(x+100)}{{35000}^2}\right)\right).\hfill \end{array}$$

(54)

Figure 12 shows the bathymetry and the initial waves. The simplified bathymetry is close enough with the real one. The simulation is conducted for $T=5000s$ and the dynamics of moving shoreline are shown in Figure 13. Based on the simulation, the waves will arrive at the coastal area about 46 minutes. According to the tsunami travel time map by NOAA [31], the arrival time of the first wave on the west coast of Aceh is about up to one hour which is consistent with our numerical result. From the simulation, the maximum run-up is 8.02 m. Using Equation (50), the inundation area is about 2.8 km.

6. Conclusions

A two-layer non-hydrostatic model was proposed in this paper. Our numerical dispersion relation confirms the exact linear dispersion relation. A staggered finite volume method was used to approximate the governing equations in the horizontal direction, and a Keller box method resolved the vertical velocity distribution with the non-hydrostatic pressure. The predictor–corrector technique that consists of a hydrostatic and a non-hydrostatic step works well for solving the non-hydrostatic model. Having a tridiagonal matrix structure for pressure make this computation more efficient than the standard two-layer system with the pentadiagonal matrix. The numerical solution provided accurate results for non-linear and dispersive waves in a series of test cases. The numerical model is able to describe wave run-up over a sloping beach. Moreover, It gives very good agreement with the experimental data for non-linear wave transformation and evolution over a sloping beach and fringing reef. The present model accounts for wave run-up over a real bathymetry.