A two-dimensional depth-integrated non-hydrostatic numerical model for nearshore wave propagation

Ocean Modelling 96 (2015) 187–202

Contents lists available at ScienceDirect

Ocean Modelling journal homepage: www.elsevier.com/locate/ocemod

A two-dimensional depth-integrated non-hydrostatic numerical model for nearshore wave propagation Xinhua Lu a,∗, Bingjiang Dong b, Bing Mao c, Xiaofeng Zhang a a

State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China Hydrology Bureau, Yangtze River Water Resource Commission, Wuhan 430010, China c Yangtze River Scientific Research Institute, Wuhan 430015, China b

a r t i c l e

i n f o

Article history: Received 28 January 2015 Revised 25 October 2015 Accepted 2 November 2015 Available online 10 November 2015 Keywords: Fractional step method Non-hydrostatic Pressure-correction method Shallow-water equations Shock-capturing Wave dispersivity Well-balanced

a b s t r a c t In this study, we develop a shallow-water depth-integrated non-hydrostatic numerical model (SNH model) using a hybrid finite-volume and finite-difference method. Numerical discretization is performed using the non-incremental pressure-correction method on a collocated grid. We demonstrate that an extension can easily be made from an existing finite-volume method and collocated-grid based hydrostatic shallow-water equations (SWE) model to a non-hydrostatic model. A series of benchmark tests are used to validate the proposed numerical model. Our results demonstrate that the proposed model is robust and well-balanced, and it captures the wet–dry fronts accurately. A comparison between the SNH and SWE models indicates the importance of considering the wave dispersion effect in simulations when the wave amplitude to water depth ratio is large.

1. Introduction Due to bed topography changes, currents, and interactions with coastal structures, waves that propagate from offshore regions to nearshore regions are subject to shoaling, refraction, diffraction, and breaking. The accurate prediction of these phenomena is of great importance for the study of nearshore hydrodynamics and is crucial for predicting sediment transport and bed deformation processes. Various approaches have been proposed for modeling these phenomena, e.g., by solving the Navier–Stokes equations (NSE), shallow-water equations (SWE), Boussinesq-type equations (BTE), and the SWE with non-hydrostatic pressure corrections (SNHE). All of these approaches have advantages and disadvantages. Advances in computational technology have facilitated the development of direct numerical simulation (DNS) models for solving the NSE, Reynolds-averaged version (RANS) models for solving the Reynolds-averaged NSE, and large-eddy simulation (LES) models for solving the filtered NSE. The flow details can be resolved using a DNS, RANS, or LES model, and the nonlinearity and wave dispersivity are fully maintained, which is useful for fundamental wave-related mechanistic studies (e.g., Lin et al., 1999; Hu et al., 2012). However, these models have disadvantages. First, they are computationally ex-

∗

Corresponding author. E-mail address: [email protected] (X. Lu).

http://dx.doi.org/10.1016/j.ocemod.2015.11.001 1463-5003/© 2015 Elsevier Ltd. All rights reserved.

© 2015 Elsevier Ltd. All rights reserved.

pensive and thus they are usually limited to small-scale simulations (especially with DNS and LES models). Second, these type of models generally require the implementation of complex algorithms to track the free surface. The implementation of a free surface tracking method requires additional computational time (the simulation domain also needs to be extended above the water surface) and, most importantly, there are also challenging issues related to the methods themselves, e.g., accurately conserving fluid mass and calculating spatial gradients near the free surface (Sussman and Puckett, 2000). Moreover, the irregular bed topography is non-trivial to process when developing a NSE based model. Compared with the 3D (or vertical 2D) DNS, RANS, or LES models, those for solving the depth-integrated equations are attractive due to their lower computational cost and because it is relatively easier to handle complex boundaries and they are simpler to code. The SWE can be solved effectively based on the theory of hyperbolic conservation laws. By using a shock-capturing numerical scheme, a SWE model without any tunable coefficients can predict the wave breaking phenomenon well; for instance, the wave energy dissipated in the wave breaking process and the maximum wave runup height after the wave breaks can be predicted reasonably well (Li and Raichlen, 2002; Tan and Chu, 2010; Hu et al., 2015). One of the disadvantages of the SWE models for applications in the nearshore area is the lack of representation of wave dispersivity. When short waves propagate over a long distance, significant errors can be predicted in the wave

188

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

phase and amplitude (Grilli et al., 2007; Kazolea and Delis, 2013). Unlike the SWE, the BTE can represent the weak dispersivity of waves propagating in shallow or intermediate depth water. Numerous studies have tried to improve the nonlinearity and dispersivity properties of the BTE. For instance, Gobbi and Kirby (1996) derived a fully nonlinear BTE model with high-order accuracy in terms of dispersion based on a specifically designed dependent variable (a weighted averaged velocity potential at two vertical locations). The BTE, especially the improved versions, contain high-order spatial or temporal derivatives (e.g., a fifth-order spatial derivative in the model proposed by Gobbi and Kirby (1996)), which is a challenging issue for numerical implementation. In addition, it is generally inappropriate to apply the BTE models directly to simulate wave breaking processes in the surf zone as a predefined breaking criterion is required, and it is usual to introduce an artificial viscosity term or a surface roller model to numerically dissipate the wave energy during the wave breaking process (Brocchini, 2013; Mccabe et al., 2013). To simultaneously represent the wave dispersivity satisfactorily and to avoid calculating high-order spatial derivatives, a SNH model for solving the SNHE was recently proposed by Stelling and Zijlema (2003), who decomposed the pressure into hydrostatic and nonhydrostatic parts, and by assuming appropriate vertical distributions for the non-hydrostatic pressure and vertical velocity, the SNHE was deduced by depth-integrating the RANS with suitable boundary conditions. Using the Arakawa C grid and employing a finite-difference scheme, Stelling and Zijlema (2003) obtained an elliptic equation for the non-hydrostatic pressure. The solution of this equation was then used to compute the horizontal velocity components at the new time level and finally, the water surface level was updated via the free surface kinematic boundary condition. Based on the study of Stelling and Zijlema (2003), various SNH models have been developed using either the non-incremental pressure-correction method (NPCM), incremental pressure-correction method (IPCM) (Guermond et al., 2006), or (fully) fractional step method (FSM) (Chorin, 1968). Note that the names of these solution methods may vary and they are elusive in previous studies, we here follow the definitions presented in Guermond et al. (2006), where the NPCM and IPCM are two versions of the pressure-correction method (Goda, 1979; van Kan, 1986). When solving the SNHE, both the NPCM and IPCM split the overall solution into two sub-steps, where the former solves the SWE (without considering the non-hydrostatic pressure effect) in the first substep and computes the non-hydrostatic pressure in the second substep, whereas the latter includes the non-hydrostatic pressure effect (approximated) to determine the intermediate velocity components in the first sub-step, as well as computing the pressure difference between the last and new time levels in the second sub-step. The FSM is very similar to the NPCM when solving the SNHE; these two can be distinguished that the NPCM applies a correction to the water depth (or surface elevation) after the non-hydrostatic pressure is solved in the second sub-step, whereas the FSM does not implement this procedure. In the following, we refer to the SNH models using the NPCM, IPCM, and FSM as SNH-NPC, SNH-IPC, and SNH-FSM, respectively. Zijlema and Stelling (2005) and Yamazaki et al. (2009) developed SNH-FSM models using finite-difference schemes with the same grid layout technique as that employed by Stelling and Zijlema (2003). Zijlema and Stelling (2005) also developed a SNH-IPC model and their numerical experiments showed that the SNH-IPC model maintained the wave shape well, while significant wave damping was predicted by their SNH-FSM model. Walters (2005) and Wei and Jia (2013) developed similar SNH models using finite-element schemes, where both employed a semi-implicit method to obtain provisional solutions in the first sub-step. The two models developed by Walters (2005) and Wei and Jia (2013) differed in that the former was a SNHFSM model whereas the latter was a SNH-NPC model. They also differed in their treatment of the convection terms as the former used a Lagrangian advection method to enhance the numerical stability,

whereas the latter adopted a bounded upwind scheme to avoid spurious oscillations. Cui et al. (2012) recently established a SNH-NPC model based on the SWE solver developed by Cui et al. (2010). In Cui et al. (2010), the finite-element algorithm for solving the SWE was implemented in the model called TsunAWI (Behrens, 2008) and it was converted into a finite-volume scheme with specifically defined finite volumes. This SWE model generally performs well, but minor non-physical oscillations in the solution can be predicted around discontinuities (Cui et al., 2010). In addition, the well-balanced property of this SWE model was not verified. Considering that most SNH models were developed based on SWE solvers using the finite-difference or finite-element schemes, while numerous SWE solvers have been developed based on the finite-volume method, and given that the state-of-the-art of the SWE solvers using the finite-volume method with schemes designed based on the hyperbolic conservation laws are robust and they can capture shocks and discontinuities well, in the present study, we aimed to develop a SNH model based on a SWE solver using the finitevolume method. In particular, we considered the design of a numerical scheme with the well-balanced property. We compared the performance of our proposed model with existing finite-difference or finite-element method-based SNH models. The remainder of this paper is organized as follows. In Section 2, we describe the governing equations and the numerical algorithms for the SNH model. In Section 3, numerical tests are used to verify various properties of the numerical model. Finally, conclusions and discussions are presented in Section 4. 2. Numerical model 2.1. Governing equations The 2D SNHE can be derived from the following 3D RANS

