Development of Two-Dimensional Non-Hydrostatic Wave Model Based on Central-Upwind Scheme

In this study, a two-dimensional depth-integrated non-hydrostatic wave model is developed. The model solves the governing equations with hydrostatic and non-hydrostatic pressure separately. The velocities under hydrostatic pressure conditions are firstly obtained and then modified using the biconjugate gradient stabilized method. The hydrostatic front approximation (HFA) method is used to deal with the wave breaking issue, and after the wave breaks, the non-hydrostatic model is transformed into the hydrostatic shallow water model, where the non-hydrostatic pressure and vertical velocity are set to zero. Several analytical solutions and laboratory experiments are used to verify the accuracy and robustness of the developed model. In general, the numerical simulations are in good agreement with the theoretical or experimental results, which indicates that the model is able to simulate large-scale wave motions in practical engineering applications.


Introduction
In the past several decades, due to the rapid development of coastal zones as well as the occurrence of a large number of coastal natural disasters (such as wind storm and tsunami), the research on wave propagations in coastal areas has attracted more attention all over the world. With the development of computer technology and the better understanding of wave mechanisms, the wave model is widely used to simulate the evolution of nonlinear and dispersive waves from deep water to shallow water. A high-precision wave model can accurately predict the wave motions in the coastal areas, so as to effectively understand the change of coastal topography, protect coastal buildings and reduce the loss of life and property caused by coastal natural disasters. It is, therefore, of great practical significance to establish a high-precision and accurate wave model. The depth-integrated models based on Boussinesq-type equations (BTEs) are widely used to simulate wave motions [1][2][3][4][5][6][7]. However, since there are several high-order partial derivative terms included, the discretization of the BTEs is very complex, and the computational cost is expensive [8]. Furthermore, the BTEs are derived under the assumption of no rotation and no viscosity and are unable to simulate the interaction of waves with rotational currents [9]. A semi-empirical method or additional terms may be needed to deal with the issue of wave breaking [10].
In order to overcome the shortcomings of the two-dimensional depth-integrated wave model, an alternative model, based on the three-dimensional Reynolds-averaged Navier-Stokes (RANS) equation, is proposed to simulate wave motions under various conditions, including deep water (RANS) model. The governing equations including the mass and momentum conservation for the three-dimensional RANS model can be expressed as: where u, v, and w are the mean velocities in the x, y, and z directions shown in Figure 1, respectively; t is time; ρ is the fluid density; g represents the gravitational acceleration. Assume that the Reynolds stresses are dominant over the fluid shear stresses and τ ij (i, j = x, y, z) are the turbulent stresses. P is the mean pressure. The pressure can be divided as the hydrostatic part p and non-hydrostatic part Γ, which is: where p linearly increases with water depth. According to the assumption made in Stelling and Zijlema [9], the non-hydrostatic Γ and vertical velocity w changes linearly with water depth. The depth-averaged non-hydrostatic pressure Γ and vertical velocity w is thus shown as follows: where subscripts s and b mean surface and bottom, respectively. In general, non-hydrostatic pressure at the surface Γ s is set to 0.
Assuming that the bed elevation is fixed (∂z b /∂t = 0), then the kinematic surface and bottom boundary conditions are: where η is the water surface fluctuation; z b is bottom elevation, which is assumed constant for fixed bottom conditions. The relationship between η, h, and z b is shown in Figure 1. By inserting Equation (5) into Equations (2)-(4), integrating the equations along the vertical direction, and introducing Equations (6) to (9) and the turbulent stress at the bottom boundary (τ xz ) b = ρgu u 2 + v 2 n 2 /h 1/3 (τ yz ) b = ρgv u 2 + v 2 n 2 /h 1/3 , the depth-averaged governing equations for the two-dimensional non-hydrostatic model are: ∂hu ∂t ∂hv ∂t where overbar "-" denotes the average quantity, and n is the Manning's coefficient. For convenience, the overbar "-" symbol will be dropped afterward, and u and v represent the averaged velocity in the x and y directions in the following section.
In order to present this numerical method clearly and compactly, Equations (10)- (12) are rewritten in vector forms: ∂q ∂t where q is a vector of the conserved variables; f and g are the flux vectors in the x and y direction, respectively; s is the source term vectors including the bed slope term s h , bed friction term s f , and non-hydrostatic pressure term s p .