∂ u ∂v ∂ w + + = 0, ∂x ∂y ∂z

(1a)

  ∂τxy ∂τxz 1 ∂P 1 ∂τxx ∂ u ∂ uu ∂ uv ∂ uw + + + =− + + + , ∂t ∂x ∂y ∂z ρ ∂x ρ ∂x ∂y ∂z (1b)

  ∂τyy ∂τyz 1 ∂P 1 ∂τyx ∂v ∂ uv ∂vv ∂vw + + + =− + + + , ∂t ∂x ∂y ∂z ρ ∂y ρ ∂x ∂y ∂z (1c)

∂ w ∂ uw ∂vw ∂ ww + + + ∂t ∂x ∂y ∂z   ∂τzy ∂τzz 1 ∂P 1 ∂τzx =− + + + − g, ρ ∂z ρ ∂x ∂y ∂z

(1d)

where u, v, and w are the velocity components in the horizontal x-, y-, and vertical z-directions, respectively; t denotes time; ρ is the flow density; g is the acceleration due to gravity; and τxl xr (l, r = 1, 2, 3) denote the shear stresses, where x1 , x2 , and x3 denote the x-, y-, and z-coordinates, respectively. P denotes the pressure, which is decomposed into hydrostatic and non-hydrostatic components as

P = p + ,

(2)

where p is the hydrostatic pressure defined as ∂∂ pz = −ρ g and  is the

non-hydrostatic pressure. Following Stelling and Zijlema (2003), we assume that  and w vary linearly with water depth; hence, their depth-averaged values  + w +w are  = s 2 b and w = s 2 b . In the following, the ( · ) denotes the depth-averaging operator and the subscripts “s” and “b” denote the

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

values at the free surface and water bottom, respectively. For the nonhydrostatic pressure, it is usual to employ s = 0. Based on the linear assumptions for  and w, and using the decomposition given in Eq. (2) and the following kinematic boundary conditions at the free surface and water bottom

ws =

dη = dt

∂η ∂η ∂η + us + vs , ∂t ∂x ∂y

wb =

dzb ∂z ∂z = ub b + vb b , dt ∂x ∂y

where the vectors are defined by

⎡ ⎡ ⎤ η ⎢ ⎢ ⎥ U = ⎣ qx ⎦, F = ⎢ uqx + ⎣ qy

⎡

(3a)

(3b)

⎢ ⎢ ⎣

∂ qx ∂ uqx ∂vqx ∂η gn qx |q| ∂ b + + = −gh − − Kx b , −I ∂t ∂x ∂y ∂x ∂x h7/3

(4b)

∂ qy ∂ uqy ∂vqy ∂η gn2 qy |q| ∂ b + + = −gh − − Ky b , −I ∂t ∂x ∂y ∂y ∂y h7/3

(4c)

⎡

⎤

0

2 ⎢ ⎥ ∂ zb ⎥ ⎥ ⎢ − gn qx |q| ⎥ +⎢ ⎥. 7/3 h ∂x ⎥ ∂ zb ⎦ ⎣ gn2 q |q| ⎦ y −gη − ∂y h7/3

(8) Liang and Borthwick (2009) theoretically proved that Eq. (7) with vectors defined by Eq. (8) is exactly well-balanced over wet domains when the Godunov-type schemes are used. An explicit second-order strong stability-preserving total variation diminishing (TVD) Runge–Kutta scheme (Shu and Osher, 1988) for the discretization of Eq. (7) gives

( 1)

∂ w b = , ∂t ρh

(4d)

U i, j =

U ki, j

∂(η+z )

where I = 2hρ and Kxl = 21ρ ∂ x b (l = 1, 2); η and zb denote the wal ter surface and bed elevations above the datum; h = η − zb denotes the water depth; qx = u · h and qy = v · h are the discharges per unit



width in the horizontal directions; |q| = qx + qy and n is the Manning’s roughness coefficient. It should be noted that when deriving Eq. (4a)–(4d), the acceleration terms on the left-hand side of Eq. (1d) are omitted following Stelling and Zijlema (2003). The shear stresses are also ignored except that at the water bottom, the bed shear stress (friction) is closed with a quadratic expression defined by the second term on the right-hand side of Eq. (4b) and (4c). 2

In this section, we propose a SNH-NPC model based on a hybrid finite-volume and finite-difference method. The SWE model is used as the basis for developing the SNH-NPC model, so we first describe the numerical algorithms related to the SWE model. 2.2.1. SWE model The SWE model is developed based on the following SaintVenant equations (Eq. (4a)–(4c) without the non-hydrostatic pressure terms)

∂η ∂ qx ∂ qy + + = 0, ∂t ∂x ∂y

(5a)

∂ qx ∂ uqx ∂vqx ∂η gn2 qx |q| + + = −gh − , ∂t ∂x ∂y ∂x h7/3

(5b)

∂ qy ∂ uqy ∂vqy ∂η gn2 qy |q| + + = −gh − . ∂t ∂x ∂y ∂y h7/3

(5c)

Note that the overbar on the variables is dropped in the following for simplicity. Following Liang and Borthwick (2009), the surface gradient term −gh ∂∂η (l = 1, 2) on the right-hand side of Eq. (5b) and (5c) is reforx l

(6)

By virtue of Eq. (6), Eq. (5) can be recast into a system in conservation form as (7)

+ t

+

x

Gki, j−1/2 − Gki, j+1/2

y



+(

)

S0 ki, j

2 +

t