Numerical Methods
The calculation of the two-dimensional non-hydrostatic wave model is divided into two steps. The first step is to solve the governing equations for hydrostatic pressure, and the non-hydrostatic pressure is then calculated.

Calculation of the Hydrostatic Pressure Term
The governing equations for hydrostatic pressure means that the source term s in Equation (14) only include the bed slope term s h and bed friction term s f . Therefore, the governing equations are actually the nonlinear shallow water equations (NLSWEs). In this study, the mathematical model is discretized using a finite volume method based on a rectangular mesh with a non-staggered (collocated) grid system, as shown in Figure 2. Integrating Equation (14) over the (i, j) control volume, applying the Green's theorem, and using the second order Runge-Kutta method for the time derivative, the discretization of Equation (14) is: where superscript o is the o-th time step; subscripts i, j are the grid indexes; ∆t is the time step, and the operator D(q i,j ) is defined as: where ∆x and ∆y are the cell lengths in x and y directions, f i+1/2,j and f i−1/2,j are the fluxes at faces (i + 1/2, j) and (i − 1/2, j). g i,j+1/2 and g i,j−1/2 are the fluxes at faces (i, j + 1/2) and (i, j − 1/2). s hi,j and s f i,j represent the source terms evaluated at the cell center.
To preserve the well-balanced property, the nonnegative reconstruction procedure proposed by Liang [30] is used in the present model (see the detail in Appendix A). Moreover, the model adopts a MUSCL interpolation with minmod limiter to achieve second order accuracy in space. Further details of MUSCL interpolation are given in Appendix A. To solve Equation (15), the central-upwind scheme is used to calculate convective fluxes at cell interfaces, which is written as: where superscript L and R represent the reconstructed states at the left and right side of the cell interface, respectively, and a + and a − denote one-sided local wave speeds, which are estimated as: where V is the velocity at the cell interface. For the treatment of wet/dry front, the bed slope term s h and bed friction term s f , one can refer to Appendix A or Wu et al. [29] for more details.

Calculation of the Non-Hydrostatic Pressure Term
The calculation of the non-hydrostatic pressure term s p is to modify the velocity obtained from the hydrostatic pressure calculation. The equations for velocity corrections are: where superscript m represents the results from the hydrostatic calculation. From the above equations, it can be known that the key to obtain the velocity under non-hydrostatic conditions is to calculate the non-hydrostatic pressure (Γ b ) o+1 i,j at the new time step o + 1. In Equation (22) i,j can be calculated according to Equation (9), which is: By re-discretizing Equations (20) and (21) at the cell interfaces, the following equations can be obtained: Assumming that After some mathematical manipulations, the velocity in x direction at interface (i + 1/2, j) and (i − 1/2, j) can be written as The similar procedure is applied to obtain v o+1 i,j+1/2 and v o+1 i,j−1/2 .
Bottom non-hydrostatic pressure at the new time step can be evaluated through Equation (1), and the corresponding implicit discretized equations are: Inserting Equations (23) and (26) into Equation (27), the linear equations (Poisson equations) for Γ b are given below: where the coefficients are: In this study, the biconjugate gradient stabilized method is adopted to solve (Γ b ) o+1 i,j in Equation (28), and (Γ b ) o+1 i,j is then inserted into Equations (20)- (22) to obtain the velocity u o+1 , v o+1 and w o+1 s .

Stability Criterion and Boundary Conditions
Since the current numerical scheme is explicit, its stability is determined by the Courant-Friedrichs-Lewy (CFL) criterion. In order to ensure the stability of the model at each time step, the Courant number N CFL has to be kept less than 0.25 [29]. In this study, the Courant number N CFL is set to be less than 0.20 for all simulations, which is where a is given by: The ghost-cell is used for the boundary conditions. For the open boundary conditions, flow variables in the ghost cells are usually determined by solving a boundary Riemann problem, which is dependent on the local flow regime. For wall boundary conditions, the velocity normal to the boundary and the water surface gradient are both set to zero at the boundary. For wave inlet boundaries, the variables on the ghost cells are set according to the corresponding analytical formula. For no reflecting boundaries, the sponge layer technique introduced by Larsen and Dancy [31] is employed to fully absorb the incident wave. The following damping coefficient function is used: where α s and r s are two parameters. n is the number of grid. As suggested by Chen et al. [32], α s = 2, r s = 0.88 − 0.92 and n = 50 − 100 are used in this study.