1) 1) F (i−1/2, − F (i+1/2, j j

+

x

2

) ) G(i,1j−1/2 − G(i,1j+1/2

, (9a)

y

 ( 1)

+ (S0 )i, j

( 2) k m Um i, j = U i, j + t S f (U i, j , U i, j ).

,

(9b)

(9c)

where the subscripts “i” and “j” denote the cell indexes; the superscripts “k”, “(1)”, “(2)”, and “m” denote the values at the last time step, at the time when the first and the second Runge–Kutta stage and the whole hydrostatic step calculations are accomplished, respectively. The friction source term vector in Eq. (9c) is discretized using a point implicit method (Fiedler and Ramirez, 2000). To evaluate the numerical fluxes in Eq. (9), the HLLC approximate Riemann solver (Toro et al., 1994) is invoked; for instance, in the xdirection, we have

F i−1/2, j =

⎧ L F i−1/2, j ⎪ ⎪ ⎨F *L

L if 0  ξi−1/2, , j

L M ξi−1/2,  0  ξi−1/2, , j j M R if ξi−1/2, j  0  ξi−1/2, , j R if ξi−1/2,  0, j

if

i−1/2, j F *R i−1/2, j F Ri−1/2, j

⎪ ⎪ ⎩

(10)

where F L = F (U L ) and F R = F (U R ). The superscripts “L”, “M”, and “R” denote that the variables are related to the left, middle (contact), and right waves at the cell interfaces, respectively. ξ L , ξ M , and ξ R denote the wave speeds, where considering the wetting and drying processes, they may be estimated as (Fraccarollo and Toro, 1995)



ξ = L



uR −2 ghR min uL −

 ξR =

mulated as

 2  g ∂ η − 2zb η ∂η ∂z =− − gη b . ∂ xl 2 ∂ xl ∂ xl

F ki−1/2, j − F ki+1/2, j

U ki, j + U (i,1j)

U (i,2j) =

2;

2.2. Numerical scheme

∂U ∂ F ∂ G + + = S, ∂t ∂x ∂y

⎤

qy

(4a) 2

−gh

⎡

⎥ ⎥ ⎢ g vq x ⎥, ⎢ η(η − 2zb )⎥ ⎦ ⎦, G = ⎣ 2 g vqy + η(η − 2zb ) uqy 2

⎤

0

⎤

qx

S = S0 + S f = ⎢−gη

the SNHE is obtained as follows by depth-integrating Eq. (1a)–(1d),

∂η ∂ qx ∂ qy + + = 0, ∂t ∂x ∂y

189



ghL , u∗ −



gh∗

L  if h = 0,



uL +  2 ghL max u + R



ghR , u∗

+



gh∗

if hL > 0,

R  if h = 0,

if hR > 0,

    ξ L hR uR − ξ R − ξ R hL uL − ξ L ξ = , hR (uR − ξ R ) − hL (uL − ξ L ) M

(11a)

(11b)

(11c)

where u∗ and h∗ , respectively, are calculated by

u∗ =

   1 L u + uR + ghL − ghR , 2

(12)

190

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

and

h∗ =

 

1 1 g 2

ghL +





ghR +

2

 1 L u − uR 4

.

The fluxes in the middle region are estimated as

F ∗ = [F1∗ F2∗ F3∗ ]T =



L

L

R

L

R

R

L

(19c)

(14)

 (qy )k+1 − (qy )m ∂ bk+1  i, j i, j k+1 − (Ky )m = −Ii,mj i, j (b )i, j , t ∂ y i, j



ξ F −ξ F +ξ ξ U −U , ξR − ξL R

(19b)

(13)

 (qx )k+1 − (qx )m ∂ bk+1  i, j i, j m k+1 = −Ii, j − (Kx )m i, j (b )i, j , t ∂ x i, j

and F *L and F *R in Eq. (10) are calculated from

⎡

F *L = ⎣

F1∗ F2∗

⎤

⎡

⎦, F *R = ⎣

L ∗ F1

F1∗ F2∗

wk+1 − wki, j i, j

⎤

⎦.

t (15)

R ∗ F1

v

v

Note that the superscript “T ” in Eq. (14) denotes the matrix transpose. The numerical flux in the y-direction can be calculated in a similar manner to that in the x-direction (Eqs. (10)–(15)). To reconstruct the Riemann states (η, h, qx , qy ) that appear in Eqs. (11)–(15), we first apply the Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) piecewise linear reconstruction method to achieve second-order accuracy in space (van Leer, 1979), and the minmod slope limiter (Roe and Baines, 1982) is used to prevent numerical oscillations in the vicinity of large gradients. The reconstructed Riemann states are then modified using a hydrostaticreconstruction approach (Audusse et al., 2004) to ensure nonnegative water depth solutions. In the method of Audusse et al. (2004), a single-valued bed elevation at the cell interface is defined as follows, e.g., in the x-direction, L L R R (zb )i−1/2, j = max (ηi−1/2, j − hi−1/2, j , ηi−1/2, j − hi−1/2, j ).

R hRi−1/2, j = max (ηi−1/2, j − (zb )i−1/2, j , 0),

(17)

and then the face value of η is modified to ensure that η = h + zb . Numerical balance is achieved by solving the reformulated equation system (Eqs. (7) and (8)), and thus no special treatment is needed to ∂z process the source term −gη ∂ xb (l = 1, 2) in Eq. (8); for instance, in the x-direction,

l

 R L ηi−1/2, + ηi+1/2, ∂ zb  j j (zb )i+1/2, j − (zb )i−1/2, j . −gη = −g ∂x  2 x

where qxl = qk+1 − qm xl xl (l = 1, 2). In particular, we need to consider the discretizations of Eq. (19b) ∂ and (19c). Simple central difference approximations for ∂ xb (l = 1, 2) l

in Eq. (19b) and (19c) would decouple the velocity (discharge per unit width) and the non-hydrostatic pressure, which may lead to a prediction with an incorrect checkerboard pressure field (Yu et al., 2002). Thus, to derive the pressure elliptic equation with an enhanced velocity–pressure coupling, we first rewrite Eq. (19b) and (19c) in the following form, so the discretizations are performed at the cell interfaces rather than at the cell centers (as expressed in Eq. (19b) and (19c)); for instance, at the cell interfaces (i − 1/2, j) and (i, j − 1/2), respectively, we have,

(qx )k+1 − (qx )m (b )k+1 − (b )k+1 i−1/2, j i−1/2, j i, j i−1, j m = −Ii−1/2, j t x k+1 −(Kx )m i−1/2, j (b )i−1/2, j ,

(20a)

(qy )k+1 − (qy )m (b )k+1 − (b )k+1 i, j−1/2 i, j−1/2 i, j i, j−1 = −Ii,mj−1/2 t y k+1 −(Ky )m i, j−1/2 (b )i, j−1/2 ,

(20b)

and at the cell interfaces (i + 1/2, j) and (i, j + 1/2), the following holds

(qx )k+1 − (qx )m (b )k+1 − (b )k+1 i+1/2, j i+1/2, j i+1, j i, j m = −Ii+1/2, j t x k+1 −(Kx )m i+1/2, j (b )i+1/2, j ,

(21a)

(qy )k+1 − (qy )m (b )k+1 − (b )k+1 i, j+1/2 i, j+1/2 i, j+1 i, j = −Ii,mj+1/2 t y k+1 −(Ky )m i, j+1/2 (b )i, j+1/2 .

(21b)

(18)

i, j

It should be noted that the SWE solver presented above predicts no spurious flow under stationary flow conditions over the wet domain. However, in the dry area or when a wet cell is adjacent to a dry cell, non-physical flux can still be predicted. In the present study, cells for which the depth is less than 1 × 10−10 m are denoted as dry cells and their discharge per unit width is set to zero. When the bed elevation of a dry cell is higher than the water surface elevation in its adjacent cell, the flux at the relevant cell interface is forced to be zero via a modification of the local bed slope (Brufau et al., 2002). In this manner, no spurious flow is generated over the dry domain or wet–dry adjacent area, and thus the SWE model described in this section is well-balanced for simulations involving wetting and drying processes over complex bed topography. 2.2.2. SNH-NPC model For a SNH-NPC model, the provisional solution Um obtained from the hydrostatic step with the SWE solver described in Section 2.2.1 is updated with the considered non-hydrostatic pressure terms by solving the following semi-discretized equations

  ηi,k+1 − ηi,mj ∂ qy  ∂ qx  j + + = 0, t ∂ x i, j ∂ y i, j

(19d)

(16)

Subsequently, the reconstructed water depth at the cell interface is redefined as L hLi−1/2, j = max (ηi−1/2, j − (zb )i−1/2, j , 0),

(b )k+1 i, j , ρ hm i, j

=

(19a)

By subtracting Eq. (20a) from Eq. (21a) and by assuming that

qk+1 x

≈ uk+1 hm ,

(22)

we obtain the velocity gradient in the x-direction as



k+1 − uk+1 uk+1 um − um ∂ u  i+1/2, j i−1/2, j i+1/2, j i−1/2, j = =  ∂ x i, j x x  m  t Ii+1/2, j  − − (b )k+1 (b )k+1 i+1, j i, j x2 hm i+1/2, j + − +

x(Kx )m  i+1/2, j 2hm i+1/2, j m Ii−1/2, j

hm i−1/2, j

+ (b )k+1 (b )k+1 i, j i+1, j



− (b )k+1 (b )k+1 i, j i−1, j



  k+1 + (  ) . (b )k+1 b i, j i−1, j

x(Kx )m  i−1/2, j 2hm i−1/2, j



(23)

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

The velocity gradient in the y-direction, ∂v ∂ y , can be discretized in a similar manner, and that in the z-direction is calculated by discretizing Eq. (19d) as

(ws )k+1 + (wb )k+1 − (ws )ki, j − (wb )ki, j (b )k+1 i, j i, j i, j . = 2t ρ hm i, j

(24)

Note that the error due to the assumption made in Eq. (22) is rather small because the change in η (and thus in h) is negligible in the non-hydrostatic step calculations (as discussed later in Appendix A). If we discretize Eq. (1a) implicitly as follows

uk+1 − uk+1 i+1/2, j i−1/2, j

x

+

vk+1 − vk+1 (ws )k+1 − (wb )k+1 i, j+1/2 i, j−1/2 i, j i, j + = 0, y hm i, j (25)

then by substituting the discretized spatial velocity gradients (Eqs. (23) and (24)), which are functions of the non-hydrostatic pressure at the new time level, into Eq. (25) and after some simple mathematical manipulations, we obtain the following elliptic equation for the non-hydrostatic pressure at the new time level

Ai, j (b )k+1 + Bi, j (b )k+1 + Ci, j (b )k+1 + Di, j (b )k+1 i−1, j i, j i+1, j i, j−1 + Ei, j (b )k+1 = Si, j , i, j+1

(26)

where the coefficients are given by

t t (Kx )i−1/2, j + , 2x hm 2 ρ x 2 i−1/2, j m

Ai, j = − Bi, j =

2t t t + + 2 2 ρ  x ρ y 2 ρ(hm ) i, j   m m t (Kx )i+1/2, j (Kx )i−1/2, j − − 2x hm hm i+1/2, j i−1/2, j   m m ( K ) ( K y i, j+1/2 y )i, j−1/2 t − − , 2 y hm hm i, j+1/2 i, j−1/2

t t (Kx )i+1/2, j − , 2x hm 2 ρ x 2 i+1/2, j

Ei, j = − Si, j = − +

vm − vm i, j+1/2 i, j−1/2 x y − (ws )ki, j − (wb )ki, j 2(wb )k+1 i, j hm i, j

−

each mesh point. Finally, from Eq. (19a), we obtain ηk+1 ≡ ηm , which means that the water surface level remains steady at the new time = wk+1 =0 level. Note that from Eqs. (24) and (28), we obtain wk+1 s b all over the domain, and thus the source term S remains zero and the solution of Eq. (26) does not change in the next time step calculations. This analysis indicates that no spurious flow is induced through the non-hydrostatic step calculations. The hydrostatic solver presented in Section 2.2.1 is also well-balanced, so the proposed SNH-NPC model is numerically well-balanced.

2.2.3. Boundary conditions and stability criteria To process the boundary conditions, ghost cells are created outside the simulation domain. In this manner, the numerical fluxes at the boundaries can be calculated in the same way as those at the interfaces between interior cells. For all of the numerical tests considered in this study, an open boundary condition is applied at all of the domain boundaries. This means, all of the physical variables (except b at the inlet boundary) at the ghost cells are zeroth-order extrapolated from those at the interior cells. In this study, the b at the ghost cells at the inlet boundary is set to zero; for waves of a particular type to be generated at the inlet boundary, b can be prescribed at the relevant ghost cells based on some theoretical analysis (e.g., see Cui et al. (2012)). Since explicit schemes are used for the SWE model, the Courant– Friedrichs–Lewy condition should be satisfied to ensure the numerical stability

⎞ x 

|uki, j | +

, ghki, j

y 

|vki, j | +

⎠,

(29)

ghki, j

.

3. Numerical tests (27)

Note that in Eqs. (26) and (27), Kxl (l = 1, 2) at the cell interface are calculated using a central difference scheme, where the water depth at the cell interface is linearly interpolated from those at the cell centers. According to previous studies (Walters, 2005; Wei and Jia, 2013), wb at the new time level in Eq. (27) can be estimated from Eq. (3b), which in a discretized form is expressed as

(wb )k+1 ≈ um i, j i, j

Then, from Eqs. (19b) and (19c), we have qk+1 = qxl = 0 (l = 1, 2) at xl

where the coefficient Cr should be specified in the range of 0 < Cr ≤ 1.

m t t (Ky )i, j+1/2 − , 2y hm 2 ρ y 2 i, j+1/2

um − um i+1/2, j i−1/2, j

∂ qx

of the term ∂ x l (l = 1, 2) in Eq. (19a), our numerical experience has l shown that the difference in the results obtained using the first-order upwind and the second-order central difference schemes is negligible, so a second-order central difference scheme is employed in the present study. Under an initially quiescent flow condition, the source term S in Eq. (26) is equal to zero at each mesh point and the unique solution for Eq. (26) is thus bk+1 ≡ 0 all over the simulation domain.

t = Cr × min ⎝

t t (Ky )i, j−1/2 =− + , 2y hm 2 ρ y 2 i, j−1/2 m

Di, j

and ηk+1 is then updated by Eq. (19a). Note that for the discretization

⎛

m

Ci, j = −

191

(zb )i+1, j − (zb )i−1, j (zb )i, j+1 − (zb )i, j−1 + vm . i, j 2 x 2 y

In this section, 1D and 2D numerical tests are used to verify the performance of the proposed numerical model. For the 1D simulations, Ny = 1 is used, where Nxl (l = 1, 2) denote the mesh-cell numbers in the xl directions. In this section, unless stated otherwise, the SWE and the SNH-NPC models refer to those developed in Sections 2.2.1 and 2.2, respectively. In all of the simulations, g = 9.81 m/s2 and ρ = 1000 kg/m3 .

3.1. 1D solitary waves propagating over a constant water depth channel (28)

To solve the elliptic equation expressed in Eq. (26), we employ the biconjugate gradient stabilized method (van de Vorst, 1992) and the incomplete LU preconditioner (Saad, 2003) is utilized to accelerate k+1 the convergence of the iteration. After obtaining bk+1 , qk+1 x , qy ,

are computed from Eqs. (19b), (19c), and (24), respectively, and wk+1 s

In this section, we simulate 1D solitary waves propagating over a 1000 m-long frictionless channel. Similar tests were used by Stelling and Zijlema (2003) and Ma et al. (2012) to evaluate the wave dispersivity representation and mass conservation property of their SNH or RANS models. The initial water surface elevation and the velocity components are prescribed according to the analytic solution derived

192

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

11

10.8

η (m)

Table 1 Wave height computed by different models with Nx = 600 for a solitary wave propagating over a constant water depth channel test with H/d = 0.2. The SNH-W05 and SNH-SZ03 model results are extracted from their publications.

Nx=125 Nx=250 Nx=500 Nx=1000 Nx=2000 Theory

10.6

10.4

10.2

10 600

650

700

750

800

850

x (m) Fig. 1. Water surface profiles computed by the SNH-NPC model using different meshcell numbers at t = 50 s for a solitary wave propagating over a constant water depth channel test with H/d = 0.1.

by Keulegan and Patterson (1940) based on the potential theory:

!

η(x, t ) = d+Hsech

2

u(x, t ) =

"

3H (x − ct − X0 ) , 4d 3

c [η(x, t ) − d] , h(x, t )

ws (x, t ) = −η(x, t )

∂ u(x, t ) , ∂x

(30a)

(30b)

(30c)

where d denotes the still water depth; H is the initial wave height of the solitary wave; c denotes the wave phase velocity; and X0 is the initial center of the solitary wave, which is taken as 200 m in the simulation. Without considering the bed friction, theoretically  the solitary wave should travel with a constant phase speed c = g(H + d) and

Time

SNH-NPC (m)

SNH-W05 (m)

SNH-SZ03 (m)

t ≈ 20 s t ≈ 40 s

1.966 1.952

1.888 1.836

1.900 1.888

maintain its shape, where all of the physical variables are unchanged according to Eq. (30). We consider two test cases with different wave height to water depth ratios, H/d = 0.1 and H/d = 0.2. For both test cases, we use d = 10 m and Cr = 0.5. The spatial resolution is selected based on a grid-independence study. Simulations using the SNH-NPC model on different meshes with Nx = 125, 250, 500, 1000, and 2000 are implemented and the computed water surface profiles at t = 50 s are presented in Fig. 1. It can be seen that the results obtained with Nx = 1000 and Nx = 2000 are close to each other, and a mesh with Nx = 2000 is used to ensure the numerical accuracy. Figs. 2 and 3 show the numerical solutions computed by the SWE and SNH-NPC models with Nx = 2000 for the two test cases. The predicted water surface level, horizontal velocity, and vertical velocity at the free surface at t = 0 s, 25 s, and 50 s are compared with the analytic solutions. It can be seen that the shape of the solitary wave and the velocity components are maintained accurately by the SNH-NPC model for both test cases. However, without considering the wave dispersion effect, the shape of the solitary wave gradually skews and steepens at the wave front, as predicted by the SWE model shown in Fig. 2(a), which is even worse for the higher H/d ratio test case, as shown in Fig. 3(a). At t = 50 s, the wave height predicted by the SWE model decreases to about 70% of its theoretical value for the H/d = 0.2 test case. In Table 1, we list the wave height at two time instances predicted by the SNH-NPC model proposed in the present study, and those given by Walters (2005) using a finite-element scheme-based SNH-FSM model and by Stelling and Zijlema (2003) using a one-step velocity–pressure coupling and finite-difference scheme-based SNH model, where the latter two models are abbreviated as SNH-W05

Fig. 2. Comparison of the water surface level (a), horizontal velocity (b), and vertical velocity at the free surface (c) at t = 0 s, t = 25 s, and t = 50 s computed using the SNH-NPC model (——) and SWE model (– – –) for a solitary wave propagating over a constant water depth channel test with H/d = 0.1. The analytic solutions are given in symbols ().

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

193

Fig. 3. Comparison of the water surface level (a), horizontal velocity (b), and vertical velocity at the free surface (c) at t = 0 s, t = 25 s, and t = 50 s computed using the SNH-NPC model (——) and SWE model (– – –) for a solitary wave propagating over a constant water depth channel test with H/d = 0.2. The analytic solutions are given in symbols ().

and SNH-SZ03, respectively, in Table 1. For all of the results listed in Table 1, the same grid resolution (Nx = 600) is used. Note that the results for the SNH-W05 model are available and given at t ≈ 18.4 s and t ≈ 36.9 s, whereas for the other two models, the results are given at t = 20 s and t = 40 s. It can be seen that the wave height predicted by the SNH-W05 model decreases by more than 0.11 m during the first 20 s of the simulation, where at t ≈ 36.9 s, compared with the initial wave profile, the wave height drops by about 0.16 m (which should be larger at t = 40 s). Among the three models, our proposed SNH-NPC model is the most accurate at maintaining the wave shape, where its wave height reduction is less than half of that given by the SNH-SZ03 model (e.g., a 0.05 m reduction is predicted by the proposed model versus a 0.11 m reduction by the latter model at t = 40 s). We mention that Cui et al. (2012) also obtained accurate results by using their SNH-NPC model for this test, but the spatial resolution employed in their simulation was not reported and thus a direct comparison cannot be made.

Fig. 4. Sketch of a solitary wave propagating over a sloped beach onto a shelf test. The gray shaded area shows the bed topography. The blue shaded area with crossed lines shows the sponge layer area used in the numerical simulation. The dash-dotted line (– · – · –) marks the initial center of the solitary wave. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

3.2. A 1D solitary wave propagating over a beach onto a shelf A solitary wave propagating over a beach onto a shelf was studied experimentally and analytically by Madsen and Mei (1969). The experimental setup is shown in Fig. 4. In the experiment, the still water depth changed from d = 0.0762 m at the left to d1 = d/2 at the right of the flume, with a 1: 20 sloped beach in between, ranging from x = 2 m to x = 2.76 m. The initial wave height to water depth ratio was H/d = 0.12. Wave gauges were used in the experiment to record the water surface elevations at locations A (x = 1.59 m), B (x = 2.76 m), C (x = 3.65 m), and D (x = 4.51 m). This test was used by Guo et al. (2013) to validate the accuracy of their SNH-NPC model. In the simulation, the solitary wave is initially centered at x = X0 = −0.8 m, and the water level and the velocity components are prescribed according to Eq. (30). At the eastern boundary, a sponge layer in the region x ∈ [6.0 m, 9.0 m] is used to minimize the wave reflection. A uniform grid is used with x = 0.01 m and a dynamic time step is employed with Cr = 0.5. Fig. 5 shows a sequence of time snapshots for the water surface elevation computed by the SWE and SNH-NPC models. Note that the gauge locations in the experiment are marked in Fig. 5 by the vertical dash-dotted lines. To facilitate the comparison, a relative time t − tr is used, where tr denotes the time when the crest

of the solitary wave passes over Gauge A. As shown in Fig. 5(a), the shape of the solitary wave is well-preserved by the SNH-NPC model when the wave propagates in the offshore constant water depth area (x < 2.0 m). Without considering the wave dispersion effect, the SWE model predicts a skewed water surface profile, as already found in the previous test in Section 3.1. When the wave propagates over the sloped beach, as shown in Fig. 5(b), the wave front steepens and the wave amplitude increases slightly, as predicted by the SNH-NPC model. At the instance t − tr = 2.36 s (Fig. 5(c)), a second wave crest is observed in the water surface profile (indicated by the arrow in Fig. 5(c)), and it finally develops to form a small-amplitude solitary wave, as shown in Fig. 5(d). This phenomenon of the disintegration of a solitary wave passing a sloped beach or a vertical step onto a shallow shelf was also reported previously through experimental studies (Madsen and Mei, 1969; Street et al., 1968) and numerical investigations based on the 3D RANS models (e.g., Liu and Cheng (2001)). The SWE model cannot reproduce the solitary wave disintegration process. As shown in Fig. 5(b)–(d), the SWE model predicts a bore-like wave shape with a decreasing wave amplitude when the wave propagates onto the shelf, but the experiment data show that the peak wave amplitude increases to exceed its initial value at later times, which is predicted well by the SNH-NPC model.

194

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

Fig. 5. Computed water surface elevation at (a) t − tr = 0 s, (b) t − tr = 1.2 s, (c) t − tr = 2.36 s, and (d) t − tr = 3.4 s. The solid line (——) and dashed line (– – –) represent the numerical solutions obtained by the SNH-NPC and SWE models, respectively. The vertical dash-dotted lines (– · – · –) denote the gauge locations in the experiment reported by Madsen and Mei (1969).

Fig. 6. Time history of the computed water surface elevation at gauge locations A (a), B (b), C (c), and D (d). The solid (——) and dashed (– – –) lines represent the solutions obtained by the SNH-NPC and SWE models, respectively. The dash-dotted line (– · – · –) represents the numerical solution given by Guo et al. (2013). The symbols () show the experimental data reported by Madsen and Mei (1969).

Fig. 6 compares the time series of the computed water surface level by the SWE and SNH-NPC models with the data measured by Madsen and Mei (1969). We also present the numerical results computed by Guo et al. (2013) using a finite-difference scheme-based SNH-NPC model with the same spatial resolution as that used in our simulation. The model proposed by Guo et al. (2013) is abbreviated as SNH-GKJ13 hereafter. It should be noted that after carefully checking, we found that the numerical result presented by Guo et al. (2013) at Gauge D (Fig. 6(d)) might be erroneous in their work since the predicted wave phase at this wave gauge location differs significantly (but not at the other gauges) from the measured data as well as our numerical results (one should also note that the times of the experimental data for Gauges B and C might also be incorrectly presented in their work). Thus, in Fig. 6(d), we present the corrected numerical prediction obtained by the SNH-GKJ13 model at Gauge D after adjusting the wave crest passing time to match our simulation results. As shown in Fig. 6, when the SWE model is used, a substantial difference is found in the wave phase compared with the measured data, where the wave arrives earlier than that in the experiment at all four wave gauge locations. The SWE model also predicts a much smaller wave amplitude than the measured value at later times when the wave travels onto the shelf. A comparison of the SNHNPC and the SNH-GKJ13 models shows that the two models perform similarly, and the predicted water surface elevations are generally in phase with the measured data. Slight differences between simulated and measured wave phases at Gauges C and D are predicted by both the SNH-NPC and SNH-GKJ13 models (see Fig. 6(c) and (d)) when the

Table 2 Values of Eη for different models at different wave gauge locations for a solitary wave propagating over a sloped beach onto a shelf. Location

SWE ( × 10−2 )

SNH-NPC ( × 10−2 )

SNH-GKJ13 ( × 10−2 )

A B C D

2.21 2.63 2.54 2.25

0.73 0.58 1.71 2.15

0.93 0.85 2.05 2.63

disintegration of the wave occurs and the small-amplitude wave forms (t − tr  2.36 s, cf. Fig. 5(c) and (d)). To quantitatively assess the accuracy of the numerical models, the following global error estimator is used (Lu et al., 2012) Nm 1 # Eφ = Nm

im =1

   φp − φ m     φm  ,

(31)

im

where φ can be any of the physical variables; the subscripts “p” and “m” denote the predicted result and the measured data, respectively; and Nm denotes the number of measured data. In general, a small value of Eφ indicates a better performance of the numerical model in terms of numerical accuracy. The values of Eη for the SWE, SNH-NPC, and SNH-GKJ13 models at different wave gauge locations are listed in Table 2, which show that the performance of the SNH-NPC model is similar to the SNH-GKJ13 model, but remains slightly better.

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

195

Fig. 7. Computed water surface profiles at t˜ = 30 (a), t˜ = 40 (b), t˜ = 45 (c), t˜ = 50 (d), t˜ = 60 (e), and t˜ = 70 (f). The thick solid line (——) and the dashed line (– – –) represent the numerical predictions obtained by the SNH-NPC and SWE models, respectively. The dash-dotted lines show the numerical results given by Titov and Synolakis (1995). The symbols () are the experimental data reported by Synolakis (1986). The thin solid line represents the bed topography. For each sub-figure, a magnified view in the lower panel is given for the region of x/d ∈ [−8, 2]. .

3.3. A 1D nonbreaking wave running up a plane beach In this test, we simulate a solitary wave running up and down a plane beach. This test was studied experimentally by Synolakis (1986). The simulation domain comprises a horizontal bottom region with bed elevation zb = 0 m, followed by a plane beach with slope 1: 19.85. The origin (x = 0 m) is defined at a location where the line z = d intersects the beach surface. In the simulation, the initial conditions are prescribed according to Eq. (30). Initially,  the solitary wave is centered at X0 /d = −19.85 −

4d 3H arccosh

1 0.05

with H/d = 0.0185. The

simulation domain (x ∈ [−100 d, 15 d]) is discretized with a uniform grid of x = 0.125 d. The time step is constrained with Cr = 0.7. Fig. 7 illustrates the water surface profiles computed by the SWE and  SNH-NPC models at various normalized time instances (t˜ = t gd). We also show the numerical predictions obtained by Titov and Synolakis (1995) using a finite-difference scheme-based SWE model

(SWE-TS95 in the following). To facilitate the comparison, magnified views are shown around the toe of the beach at the bottom of each sub-figure in Fig. 7. As shown in Fig. 7, the results computed by the SWE and SNH-NPC models are indistinguishable, and they all agree well with the measured data. The tiny difference between the numerical predictions obtained by the SWE and SNH-NPC models indicates that the wave dispersion effect is insignificant in this test. It is obvious that compared with the observed data, the SWE and the SNH-NPC models developed in this study yield more favorable results than the SWE-TS95 model, in particular the former two models are superior to the latter at capturing the wet–dry fronts. 3.4. A 1D breaking solitary wave running up a plane beach A breaking solitary wave propagating over a 1: 19.85 sloped beach is studied to further assess the performance of the proposed numerical model. All of the parameters are taken the same as in Section 3.3,

196

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

Fig. 8. Snapshots of the water surface profile at t˜ = 10 (a), t˜ = 15 (b), t˜ = 20 (c), t˜ = 25 (d), t˜ = 30 (e), t˜ = 35 (f), t˜ = 40 (g), t˜ = 45 (h), t˜ = 50 (i), t˜ = 55 (j), t˜ = 60 (k), and t˜ = 65 (l). The thick solid line (——) and the dashed line (– – –) represent the numerical predictions obtained by the SNH-NPC and SWE models, respectively. The dash-dotted lines (– · – · –) show the numerical results reported by Zelt (1991). The symbols () are the experimental data reported by Synolakis (1986). The thin solid line represents the bed topography.

except that H/d = 0.28 is considered. This test was used by Zelt (1991) and Li and Raichlen (2002) to evaluate their BTE and SWE models, respectively. Fig. 8 presents the water surface profiles computed by the SWE and SNH-NPC models. We also show the numerical results reported by Zelt (1991) via solving the BTE using a Lagrangian finite-element method, where this model is abbreviated as BTE-Z91 in the following. The numerical results are compared with the measured data reported by Synolakis (1986). Note that the numerical results for the BTE-Z91 model at t˜ = 35 and t˜ = 40 are unavailable, so the relevant results are not plotted in Fig. 8. As shown in Fig. 8(a) and (b), as the wave passes the toe of the beach (x/d = −19.85), the wave profile loses symmetry and the wave front steepens due to the shoaling effect. The wave amplitude predicted at this stage by the SWE model is lower than that observed in the experiment. In fact, at this stage, the SWE model already simulates a broken wave profile, whereas the wave broke between t˜ = 20

and t˜ = 25 in the experiment, so the wave phase is simulated incorrectly by this model at this stage. By considering the wave dispersion effect, the SNH-NPC model gives an accurate prediction of the wave phase speed at this stage and the simulated wave amplitude matches with the measured data fairly well. At the instance t˜ = 20 (Fig. 8(c)), the wave amplitude becomes high. The SWE model still predicts a faster wave phase speed than the measured data. The results obtained by the SNH-NPC and BTE-Z91 models are in phase with the observed data, where the wave phase is predicted more accurately by the former model. The simulated wave amplitude is slightly higher for the SNH-NPC model and lower for the BTE-Z91 model with respect to the measurements. The broken wave climbs the sloped beach, as shown in Fig. 8(d)– (f). All three numerical models (SWE, SNH-NPC, and BTE-Z91) give relatively accurate solutions at this stage, although BTE-Z91 is not as accurate as the other two models in describing the wetting and drying processes (see Fig. 8(d)). The wave then retreats and the formation

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

197

Table 3 Values of Eη for different models at different time instances for the breaking wave propagating over a plane beach. Time

SWE (×10−2 )

SNH-NPC (×10−2 )

BTE-Zelt (×10−2 )

t˜ = 10 t˜ = 15 t˜ = 20 t˜ = 25 t˜ = 30 t˜ = 45 t˜ = 50 t˜ = 55 t˜ = 60 t˜ = 65

1.85 2.34 2.29 0.92 0.62 1.06 1.66 1.56 1.60 1.09

0.74 1.19 1.70 0.97 0.79 1.16 1.61 1.38 1.65 1.07

1.00 1.41 1.52 1.98 1.16 1.41 1.70 1.60 2.57 2.32

of the backwash bore can be seen clearly (Fig. 8(g)–(l)). The numerical results computed by the SWE and the SNH-NPC models match reasonably well with the experimental data during the wave retreating process. A comparison of the experimental data for water surface elevations between the instances t˜ = 55 and t˜ = 60 (Fig. 8(j) and (k)) show that the backwash bore seems to travel onshore (rather than offshore) during this period in the experiment, which implies a reverse wave breaking. This phenomenon was reported by Dalrymple and Rogers (2006) in an investigation of wave breaking on a 1: 13.5 plane beach, and it might also have occurred during this period in the experiment by Synolakis (1986). Various numerical schemes are unable to simulate the reverse wave breaking process well. For instance, even when sufficient artificial viscosity is introduced to enhance the numerical stability, an oscillation in the solution is still predicted by the BTE-Z91 model (see the magnified views in Fig. 8(i)–(l)). In addition, the wave phase may be predicted inaccurately, as shown by the numerical results reported by Li and Raichlen (2002) via solving the SWE using a weighted essentially non-oscillatory scheme. The values of Eη calculated for the SWE, SNH-NPC, and BTE-Z91 models are listed in Table 3, which demonstrates that the SNH-NPC model performs much better than the SWE model when the wave dispersivity and nonlinearity are both important (t˜  20). By contrast, in the wave shoaling and retreating processes after the wave breaks (t˜  25), the wave dispersion effect is of no importance. If we compare the results obtained by the SNH-NPC and BTE-Z91 models, it is seen that except for t˜ = 20, the former model is generally more accurate than the latter. It should be noted that Zelt (1991) used an empirical artificial viscous term to consider energy dissipation in the wave breaking process and the value of the empirical coefficient was varied by a factor of eight to select the value that agreed best with the laboratory data. However, no tunable coefficient is required in the present SNH-NPC model. Therefore, our proposed SNH-NPC model might be more convenient and suitable for describing surf zone dynamics. 3.5. 2D solitary waves running up a circular island We extend the validation of the numerical model to twodimensional problems by considering solitary waves running up a circular island. To study this problem, Briggs et al. (1995) conducted experiments in a large wave basin (25 m-long by 30 m-wide by 0.6 m-deep). A schematic view of this experiment is shown in Fig. 9. The island with a height of 0.625 m is centered at the point (12.96 m, 13.80 m) and its shape is a truncated circular cone with a diameter of 7.2 m at the base and 2.2 m at the crest. In total, 27 capacitance wave gauges (some are shown in Fig. 9) were used to measure the water surface elevations. The laboratory experiments were performed with different still offshore water depths (d = 0.32 m and 0.42 m) and relative wave heights (H/d = 0.045, 0.096 and 0.181). In this test, we only consider the cases with the smaller still water depth (d = 0.32 m) and the higher water height to depth ratios (H/d = 0.096 and 0.181), which are more challenging for numerical simulations.

Fig. 9. Sketch of the wave basin and the gauge locations for the solitary wave running up a circular island, where (a) and (b) are the plan and side views, respectively.

In the simulations, the solitary waves are initially centered at x = 0 m and the initial conditions are prescribed according to Eq. (30). Note that the length of the simulation domain on the western boundary is extended by 5 m to minimize the reflection following Ai and Jin (2012). The simulation domain is discretized with mesh sizes of x = y = 0.05 m, and Cr = 0.25 is used. The roughness coefficient is set to n = 0.01 since the island and the basin in the experiment were built of concrete with a smooth bed surface. The data measured at Gauge 2 provides a reference for adjusting the timing of the computed waveforms. Fig. 10 shows the water surface profiles simulated by the SNHNPC model at various time instances for the H/d = 0.181 case. It can be seen that the solitary wave propagates from the left of the domain to reach and inundate the circular island (Fig. 10(a)). The wave then diffracts around the island to form two edge waves (Fig. 10(b)), which later collide with each other at the lee side of the island (Fig. 10(c)). After the collision, the edge waves propagate toward the eastern boundary and a very complex water surface profile remains (Fig. 10(d)). In Fig. 11, we plot the time series of the water surface elevation computed by the SNH-NPC model at different gauge locations for H/d = 0.096 (Fig. 11(a)–(e)) and H/d = 0.181 (Fig. 11(f)–(j)). The data measured by Briggs et al. (1995) are also shown. Previous studies of these two test cases generally require the calibration of the incoming wave amplitude to match the measured data, e.g., at Gauge 2 in Yamazaki et al. (2009) and at Gauge 4 in Cui et al. (2012). However, no calibration is needed in our current simulations. One can see from Fig. 11(a) and (f) that the computed wave amplitudes at Gauge 2 match fairly well with the observed data, which once again signifies that the proposed model can maintain the wave shape well in constant water depth channel simulations (also see Section 3.1). For both test cases, overall good agreement is provided between the computed water surface elevations and the measured data, and our numerical results also match well with the numerical predictions obtained by Ai and Jin (2012) using a 3D non-hydrostatic model. A small deviation is simulated in the wave phase around the wave peak at Gauge 22 between the numerical result and the measured data, which is also predicted by other SNH models (e.g., Yamazaki et al. (2009); Cui et al. (2012); Wei and Jia (2014)). We remark that oscillations in solutions can easily be generated (especially at Gauge 22), as shown by the numerical results for the H/d = 0.181 test case by Yamazaki et al. (2009)

198

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

Fig. 10. Snapshots of the water surface profile computed for the solitary wave running up a circular island with H/d = 0.181 at t = 30.23 s (a), t = 32.28 s (b), t = 33.52 s (c), and t = 35.52 s (d). The color shows the distribution of the computed water depth. Note that different angles of view are used in top and bottom panels.

Fig. 11. Time series of the water surface profile at different gauge locations in the test case of the solitary wave running up a circular island for H/d = 0.096 (a)–(e) and H/d = 0.181 (f)–(j). The thick solid lines (——) and the symbols () show the numerical results computed by the SNH-NPC model and the measured data reported by Briggs et al. (1995), respectively.

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

(a)

199

(b)

Fig. 12. Water surface profile at t = 20 s (a) and the time history of the maximum absolute value of the velocity components in the simulation domain (b) computed by the SNH-NPC model for the well-balanced test. In (b), the velocity units are m/s.

(a) (b)

(b) (b)

11.2

SNH-NPC Nx=125 SNH-IPC Nx=125 SNH-NPC Nx=2000 SNH-IPC Nx=2000

11.2

11

η (m)

η (m)

11

10.8

10.8

10.6

10.4

SNH-FSM Nx=125 SNH-IPC Nx=125 SNH-FSM Nx=2000 SNH-IPC Nx=2000

10.6

700

720

x (m)

740

10.4

700

720

x (m)

740

Fig. 13. Comparison of the water surface profiles computed at t = 50 s for a solitary wave propagating over a constant water depth channel with H/d = 0.1: (a) SNH-NPC versus SNH-IPC; and (b) SNH-FSM versus SNH-IPC.

Table 4 Errors in the maximum wave amplitude and the corresponding occurrence time simulated by different models in the test case of the solitary wave running up a circular island for H/d = 0.181. Time

SNH-NPC

SNH-WJ14

SNH-YKC09

Hmax (m) tmax (s)

0.014 0.156

−0.013 0.340

0.029 0.240

elevation η = d is prescribed all over the domain. Fig. 12(a) and (b) shows the water surface profile computed at t = 20 s and the maximum absolute value of the computed velocity components in the simulation domain, respectively. It is clear that the steady stationary flow state is preserved, as expected (see also the analysis given in Section 2.2.2).

4. Conclusion and discussion

and for both test cases by Wei and Jia (2014); however, there are no oscillations in solutions of the proposed SNH-NPC model. Table 4 lists the errors in the predicted maximum wave amplitude and the corresponding occurrence time at Gauge 22 for various SNH models. The models considered include the proposed SNHNPC model, the SNH-NPC model developed by Wei and Jia (2014) using a finite-element scheme, and the SNH-FSM model developed by Yamazaki et al. (2009) using a finite-difference scheme, where the latter two models are abbreviated as SNH-WJ14 and SNH-YKC09 in Table 4, respectively. In Table 4, Hmax = Hmax, p − Hmax, m and tmax = tmax, p − tmax, m , where Hmax and tmax denote the maximum wave amplitude at Gauge 22 and the corresponding occurrence time, respectively; and the subscripts “p” and “m” denote the numerical prediction and the measured data, respectively. Note that the three numerical models use the same grid resolution in the simulations (x = y = 0.05 m). As shown in Table 4, in terms of the quantitative accuracy when predicting the maximum wave amplitude and the wave phase, the proposed SNH-NPC model gives the most accurate results. We also test the well-balanced property of the SNH-NPC model. In this test, all of the simulation parameters are the same as above except that initially, the flow is at rest and a constant water surface

In this study, we proposed a 2D SNH-NPC model based on the NPCM solving the shallow-water equations with non-hydrostatic pressure corrections. The HLLC approximate Riemann solver is used to calculate the numerical flux at the cell interface in the first (hydrostatic) sub-step calculations. We assessed the performance of the proposed numerical model using various benchmark tests, where the results showed that the model is accurate, robust, and wellbalanced. Comparisons against the analytic solutions, experimental data, and available numerical results demonstrated that the proposed model describes the wave dispersivity well and it accurately captures the wet–dry fronts. By using a shock-capturing numerical scheme, the proposed model simulates the wave breaking and reverse wave breaking processes reasonably well without requiring any tunable empirical coefficients, which are usually required by models solving the BTE. It should be noted that some inaccuracy was simulated in wave phase or wave amplitude (in the second, fourth, and fifth tests) using the proposed SNH-NPC model (as well as other SNH models known from the literature). We show here that this inaccuracy is not caused by the splitting error in the proposed model. In Appendix A, we present a SNH-IPC model based on the same SWE solver as described in Section 2.2.1. Fig. 13(a) and (b) plot, respectively, the water surface profiles computed by the SNH-NPC and SNH-IPC models, as well as

200

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

by the SNH-FSM and SNH-IPC models, for the solitary wave propagating over a constant water depth channel test with H/d = 0.1 (see Section 3.1). It is obvious that the differences in the results simulated by these three models are negligible, regardless of the mesh-cell numbers used. Similar comparisons were made for the other numerical tests considered in this study, and lead to the same conclusions (the relevant results are not shown). We make some further remarks about these comparisons. First, the results by the SNH-NPC and SNHFSM models agree well, thereby suggesting that the assumption described by Eq. (22) is reasonable and that the correction of the intermediate value of η (obtained in the hydrostatic sub-step) omitted in a SNH-FSM model is small (e.g., see Walters (2005); Zijlema and Stelling (2005)). Second, for the SNH-NPC and SNH-FSM models, the hydrostatic step is of the second-order accuracy in time, whereas the order of accuracy for the whole scheme may be reduced due to the splitting error (on the order O(t )) since the advection and pressure gradient terms do not commute (Zijlema and Stelling, 2005). Nevertheless, the magnitude of the splitting error for these two models is found to be negligible, as suggested by comparisons here. It is interesting to mention that Casulli and Stelling (1998), modeling short-wave motions with a FSM-based 3D non-hydrostatic model, reported significant wave damping, which was thought to be caused by the splitting error (Stelling and Zijlema, 2003). However, as a rather small value of t (Cr ≈ 0.02) was used by them, the splitting error (on the order of O(t )) is unlikely to be important and we think that their reported damping might be due to the inappropriate treatment of the pressure boundary condition at the free surface, because the zero pressure condition was enforced at the uppermost cell centers (half vertical grid below the free surface). In this sense, the Kellerbox scheme employed by Stelling and Zijlema (2003) (also used in the present study), which enforces an exact zero pressure condition at the free surface, appears to be very promising for representing the wave dispersivity and for maintaining the wave shape. While the inaccuracy in the predicted wave phase or wave amplitude was not caused by the splitting error in the model developed, it still can be caused by underestimation of wave dispersion in a SNH model (Stelling and Zijlema, 2003). In addition, the flow motions near the free surface are violent and complex when the wave breaks (Hu et al., 2012; Lu et al., 2013), so simple assumptions of linear vertical profile for the non-hydrostatic pressure and vertical velocity might be inappropriate. Using more layers in the vertical direction might lead to more accurate predictions of the wave phase and the wave breaking process; we refer the reader to Stelling and Zijlema (2003) and Bai and Cheung (2013) for studies using multi-layer models. It should be noted that only one particular SWE solver based on the finite-volume method was considered here for the SNH-NPC model. Any other finite-volume method-based SWE solver can be used for that relying on the methods described in this study. Since the numerical scheme here is based on the collocated grid, it can easily be extended to unstructured meshes.

using the θ -method proposed by Casulli (1990)

1 t

$

t+t

t

b dt ≈ θ bk+1 + (1 − θ )bk = bk + θ b ,

(A.1)

where b = bk+1 − bk . Then, the time-integration of the b related terms in Eq. (4), for instance, in Eq. (4b), can be approximated as

1 t

$

t+t



t

 ∂ bk ∂ b −I − Kxl b dt ≈ −Ik − Kxkl bk ∂ xl ∂ xl ∂ b −θ Ik − θ Kxkl b . ∂ xl

(A.2)

Hence, we can obtain the intermediate solution Um by using the SWE model described in Section 2.2.1 with a new source term given by

⎡

0

⎢ ⎢

⎤

⎥ ∂ zb gn2 qx |q| k ∂ bk ⎥ − − Kxk bk ⎥. −I 7/3 ∂x ∂x h ⎥ ⎣ ⎦ ∂ zb gn2 qy |q| k ∂ bk − − Kyk bk − gη −I 7/3 ∂y ∂y h

− gη S=⎢ ⎢

(A.3)

After Um is obtained, the intermediate vertical velocity wm is evalb uated by Eq. (28) and wm s is computed from

(ws )m + (wb )m − (ws )ki, j − (wb )ki, j (b )ki, j i, j i, j . = 2t ρ hm i, j

(A.4)

In the second step of the SNH-IPC method, Um and the vertical velocity are corrected to achieve a divergence-free velocity field by solving the following semi-discretized equations

(qx )k+1 − (qx )m (b )i+1, j − (b )i, j i+1/2, j i+1/2, j k = − θ Ii+1/2, j t x − θ (Kx )ki+1/2, j (b )i+1/2, j ,

(A.5a)

(qy )k+1 − (qy )m (b )i, j+1 − (b )i, j i, j+1/2 i, j+1/2 = − θ Ii,k j+1/2 t y − θ (Ky )ki, j+1/2 (b )i, j+1/2 ,

(A.5b)

Acknowledgments This study is part of a research project sponsored by the National Natural Science Foundation of China (Grant no. 51409195), the China Postdoctoral Science Foundation (Grant no. 2014M550408), and by the Fundamental Research Funds for the Central Universities (Grant No. 2042014kf0068). We are very grateful to the anonymous reviewers for their constructive comments and suggestions, which significantly improved the quality of this manuscript.

wk+1 − wm i, j i, j

t

=θ

(b )i, j , ρ hm i, j

(A.5c)

as well as the discretized continuity equation given by Eq. (25). Eqs. (A.5a)–(A.5c) share the same form as Eqs. (21a), (21b), and (24), so we can derive the following elliptic equation for the pressure correction b by procedures similar to those described in Section 2.2.2

Appendix A. SNH-IPC model In this section, we briefly describe a SNH-IPC model for solving the SNHE. In this model, the b can be integrated in time as follows

Aci, j (b )i−1, j + Bci, j (b )i, j + Ci,c j (b )i+1, j + Dci, j (b )i, j−1 + Ei,c j (b )i, j+1 = Si,c j ,

(A.6)

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

with the coefficients given by k − x(Kx )ki−1/2, j θ t 2Ii−1/2, j , hm 2x2 i−1/2, j

Aci, j = − Bci, j = +

2θ t θ t + 2 2  x2 ρ(hm ) i, j



k 2Ii+1/2, − x(Kx )ki+1/2, j j

hm i+1/2, j

+

k 2Ii−1/2, + x(Kx )ki−1/2, j j



hm i−1/2, j

  θ t 2Ii,k j+1/2 − y(Ky )ki, j+1/2 2Ii,k j−1/2 + y(Ky )ki, j−1/2 + , hm hm 2y2 i, j+1/2 i, j−1/2

Ci,c j = −

k + x(Kx )ki+1/2, j θ t 2Ii+1/2, j , 2 hm 2x i+1/2, j

Dci, j = −

θ t 2Ii,k j−1/2 − y(Ky )ki, j−1/2 , hm 2y2 i, j−1/2

Ei,c j = −

θ t 2Ii,k j+1/2 + y(Ky )ki, j+1/2 , hm 2y2 i, j+1/2

Si,c j = −

− um um i+1/2, j i−1/2, j

x

−

2(wb )i,k+1 − (ws )m − (wb )m vm − vm i, j i, j i, j+1/2 i, j−1/2 j + . y hm i, j

(A.7) Eq. (A.6) can be solved by employing the same method used to k+1 k+1 are comsolve Eq. (26). After b is obtained, qk+1 x , qy , and ws k+1 is updated by Eq. (19a). Eqs. (A.6) puted using Eq. (A.5). Finally, η and (26) are of the same form, so we can easily prove that the SNHIPC model is also well-balanced by employing the analysis method described in Section 2.2.2. Note that we also tested the following expression to approximate the left part of Eq. (A.2)

 ∂ b − Kxl b dt ∂ xl t  ∂ bk 1  k  1 ≈ − Ik + I(1) − K + Kx(l1) bk 2 ∂ xl 2 xl  ∂ b θ  k  θ − Ik + I(1) − K + Kx(l1) b , 2 ∂ xl 2 xl

1 t

$

t+t



−I

(A.8)

and thus the latest values of h and η were used in the calculations of I and Kxl (l = 1, 2) during each Runge–Kutta stage calculation, but we found that the difference in results was negligible. Thus, we employed Eq. (A.2) in this study since I and Kxl only need to be calculated once during each time step. The parameter θ is set to 0.6 in this study for numerical stability. References Ai, C., Jin, S., 2012. A multi-layer non-hydrostatic model for wave breaking and run-up. Coast. Eng. 62, 1–8. Audusse, E., Bouchut, F., Bristeau, M.-O., Klein, R., Perthame, B., 2004. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM J. Sci. Comput. 25 (6), 2050–2065. Bai, Y., Cheung, K.F., 2013. Dispersion and nonlinearity of multi-layer non-hydrostatic free-surface flow. J. Fluid Mech. 726, 226–260. Behrens, J., 2008. Tsunami-unstructured mesh finite element model for the computation of tsunami scenarios with inundation. In: Proceedings of NAFEMS Seminar, Wiesbaden, Germany, pp. 1–9. Briggs, M.J., Synolakis, C.E., Harkins, G.S., Green, D.R., 1995. Laboratory experiments of tsunami runup on a circular island. Pure Appl. Geophys. 144 (3-4), 569–593. Brocchini, M., 2013. A reasoned overview on Boussinesq-type models: the interplay between physics, mathematics and numerics. Proc. R. Soc. A 469 (2160), 1–27. Brufau, P., Vázquez-Cendón, M.E., García-Navarro, P., 2002. A numerical model for the flooding and drying of irregular domains. Int. J. Numer. Methods Fluids 39 (3), 247– 275. Casulli, V., 1990. Semi-implicit finite difference methods for the two-dimensional shallow water equations. J. Comput. Phys. 86, 56–74. Casulli, V., Stelling, G.S., 1998. Numerical simulation of 3D quasi-hydrostatic, freesurface flows. J. Hydraul. Eng., ASCE 124 (7), 678–686. Chorin, A.J., 1968. Numerical solution of the Navier–Stokes equations. Math. Comput. 22, 745–762.

201

Cui, H., Pietrzak, J.D., Stelling, G.S., 2010. A finite volume analogue of the P1NC-P1 finite element: with accurate flooding and drying. Ocean Modell. 35, 16–30. Cui, H., Pietrzak, J.D., Stelling, G.S., 2012. Improved efficiency of a non-hydrostatic, unstructured grid, finite volume model. Ocean Modell. 54-55, 55–67. Dalrymple, R.A., Rogers, B.D., 2006. Numerical modeling of water waves with the SPH method. Coast. Eng. 53, 141–147. Fiedler, F.R., Ramirez, J.A., 2000. A numerical method for simulating discontinuous shallow flow over an infiltrating surface. Int. J. Numer. Methods Fluids 32, 219– 239. Fraccarollo, L., Toro, E.F., 1995. Experimental and numerical assessment of the shallow water model for two-dimensional dam-break type problems. J. Hydraulic Res. 33 (6), 843–864. Gobbi, M.F., Kirby, J.T., 1996. A fourth order Boussinesq-type wave model. In: Proceedings of 25th Conference on Coastal Engineering, Orlando, Florida, pp. 1116–1129. Goda, K., 1979. A multistep technique with implicit difference schemes for calculating two- or three-dimensional cavity flows. J. Comput. Phys. 30, 76–95. Grilli, S.T., Ioualalen, M., Asavanant, J., Shi, F., Kirby, J.T., Watts, P., 2007. Source constraints and model simulation of the December 26, 2004, Indian Ocean tsunami. J. Waterway, Port, Coast. Ocean Eng. 133 (6), 414–428. Guermond, J., Minev, P., Shen, J., 2006. An overview of projection methods for incompressible flows. Comput. Methods Appl. Mech. Eng. 195, 6011–6045. Guo, X., Kang, L., Jiang, T., 2013. A new depth-integrated non-hydrostatic model for free surface flows. Sci. China Technol. Sci. 56 (4), 824–830. Hu, P., Li, W., He, Z., Pähtz, T., Yue, Z., 2015. Well-balanced and flexible morphological modeling of swash hydrodynamics and sediment transport. Coast. Eng. 96, 27–37. Hu, Y., Guo, X., Lu, X., Liu, Y., Dalrymple, R.A., Shen, L., 2012. Idealized numerical simulation of breaking water wave propagating over a viscous mud layer. Phys. Fluids 24, 112104. Kazolea, M., Delis, A.I., 2013. A well-balanced shock-capturing hybrid finite volume– finite difference numerical scheme for extended 1D Boussinesq models. Appl. Numer. Math. 67, 167–186. Keulegan, G.H., Patterson, G.W., 1940. Mathematical theory of irrotional translation waves. J. Res. Natl. Bureau Stand. 24, 47–101. Li, Y., Raichlen, F., 2002. Non-breaking and breaking solitary wave run-up. J. Fluid Mech. 456, 295–318. Liang, Q., Borthwick, A.G.L., 2009. Adaptive quadtree simulation of shallow flows with wet–dry fronts over complex topography. Comput. Fluids 38, 221–234. Lin, P., Chang, K.-A., Liu, P.L.-F., 1999. Runup and rundown of solitary waves of sloping beaches. J. Waterway, Port, Coast. Ocean Eng. 125 (5), 247–255. Liu, P.L.-F., Cheng, Y., 2001. A numerical study of the evolution of a solitary wave over a shelf. Phys. Fluids 13 (6), 1660–1667. Lu, X., Dong, B., Mao, B., Zhang, X., 2012. Unstructured mixed grid and SIMPLE algorithm based model for 2D-SWE. Proc. Eng. 28, 117–121. Lu, X., Zhang, X., Dong, B., Liu, H., Mao, B., 2013. Physical based numerical schemes for the discretization of the sediment settling term. C. R. Mécanique 341 (7), 581–591. Ma, G., Shi, F., Kirby, J.T., 2012. Shock-capturing non-hydrostatic model for fully dispersive surface wave processes. Ocean Modell. 43-44, 22–35. Madsen, O.S., Mei, C.C., 1969. The transformation of a solitary wave over an uneven bottom. J. Fluid Mech. 39 (4), 781–791. Mccabe, M.V., Stansby, P.K., Apsley, D.D., 2013. Random wave runup and overtopping a steep sea wall: shallow-water and Boussinesq modelling with generalised breaking and wall impact algorithms validated against laboratory and field measurements. Coast. Eng. 74, 33–49. Roe, P.L., Baines, M.J., 1982. Algorithms for advection and shock problems. In: 4th GAMM Conference on Numerical Methods in Fluid Mechanics, Braunschweig, Vieweg, pp. 281–290. Saad, Y., 2003. Iterative methods for sparse linear systems, 2nd ed. SIAM. Shu, C.-W., Osher, S., 1988. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J. Comput. Phys. 77 (2), 439–471. Stelling, G.S., Zijlema, M., 2003. An accurate and efficient finite-difference algorithm for non-hydrostatic free-surface flow with application to wave propagation. Int. J. Numer. Methods Fluids 43, 1–23. Street, R.L., Burges, S.J., Whitford, P.W., 1968. The behavior of solitary waves on a steeped slope. Technical Report. Stanford University, Stanford, CA. Sussman, M., Puckett, E.G., 2000. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows. J. Comput. Phys. 162, 301–337. Synolakis, C.E., 1986. The runup of long waves. California Institute of Technology, Pasadena, CA (Ph.D. thesis). Tan, L.-W., Chu, V.H., 2010. Wave runup simulations using Lagrangian blocks on Eulerian mesh. J. Hydro-environ. Res. 3, 193–200. Titov, V.V., Synolakis, C.E., 1995. Modeling of breaking and nonbreaking long-wave evolution and runup using VTCS-2. J. Waterway, Port, Coast., Ocean Eng. 121 (6), 308– 316. Toro, E.F., Spruce, M., Speares, W., 1994. Restoration of the contact surface in the HLLRiemann solver. Shock Waves 4, 25–34. van de Vorst, H.A., 1992. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Comput. 13 (2), 631– 644. van Kan, J., 1986. A second-order accurate pressure-correction scheme for viscous incompressible flow. SIAM J. Sci. Stat. Comput. 7 (3), 870–891. Walters, R.A., 2005. A semi-implicit finite element model for non-hydrostatic (dispersive) surface waves. Int. J. Numer. Methods Fluids 49, 721–737. Wei, Z., Jia, Y., 2013. A depth-integrated non-hydrostatic finite element model for wave propagation. Int. J. Numer. Methods Fluids 73, 976–1000.

202

X. Lu et al. / Ocean Modelling 96 (2015) 187–202

Wei, Z., Jia, Y., 2014. Simulation of nearshore wave processes by a depth-integrated non-hydrostatic finite element model. Coast. Eng. 83, 93–107. Yamazaki, Y., Kowalik, Z., Cheung, K.F., 2009. Depth-integrated, non-hydrostatic model for wave breaking and run-up. Int. J. Numer. Methods Fluids 61 (5), 473–497. Yu, B., Kawaguchi, Y., Tao, W., Ozoe, H., 2002. Checkerboard pressure predictions due to the underrelaxation factor and time step size for a nonstaggered grid with momentum interpolation method. Numer. Heat Transfer, Part B 41, 85–94.

Zelt, J.A., 1991. The run-up of nonbreaking and breaking solitary waves. Coast. Eng. 15, 205–246. Zijlema, M., Stelling, G.S., 2005. Further experiences with computing non-hydrostatic free-surface fows involving water waves. Int. J. Numer. Methods Fluids 48, 169–197.