Manipulation of Wave Breaking Conditions
To deal with the wave breaking conditions, the model adopts the hydrostatic front approximation (HFA) method proposed by Smit et al. [33]. When the temporal variation rate of water levels ∂η/∂t at cells is larger than α 1 gh (α 1 is an empirical coefficient and set to 0.6 as Smit et al. [33] suggested in this study), the cell is labeled as wave breaking and the hydrostatic computation is implemented. Therefore, for the wave breaking cells, the non-hydrostatic pressure Γ and vertical velocity w are set to zero.

Model Verification
The numerical method presented in the above section is tested with several analytical solutions and laboratory experiments in this section. For all tests, the gravitational acceleration g = 9.8 m/s 2 and water density ρ = 1000 kg/m 3 .

Case 1: Propagation of One-Dimensional Solitary Wave in Constant Depth
The first test case is the solitary wave propagation in constant water depth. Solitary wave propagation is extensively used to validate the dispersion characteristics of Boussinesq and non-hydrostatic numerical models. The computational domain is 1000 m long and 10 m deep with a smooth and flat bottom. The solitary wave is initially at x 0 = 200 m with wave height H = 2.0 m. According to the potential flow theorem, the corresponding water level, horizontal and vertical velocities can be obtained, which are [9]: where d is the still water depth, and c = g(H + d) is the celerity for the solitary wave. The computational domain is discretized by 2000 × 3 uniform and rectangular grids with ∆x = ∆y = 0.5 m. The time step ∆t is adjusted during the simulation based on the Courant number, which is taken as 0.2. The simulation time is 50 s. Equations (33)- (35) give the initial water level and velocities at t = 0. The left and right sides of the computation domain are set as no wave reflection boundaries. Figure 3 shows the comparisons of simulated results (surface elevations and velocities at the surface) and analytic solutions at t = 0, 25 and 50 s. The simulations agree well with the analytic results, which means that the developed model can successfully capture the wave feature, and the numerical diffusion almost has no significant impact on the simulations.
In order to evaluate the order of the accuracy of the model, four different grids with the mesh size of ∆x = 0.5 m, 1.0 m, 2.0 m and 4.0 m are considered in this study. Figure 4 gives L 2 error norms in η at t = 50 s as a function of grid size ∆x, which is defined as: where N is the total number of computational cells, η a is the analytical solution. It can be found that the numerical errors are decreased by increasing the number of grids and the accuracy of the present model is close to second-order.

Case 2: Runup of One-Dimensional Solitary Wave
The experiment for runup of one-dimensional solitary wave conducted by Synolakis [34] is used widely to test the capability of numerical models on dealing with wave breaking and wetting-drying. Figure 5 provides the setup and bottom geometry of the experiment, where H is the wave height of the solitary wave, d is the still water depth, X 1 is the location of wave crest and β (= 2.884 • ) is the beach slope. The computational domain x ∈ [−30, 100] with the corresponding bathymetry: The still water depth d is 1 m. At t = 0, the initial surface elevation is as follows: where γ = √ 3H/4d and X 1 = √ 4d/3H arcosh( √ 1/0.05). The initial horizontal velocity u is: The right-side boundary condition is free outflow, and the Manning's coefficient n for the whole computational domain is designated as 0.01 s/m 1/3 . The grid size is 0.1 × 0.1 m, and the time step is adaptive. In this study, the numerical simulations using the non-hydrostatic model and shallow water model are compared with the experimental results at H/d = 0.3. Figure 6 presents the comparisons of the surface elevation variations at different dimensionless times. At the initial stage, when the incident wave propagates to the shoal, the solitary wave gradually becomes asymmetric, i.e., the wave front steepens, and the wave leeward turns into milder, and the wave height increases as the water depth decreases. It can be seen from Figure 6 that the wave front calculated by the shallow water equation (SWE) is steeper than the experimental results, which is markedly different from the results obtained from the non-hydrostatic wave model (NHW). On the other hand, the NHW method can simulate dispersion characteristics of wave rather well, and the simulation result is in good agreement with the experimental measurements. The main reason for obvious differences in the shallow water model is that it does not contain high-order diffusion terms, leading to its inability to accurately simulate the dispersion effect for the wave motions. When t(g/d) 1/2 = 25, because of the shoaling effect, the solitary wave breaks and keeps moving up to the beach. At t(g/d) 1/2 = 40, the wave climbs to the highest point and then begins to run down. During rundown periods, the flow on the slope becomes supercritical. Since the flow is subcritical over constant depth regions, a hydraulic jump will occur at the toe of the slope. From Figure 6, it can be seen that after wave breaking, the numerical results from the SWE and NHW models are similar. This is because the hydrostatic front approximation (HFA) method is used on wave breaking issues. The essence of the method is that when the wave breaks, the non-hydrostatic model is transformed into the hydrostatic shallow water model, and the non-hydrostatic pressure Γ and vertical velocity w are set to zero. In general, the non-hydrostatic model presented in this study can successfully simulate the breaking solitary wave runup and rundown along a beach.

Case 3: Propagation of Solitary Wave over Fringing Reef
Fringing reef is a kind of coral reef widely existing in tropical and subtropical areas. It connects directly to the coastline and stretches out hundreds to several kilometers. An ideal reef consists of fore reef and reef flat. When the wave passes through the reef, under the action of the slope in front of the reef and the horizontal reef flat, the wave breaks and decays. Therefore, fringing reef plays a protective role in the coastline.
In this study, the propagation and deformation of the solitary wave over the reef is used to examine the ability of the model to deal with wave breaking and wetting and drying. Roeber and Cheung [35] carried out a laboratory experiment to study the propagation and deformation of the solitary wave over the reef. The experimental setup and corresponding geometry is shown in Figure 7. The water is 2.5 m deep and the reef flat is submerged by a depth of 0.14 m. The reef slope is 1:12 with a crest 0.065 m above the water level. The wave height of the input solitary wave is 0.75 m. When the solitary wave propagates into the shallow water, a plunging breaker is developed over the reef. Several water level gauges (red dots on Figure 7) are employed to provide temporal and spatial water level variations.
The numerical model uses a grid size ∆x = 0.05 m and a Manning coefficient of n = 0.016 s/m 1/3 to represent the finished concrete surface of the reef model. The time step ∆t is adaptive to ensure the model stability. As previous examples, the HFA method is used on wave breaking issues. Figure 8 provides the comparison between numerical results and measurements. At t(g/d) 1/2 = 63 when the solitary wave reaches the reef slope, its wave height increases, and its wave form gradually become asymmetric, i.e., the wave front becomes steeper. At t(g/d) 1/2 = 67, the wave breaks. When t(g/d) 1/2 = 70.5, the breaking wave enters the reef flat and leads to a hydraulic jump. At t(g/d) 1/2 = 83, the wave reaches the right boundary and reflects back to the left side. The reflection wave then jumps over the reef crest and propagates to the upstream. From the comparisons at different times, the computational results are in acceptable agreement with the measurements, which proves that the developed model is able to accurately simulate the processes of solitary wave propagation, breaking and reflection. Figure 9 shows good agreements between the simulated and measured time series of free surface elevation at several wave gauges along the flume centerline, which further confirm the capabilities of the present model to simulate wave motion accurately. It should be noted that the wave amplitudes at 61.5 and 65.2 m are overpredicted by our model, which are also found in the simulation results in Kazolea et al. [36]. A possible explanation is that the flow in this area is very complex and the depth-averaged wave models may not be suitable for this situation.

Case 4: Runup of Solitary Wave on a Circular Island
To test the model's capability of simulating solitary wave runup over a 3D uneven bottom, laboratory experiments conducted by Briggs et al. [37] are used. The model setup and bottom geometry are shown in Figure 10  The model is tested under three experimental conditions, which are H/d = 0.045, 0.096, and 0.181, respectively, where the water depth d (= 0.32 m) is constant. The grid size for the computational domain is ∆x = ∆y = 0.05 m with the total grid number of 340,000. Since the model is constructed with smooth concrete, the Manning n value is set to zero. The time step is adaptive, and the total simulation time is 40 s. Figure 11 gives the comparisons of the water level fluctuations between the simulation and measurements. It can be found that overall the simulated water surface fluctuations agree well with the measured data. Figure 12 provides the maximum wave runup along the truncated cone, and shows the simulated values against the measurements. It is, therefore, reasonable to conclude based on the above results, that the model is also able to reproduce the measured wave propagation and runup over a 3D uneven bottom accurately.

Case 5: Breaking and Runup of Solitary Wave over a Three-Dimensional Complex Bathymetry
In this section, the simulated case aims to test the model's capability on simulating solitary wave breaking, runup and rundown over a three-dimensional complex bathymetry, which consists of a 1:30 slope, a triangular shaped shelf with a conical island located at the offshore point of the shelf as shown in Figures 13 and 14. Laboratory experiments were conducted in a 48.8 m long, 26.5 m wide and 2.1 m deep wave basin [38]. A piston-type wave maker generates a solitary wave from the left side in Figure 13, and the other three sides are wall. In the entire wave basin, there are nine wave gauges (G1-G9) used to measure free surface elevation and three acoustic Doppler Velocimeters (ADV1-ADV3) adopted to provide velocity data, as shown in Figure 14.    The water depth in the basin is kept constant at 0.78 m, and the incident wave height is 0.39 m. The computational domain includes 500 × 260 rectangular cells (∆x = ∆y = 0.1 m), and the Manning's coefficient is assigned as 0.014 s/m 1/3 . Figure 15 shows the water level fluctuations over the entire domain at different times. As seen in Figure 15, the solitary wave propagates from the left side to the right side. At t = 5 s, the wave encounters triangular-shaped shelf and begins to climb on the shelf. When t = 6.75 s, the wave passes over the conical island, the shoaling effect causes the wave to become asymmetric, and then the wave breaking occurs. After the wave breaking, the wave keeps climbing over the shelf, reaches the right boundary, and then reflects back at t = 29 s.  Figure 16 presents the simulated water level using the NHW and SWE model, respectively. For G1, G2, G4, and G5 stations, the wave amplitude calculated by the SWE model is lower than the experimental results. By considering the wave dispersion effects, the NHW model could simulate wave amplitudes more accurately and simulated results agree with the measurement fairly well. For G3 and G6 stations, there are certain deviations between the simulations and measurements, which also exist in other studies [19][20][21]. For G7, G8, and G9 stations (near the side wall), the simulations are fairly close to the measurements, which implies that the model can also accurately simulate the wave reflection from the side wall regions. Figure 17 performs the temporal velocity variations at the x direction for three velocity measurement stations using NHW and SWE model, respectively. Both the NHW and SWE model can capture most of the features of the velocities at ADV1-ADV3. In general, the developed NHW model exhibits reasonable agreement with the measurements of the solitary wave propagating over a three-dimensional complex bathymetry. The source code of the model, which is compiled with the Visual Studio 2008 and C++ language, executed on a single CPU core. The CPU time is about 2.015 h on a desktop computer with Intel i5-3470 CPU 3.2 GHz, 8 GB memory.

Conclusions
In this study, a two-dimensional depth-integrated non-hydrostatic wave model is developed. The model solves the governing equations with hydrostatic and non-hydrostatic pressure separately. To obtain the velocities under hydrostatic pressure conditions, the finite volume method with central-upwind scheme is adopted. For including the non-hydrostatic effect, the biconjugate gradient stabilized method is adopted to solve the non-hydrostatic pressure, and then the velocities can be corrected. The HFA method is used to deal with the wave breaking issue, and after the wave breaks, the non-hydrostatic model is transformed into the hydrostatic shallow water model, where the non-hydrostatic pressure and vertical velocity are set to zero.
Several analytical solutions and laboratory experiments are used to verify the accuracy and robustness of the developed model. In general, the numerical simulations are in good agreement with the theoretical or experimental results, which indicates that the model can be used to simulate large-scale wave motions in practical engineering applications. This model can only be used for the calculation of weakly nonlinear and dispersive waves, while for the simulation of highly nonlinear and dispersive waves, a multi-layer (three-dimensional) non-hydrostatic model is needed, which is our future work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Flux and Source Term Calculation
Appendix A. 1

. Nonnegative Reconstruction of Conserved Variables
The nonnegative reconstruction procedure proposed by Liang [30] is used to ensure the well-balanced property of the model, also in presence of wet/dry fronts. The conserved variables η, hu, hv and water depth h are reconstructed at four edges of cell (i, j). The reconstructed values at the cell interface (i + 1/2, j) are: where is a pre-defined threshold value of water depth and set to be 10 −8 m in all simulations. Moreover, the velocities are forced to be 0 when the cell is dry. The bed elevation can be calculated as: As suggested by Audusse et al.