Improvement of Non-Hydrostatic Hydrodynamic Solution Using a Novel Free-Surface Boundary Condition

Abstract

Hydrodynamic models based on the RANS equation are well-established tools to simulate three-dimensional free surface flows in large aquatic ecosystems. However, when the ratio of vertical to horizontal motion scales is not small, a non-hydrostatic approximation is needed to represent these processes accurately. Increasing efforts have been made to improve the efficiency of non-hydrostatic hydrodynamic models, but these improvements require higher implementation and computational costs. In this paper, we proposed a novel free-surface boundary condition based on a fictional sublayer at the free-surface (FSFS). We applied the FSFS approach at a finite difference numerical discretization with a fractional step framework, which uses a Neumann type of boundary condition to apply a hydrostatic relation in the top layer. To evaluate the model performance, we compared the Classic Boundary Condition Approach (CBA) and the FSFS approach using two numerical experiments. The experiments tested the model’s phase error, capability in solving wave celerity and simulate non-linear wave propagation under different vertical resolution scenarios. Our results showed that the FSFS approach had a lower phase error (2 to 5 times smaller) than CBA with a little additional computational cost (ca. 7% higher). Moreover, it can better represent wave celerity and frequency dispersion with 2 times fewer layers and low mean computational cost (CBA $\delta$t = 2.62 s and FSFS $\delta$t = 1.22 s).

1. Introduction

Hydrodynamic models based on the Reynolds Averaged Navier-Stokes (RANS) equation are well-established tools to simulate three-dimensional free surface flows in large aquatic ecosystems, such as lakes, estuaries, reservoirs, and coastal zones [1,2,3,4,5]. These models usually are based on the hydrostatic assumption of the pressure distribution, which is applied satisfactorily to large shallow water ecosystems with relatively low computational cost algorithms [6]. However, when the ratio of vertical to horizontal motion scales is not small, a non-hydrostatic approximation is needed to model accurately three-dimensional free surface flows.

Few alternatives were explored to improve efficiency and reduce the computational cost of non-hydrostatic hydrodynamic models. Most of them are dedicated to improving the model’s ability to solve the elliptic equation to non-hydrostatic pressure through more suitable boundary conditions and use a vertical momentum discretization less dependent on the velocity vertical profile [7,8,9,10,11,12]. Some approaches use a hybrid model (hydrostatic and non-hydrostatic), minimizing the computational cost by identifying regions where the hydrostatic assumption may be applied [6,13,14]. Optimization in source code processing is also a feasible solution, like a parallel computational algorithm, which allocated the loop calculations on the distributed threads [6,15]. The implementation costs are associated with these issues and aspects related to assumptions adopted in models.

As for the assumptions, two different types of boundary conditions may be applied for non-hydrostatic pressure at the free-surface: (a) a homogenous Dirichlet condition, where pressure is set equal zero at the free-surface (e.g., [6,7,16,17]); and (b) a Neumann condition, usually a hydrostatic relation to approximate pressure’s value equal to zero at the free-surface (e.g., [14,18,19,20]). Both types of free-surface boundary conditions were well discussed in Bergh and Berntsen [21]. However, many numerical models solve Poisson’s equation for pressure at the center of the computational cell (e.g., [8,11,14,22,23]). This assumption, without a suitable treatment, may incorrectly estimate the non-hydrostatic boundary condition to the center of the top layer instead of at the free-surface. This assumption becomes more inaccurate as the distance between the center of the top layer and the free surface increases. Hence, a high vertical resolution (1–20 layers) is needed to achieve accurate results for different studies case [7,8,9,10,11,24].

In addition to using a suitable vertical discretization, [8,9] showed that cell-center non-hydrostatic models might have a substantial accumulated phase error, due to the wrong computations in wave celerity. Many efforts were carried out in order to properly overcome this issue, of which we can highlight: (a) an implementation of edge-based non-hydrostatic pressure (e.g., [7,25]), (b) the use of an integration method to estimate the pressure at the free-surface based on pressure in the center of the top layer (e.g., [8,10,11]), and (c) the use of a piecewise linear profile of the non-hydrostatic pressure to estimate the pressure at the free-surface (e.g., [12]). For the existing cell-centered models, the previous algorithms may demand a substantial implementation cost due to changes in the vertical momentum discretization and in the treatment needed to address the non-hydrostatic pressure at the free-surface boundary condition adequately. Thus, in this paper, we proposed a low implementation cost method to approach the free-surface non-hydrostatic pressure in hydrodynamics models properly. The technique is a novel free-surface boundary condition based on a fictional sublayer at the free-surface ($FSFS$), to reach satisfactory results without changing the momentum equations.

We applied the FSFS, and the CBA approaches at a finite difference discretization with a fractional step framework based on [23], which uses a Neumann type of boundary condition that applies a hydrostatic relation in the top layer. To evaluate the model performance, we used two widely applied numerical models benchmarks [8,11,22,23,26,27,28], selected to test our algorithm in two different purposes: (a) a standing wave in a three-dimensional closed basin to test the model’s capability in solving wave celerity under different vertical resolutions; and (b) a wave propagation over a submerged bar to validate the proposed boundary condition at the free-surface and evaluated the effect of vertical resolution under a non-linear wave propagation.

2. Methods

2.1. Governant Equations

The RANS equations are used to describe three-dimensional free-surface flows. These equations express the physical principle of volume, mass, and momentum conservation. The momentum equations for an incompressible fluid have the following form:

$$\begin{array}{c}\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}-fv=-\frac{\partial p_a}{\partial x}-g\frac{\partial \eta }{\partial x}\hfill \\ -g\frac{\partial }{\partial x}\left[{\int }_z^{\eta }\frac{\rho -{\rho }_0}{{\rho }_0}d\zeta \right]-\frac{\partial q}{\partial x}+{\nu }^h\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)+\frac{\partial }{\partial z}\left({\nu }^v\frac{\partial u}{\partial z}\right),\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill & \frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}-fu=-\frac{\partial p_a}{\partial y}-g\frac{\partial \eta }{\partial y}\hfill \\ \hfill & -g\frac{\partial }{\partial y}[{\int }_z^{\eta }\frac{\rho -{\rho }_0}{{\rho }_0}d\zeta ]-\frac{\partial q}{\partial y}+{\nu }^h\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)+\frac{\partial }{\partial z}\left({\nu }^v\frac{\partial v}{\partial z}\right),\hfill \end{array}$$

(2)

$$\frac{\partial w}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z}=-\frac{\partial q}{\partial z}+{\nu }^h(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2})+\frac{\partial }{\partial z}\left({\nu }^v\frac{\partial w}{\partial z}\right),$$

(3)

where $u(x,y,z,t)$, $v(x,y,z,t)$, and $w(x,y,z,t)$ are the velocity components in the horizontal (x and y) and vertical (z) directions $(m/s)$, respectively; ${\nu }^h$ e ${\nu }^v$ are the horizontal and vertical turbulent eddy viscosity coefficients $(m/s)$, respectively; t is the time $\left(s\right)$; $p_a(x,y,z,t)$ is the atmospheric pressure normalized by fluid density $(Pa.m^3/kg)$; $\eta$ is the free surface elevation $\left(m\right)$ from a water level References; $q(x,y,z,t)$ denotes the non-hydrostatic pressure component normalized by fluid density $(Pa.m^3/kg)$; f is the Coriolis parameter $\left(s^{-1}\right)$; and g is the gravitational acceleration $(m/s^2)$. The second and third terms on the right-hand side of Equations (1) and (2) represent the barotropic and the baroclinic contributions to the hydrostatic pressure.

When a simple hydrostatic approach is considered, Equation (3) is neglected and q is assumed to be equal to zero in Equations (1) and (2). In this case, it is assumed that vertical acceleration do not have a significant effect in the velocity field in comparison with horizontal acceleration.

The volume conservation is expressed by the incompressibility condition and the continuity equation, given by:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0.$$

(4)

Integrating Equation (4) over depth yields to the following equation:

$${\int }_{-h}^{\eta }\left[\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right]dz={\int }_{-h}^{\eta }\frac{\partial u}{\partial x}dz+{\int }_{-h}^{\eta }\frac{\partial v}{\partial y}dz+{\int }_{-h}^{\eta }\frac{\partial w}{\partial z}dz=0,$$

(5)

where “h” is the bathymetry measured from the theoretical undisturbed water surface (zero referential). Using the Leibniz integration rule, in each direction, in Equation (5) and using a kinematic condition at the free-surface leads to the following free-surface equation:

$$\frac{\partial \eta }{\partial t}+\frac{\partial }{\partial x}{\int }_{-h}^{\eta }udz+\frac{\partial }{\partial y}{\int }_{-h}^{\eta }vdz=0,$$

(6)

The boundary conditions were implemented under the assumption of “free-slip” boundaries. The Dirichlet and Neumann conditions were assigned to represent the normal and tangential velocities in the solid boundaries, respectively.

The tangential stress boundary conditions for the momentum equations (Equations (1) and (2)) are specified at the free-surface by the prescribed wind stresses, which can be approximated as:

$${\nu }^v\frac{\partial u}{\partial z}={\gamma }_T(u_a-u),{\nu }^v\frac{\partial v}{\partial z}={\gamma }_T(v_a-v),atz=\eta ,$$

(7)

where $u_a$ and $v_a$ are the horizontal wind velocity components, and ${\gamma }_T$ is a non-negative wind stress coefficient. The bottom friction is specified by:

$${\nu }^v\frac{\partial u}{\partial z}={\gamma }_{B}u,{\nu }^v\frac{\partial v}{\partial z}={\gamma }_{B}v,atz=-h,$$

(8)

where ${\gamma }_B$ is a non-negative bottom friction coefficient.

2.2. Grid and Variables Locations

The computational grid can be described as a generic unstructured orthogonal grid, having a set of non-overlapping convex $N_p$ elements, each having an arbitrary number of sides $S_i$$\ge 3,i=1,2,\dots ,N_p$ (Figure 1). Let $N_s$ be the total number of sides in the grid. The length of each side is ${\lambda }_j,j=1,2,\dots ,N_p$. The vertical faces of the i-th element are identified by an index $j_{(i,l)},l=1,2,\dots ,S_i$, so that $1\le j_{(i,l)}\le N_s$. Similarly, the two polygons which share the j-th vertical face of the grid are identified by the indices $i_{(j,1)}$ and $i_{(j,2)}$, so that $1\le i_{(j,1)}\le N_p$ and $1\le i_{(j,2)}\le N_p$. The nonzero distance between centers of two adjacent polygons which share the $j\text{-}\mathrm{th}$ side is denoted with ${\delta }_j$.

Along the vertical direction, a simple finite difference discretization, not necessarily uniform, is adopted. By denoting with $\mathsf{\Delta }{z}_{k+\frac{1}{2}}$ a given top computational cell level surface, the vertical discretization step is defined by:

$$\mathsf{\Delta }{z}_k=\mathsf{\Delta }{z}_{k+\frac{1}{2}}-\mathsf{\Delta }{z}_{k-\frac{1}{2}}k=1,2,\cdots ,N_s.$$

(9)

The three-dimensional spatial discretization consists of elements whose horizontal faces are the polygons of a given orthogonal grid, represented by the layers at $k+\frac{1}{2}$ (upper face) or $k-\frac{1}{2}$ (bottom face), whose height, for each layer, is $\mathsf{\Delta }{z}_k$. The water surface elevation, ($\eta$), is located at the barycenter of the upper horizontal face for each i-th element. The velocity component normal to each horizontal face is assumed to be constant over the face of each computational cell and is defined at the point of intersection between the face and the segment joining the centers of the two prisms which share the face. The non-hydrostatic pressure component $q_{i,k}^n$ is located at the center of the i-th computational cell, halfway between $\mathsf{\Delta }{z}_{k+\frac{1}{2}}$ and $\mathsf{\Delta }{z}_{k-\frac{1}{2}}$. Finally, the water depth $h_j$ is specified and assumed constant on each vertical face of an element.

2.3. Numerical Approximation

We used a semi-implicit method of finite volume ($\theta$-Method [29]) with an Eulerian–Lagrangian Method [30] to solve the convective and viscous terms of the RANS equations, which applies a quadratic interpolation [31] to estimate the velocity field at the end of backtracking process (multi-step backward Euler with ten sub-time steps). A fractional-step framework [23] is used to solve the pressure component by splitting into hydrostatic and non-hydrostatic parts. The algorithm procedures take the following steps:

• Definition of initial parameters, initial conditions, and boundary conditions.
• Solution of convective terms using the Eulerian–Lagrangian Method.
• Determination of the provisional free-surface elevation ($\tilde{\eta }$) through the preconditioned conjugate gradient iterations until the residual norm is smaller than a given tolerance ${ϵ}_q$.
• Numeric solution of provisional velocity field ($\tilde{u}$ and $\tilde{w}$).
• Solution of non-hydrostatic pressure (q) through the preconditioned conjugate gradient iterations until the residual norm is smaller than a given tolerance ${ϵ}_q$.
• Numeric correction of velocity field and free surface elevation.

A complete description of the numerical solution may be found in [23]. Here we focus only in described the conventional non-hydrostatic pressure discretization separately, the given boundary condition, and how applied the $FSFS$ boundary condition in the conventional solution to improve the non-hydrostatic hydrodynamic solution.

2.3.1. Non-Hydrostatic Pressure Discretization

A correction of provisional velocity field $({\tilde{u}}_{j,k}^{n+1}$,${\tilde{w}}_{j,k}^{n+1})$ are computed after including the non-hydrostatic pressure terms, specifically:

$$u_{j,k}^{n+1}={\tilde{u}}_{j,k}^{n+1}-\theta \frac{\mathsf{\Delta }t}{{\delta }_j}\left({\tilde{q}}_{i\left(j,r\right),k}^{n+1}-{\tilde{q}}_{i\left(j,l\right),k}^{n+1}\right),$$

(10)

$$w_{j,k}^{n+1}={\tilde{w}}_{j,k}^{n+1}-\theta \frac{\mathsf{\Delta }t}{\mathsf{\Delta }{z}_{i,k+\frac{1}{2}}}({\tilde{q}}_{i,k+1}^{n+1}-{\tilde{q}}_{i,k+1}^{n+1}),$$

(11)

where the vertical space increment $\mathsf{\Delta }z$ is defined as the distance between two consecutive level surfaces, except near the bottom and near the free surface where $\mathsf{\Delta }z$ is the distance between a level surface and the bottom, or free-surface, respectively. The $\tilde{q}$ term denotes the non-hydrostatic pressure term, which in combination with the provisional free-surface elevation ($\tilde{\eta }$), gives the pressure:

$$p_{j,k}^{n+1}=g\left({\tilde{\eta }}_i^{n+1}-z_k\right)+{\tilde{q}}_{i,k}^{n+1},$$

(12)

where $z_k$ is the z-coordinate of the k-th horizontal level surface, and g is the gravity acceleration. In each computational cell bellow the free-surface, the finite volume discretized incompressibility condition is taken to be:

$${\displaystyle \sum _{L=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,L\right)}\mathsf{\Delta }{z}_{j\left(i,L\right),k}^n{u}_{j\left(i,L\right),k}^{n+1}+V_i\left(w_{i,k+\frac{1}{2}}^{n+1}-w_{i,k-\frac{1}{2}}^{n+1}\right)=0k=m,m+1,\dots ,M-1,$$

(13)

where $V_i$ is the area of the i-th polygon and $S_i$ is the number of sides for the i-th element. At the free-surface, the finite difference approximation of Equation (6) considering $w_{i,m-\frac{1}{2}}^{n+1}=0$ and using the incompressibility condition (13) is:

$$\begin{array}{c}{V}_i{\eta }_i^{n+1}=V_i{\eta }_i^n-\theta \mathsf{\Delta }t{\displaystyle \sum _{l=1}^{S_i}}\left[s_{i,l}{\lambda }_{j\left(i,L\right)}\mathsf{\Delta }{z}_{j\left(i,l\right),M}^n{u}_{j\left(i,l\right),M}^{n+1}\right]+\theta \mathsf{\Delta }t{V}_i{w}_{i,M-\frac{1}{2}}^{n+1}\hfill \\ -\left(1-\theta \right)\mathsf{\Delta }t{\displaystyle \sum _{l=1}^{S_i}}\left[s_{i,l}{\lambda }_{j\left(i,l\right)}{\displaystyle \sum _{k=m}^M}\mathsf{\Delta }{z}_{j\left(i,l\right),k}^n{u}_{j\left(i,l\right),k}^n\right],\hfill \end{array}$$

(14)

where $\theta$ is the implicitness factor, $s_{i,l}$ is sign function associated with the orientation of the normal velocity defined on the l side of an element i. Assuming that the pressure at the FSFS is hydrostatic, the pressure correction term (${\tilde{q}}_{i,M}^{n+1}$) is obtained by the following hydrostatic relation:

$$p_{j,M}^{n+1}=g({\eta }_i^{n+1}-z_M)=g({\tilde{\eta }}_i^{n+1}-z_M)+{\tilde{q}}_{i,k}^{n+1},$$

(15)

hence, a substitution the $p_{j,M}^{n+1}$ in the Equation (12) yields:

$$\begin{array}{c}{V}_i{\tilde{\eta }}_i^{n+1}=g{V}_i({\eta }_i^n-{\tilde{\eta }}_i^{n+1})-g\theta \mathsf{\Delta }t{\displaystyle \sum _{L=1}^{S_i}}\left[s_{i,L}{\lambda }_{j\left(i,L\right)}\mathsf{\Delta }{z}_{j\left(i,L\right),M}^n{u}_{j\left(i,L\right),M}^{n+1}\right]+g\theta \mathsf{\Delta }t{V}_i{w}_{i,M-\frac{1}{2}}^{n+1}\hfill \\ -g\left(1-\theta \right)\mathsf{\Delta }t{\displaystyle \sum _{L=1}^{S_i}}\left[s_{i,L}{\lambda }_{j\left(i,L\right)}{\displaystyle \sum _{k=m}^M}\mathsf{\Delta }{z}_{j\left(i,L\right),k}^n{u}_{j\left(i,L\right),k}^n\right].\hfill \end{array}$$

(16)

A system of equations to solve non-hydrostatic pressure is now derived by substituting the expressions for the new velocities from (10) and (11) into (13) and (16), respectively. The following finite difference equations are obtained:

$$\begin{array}{c}g{\theta }^2\mathsf{\Delta }{t}^2\left[{\displaystyle \sum _{l=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,l\right)}\mathsf{\Delta }{z}_{j\left(i,l\right),k}^n\frac{{\tilde{q}}_{i\left[j\left(i,l\right),1\right],k}^{n+1}-{\tilde{q}}_{i\left[j\left(i,l\right),2\right],k}^{n+1}}{{\delta }_{j\left(i,l\right)},k}+V_i\left(\frac{{\tilde{q}}_{i,k}^{n+1}-{\tilde{q}}_{i,k+1}^{n+1}}{\mathsf{\Delta }{z}_{i,k+\frac{1}{2}}^n}-\frac{{\tilde{q}}_{i,k-1}^{n+1}-{\tilde{q}}_{i,k}^{n+1}}{\mathsf{\Delta }{z}_{i,k-\frac{1}{2}}^n}\right)\right]\hfill \\ =g\theta \mathsf{\Delta }t{V}_i\left({\tilde{w}}_{i,k-\frac{1}{2}}^{n+1}-{\tilde{w}}_{i,k+\frac{1}{2}}^{n+1}\right)-g\theta \mathsf{\Delta }t{\displaystyle \sum _{l=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,l\right)}\mathsf{\Delta }{z}_{j\left(i,l\right),k}^n{\tilde{u}}_{j\left(i,l\right),k}^{n+1},k=m,m+1,\dots ,M-1,\hfill \end{array}$$

(17)

$$\begin{array}{c}g{\theta }^2\mathsf{\Delta }{t}^2\left[{\displaystyle \sum _{l=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,l\right)}\mathsf{\Delta }{z}_{j\left(i,l\right),M}^n\frac{{\tilde{q}}_{i\left[j\left(i,l\right),1\right],M}^{n+1}-{\tilde{q}}_{i\left[j\left(i,l\right),2\right],M}^{n+1}}{{\delta }_{j\left(i,l\right),M}}-V_i\frac{{\tilde{q}}_{i,M-1}^{n+1}-{\tilde{q}}_{i,M}^{n+1}}{\mathsf{\Delta }{z}_{i,M-\frac{1}{2}}^n}\right]+V_i{\tilde{q}}_{i,M}^{n+1}=\hfill \\ g\theta \mathsf{\Delta }t{V}_i{\tilde{w}}_{i,M-\frac{1}{2}}^{n+1}-g\theta \mathsf{\Delta }t{\displaystyle \sum _{l=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,l\right)}\left(\mathsf{\Delta }{z}_{FSFS}\right){\tilde{u}}_{j\left(i,l\right),M}^{n+1}+g{V}_i({\eta }_i^n-{\tilde{\eta }}_i^{n+1})-\hfill \\ g(1-\theta )\mathsf{\Delta }t{\displaystyle \sum _{l=1}^{S_i}}[s_{i,L}{\lambda }_{j\left(i,l\right)}{\displaystyle \sum _{k=m}^M}\mathsf{\Delta }{z}_{j\left(i,l\right),k}^n{u}_{j\left(i,l\right),k}^n,k=M.\hfill \end{array}$$

(18)

Once the non-hydrostatic pressure terms are computed, the velocities field are corrected by Equation (10), while vertical velocity can be estimated, equivalently, by Equation (11) or by the incompressibility condition (13) by setting $w_{i,m+\frac{1}{2}}^{n+1}=0$:

$$w_{i,k+\frac{1}{2}}^{n+1}=w_{i,k-\frac{1}{2}}^{n+1}-\frac{1}{V_i}{\displaystyle \sum _{L=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,L\right)}\mathsf{\Delta }{z}_{j\left(i,L\right),k}^n{u}_{j\left(i,L\right),k}^{n+1};k=m,m+1,\dots ,M-1.$$

(19)

This equation guarantees that the resulting velocity field is exactly discrete divergence-free [23], so we used it to compute the vertical velocity components.

The final free surface elevation is obtained by the hydrostatic relation (15) as follows:

$${\eta }_i^{n+1}={\tilde{\eta }}_i^{n+1}+\frac{{\tilde{q}}_{i,M+1}^{n+1}}{g}.$$

(20)

Finally, the non-hydrostatic pressure component can be obtained by:

$$q_{i,k}^{n+1}={\tilde{q}}_{i,k}^{n+1}-{\tilde{q}}_{i,M+1}^{n+1};k=m,m+1,\dots ,M,M+1,$$

(21)

making the non-hydrostatic pressure at free-surface equal to zero (i.e., for $k=M+1$).

2.3.2. Free-Surface Boundary Condition Treatment

To apply the Fictional Sublayer at the Free-Surface ($FSFS$) into the existing computational domain is only required one additional numerical vertical layer at the top of the computational domain, which does not account to the computational domain. As the height of the $FSFS$ is assumed to be equal to zero, the model can always be guaranteed the hydrostatic relation at the free-surface, independent of the number of vertical layers. As the pressure is estimated at the center of the layer, which can be further away from the surface the smaller the number of layers in CBA, the proposed FSFS condition solves these problems and provides a more physically consistent numerical solution, since the non-hydrostatic pressure condition is set in the fictional sublayer at the free-surface. To use this method only is required a simple adaption in the numerical Equations (17) and (18), considering the position of the layer:

(i)

For bottom and middle layers (i.e., $k=m$ to $M-1$), the Equation (17) is applied using its original form;

(ii)

For the top layer ($k=M$), Equation (17) is adapted to take into account the influence of FSFS height ($\mathsf{\Delta }{z}_{FSFS}$) in $\mathsf{\Delta }{z}_{i,M+\frac{1}{2}}^n$ and $\mathsf{\Delta }{z}_{j\left(i,l\right),k}^n$:

$$\begin{array}{c}g{\theta }^2\mathsf{\Delta }{t}^2\left[{\displaystyle \sum _{l=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,l\right)}\left(\mathsf{\Delta }{z}_{j\left(i,l\right),M}^n-\mathsf{\Delta }{z}_{FSFS}\right)\frac{{\tilde{q}}_{i\left[j\left(i,l\right),1\right],M}^{n+1}-{\tilde{q}}_{i\left[j\left(i,l\right),2\right],M}^{n+1}}{{\delta }_{j\left(i,l\right)},M}\right]+\hfill \\ \left[V_i\left(\frac{{\tilde{q}}_{i,M}^{n+1}-{\tilde{q}}_{i,M+1}^{n+1}}{\mathsf{\Delta }{z}_{i,M+\frac{1}{2}}^n}-\frac{{\tilde{q}}_{i,M-1}^{n+1}-{\tilde{q}}_{i,M}^{n+1}}{\mathsf{\Delta }{z}_{i,M-\frac{1}{2}}^n}\right)\right]=g\theta \mathsf{\Delta }t{V}_i\left({\tilde{w}}_{i,M-\frac{1}{2}}^{n+1}-{\tilde{w}}_{i,M+\frac{1}{2}}^{n+1}\right)-\hfill \\ g\theta \mathsf{\Delta }t{\displaystyle \sum _{l=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,l\right)}\left(\mathsf{\Delta }{z}_{j\left(i,l\right),M}^n-\mathsf{\Delta }{z}_{FSFS}\right){\tilde{u}}_{j\left(i,l\right),M}^{n+1},k=M,\hfill \end{array}$$

(22)

$$z_{i,M+\frac{1}{2}}^n=\frac{1}{2}[(\mathsf{\Delta }{z}_{i,M}^n-\mathsf{\Delta }{z}_{FSFS})+\mathsf{\Delta }{z}_{FSFS}]=\frac{1}{2}(\mathsf{\Delta }{z}_{i,M}^n)$$

(23)

$$\mathsf{\Delta }{z}_{i,M-\frac{1}{2}}^n=\frac{1}{2}[(\mathsf{\Delta }{z}_{i,M}^n-\mathsf{\Delta }{z}_{FSFS})+\mathsf{\Delta }{z}_{i,M-1}^n];$$

(24)

(iii)

For the FSFS layer ($k=M+1$), Equation (18) is adapted to take into account the $FSFS$ height and the velocity field in the layer M. Preliminary simulations showed that making $\mathsf{\Delta }{z}_{FSFS}=0$ a stable solution is achieved for any vertical discretization:

$$\begin{array}{c}g{\theta }^2\mathsf{\Delta }{t}^2\left[{\displaystyle \sum _{l=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,l\right)}\left(\mathsf{\Delta }{z}_{FSFS}\right)\frac{{\tilde{q}}_{i\left[j\left(i,l\right),1\right],M+1}^{n+1}-{\tilde{q}}_{i\left[j\left(i,l\right),2\right],M+1}^{n+1}}{{\delta }_{j\left(i,l\right),M}}-V_i\frac{{\tilde{q}}_{i,M}^{n+1}-{\tilde{q}}_{i,M+1}^{n+1}}{\mathsf{\Delta }{z}_{i,\left(M+1\right)-\frac{1}{2}}^n}\right]+V_i{\tilde{q}}_{i,M+1}^{n+1}=\hfill \\ g\theta \mathsf{\Delta }t{V}_i{\tilde{w}}_{i,M-\frac{1}{2}}^{n+1}-g\theta \mathsf{\Delta }t{\displaystyle \sum _{l=1}^{S_i}}{s}_{i,L}{\lambda }_{j\left(i,l\right)}\left(\mathsf{\Delta }{z}_{FSFS}\right){\tilde{u}}_{j\left(i,l\right),M}^{n+1}+g{V}_i({\eta }_i^n-{\tilde{\eta }}_i^{n+1})\hfill \\ -g(1-\theta )\mathsf{\Delta }t{\displaystyle \sum _{l=1}^{S_i}}\left[s_{i,L}{\lambda }_{j\left(i,l\right)}\right]{\displaystyle \sum _{k=m}^M}\mathsf{\Delta }{z}_{j\left(i,l\right),k}^n{u}_{j\left(i,l\right),k}^n,k=M+1,\hfill \end{array}$$

(25)

where

$$\mathsf{\Delta }{z}_{i,(M+1)-\frac{1}{2}}^n=\frac{1}{2}[(\mathsf{\Delta }{z}_{i,M}^n-\mathsf{\Delta }{z}_{FSFS})+\mathsf{\Delta }{z}_{FSFS}]=\frac{1}{2}(\mathsf{\Delta }{z}_{i,M}^n).$$

(26)

2.4. Numerical Experiments

The proposed numerical approach was applied in two consolidated benchmarks, usually used to verification and validation of numerical models (e.g., [8,11,22,23,26,28]). Each numerical experiment has a different purpose, as follows:

a

Standing waves in a three-dimensional closed basin: This test was widely applied in the literature to verify the model’s ability to simulate 3D linear waves comparing the analytic solution with the numerical solution in regard to phase and amplitude representation [8,9,11,20]. We evaluated the model capability in calculating the wave celerity and frequency wave dispersion with the Classic Boundary Approach (hereafter named CBA) and with the proposed FSFS boundary condition. We used six different vertical resolutions (20~5 layers), as most of the previous studies do, and since previous analyses showed that more than 20 layers do not have substantial improvement over 20 layers non-hydrostatic solution. We compare the free surface elevation cumulative phase error, the mean one time-step computational cost, the number of wave periods, and the relation with the free surface vertical velocity after 30 s of simulation in comparison with the analytical solution. We also compared the free surface elevation results with some metrics (RMSE, BIAS, Volume Error, KGE, NSE). At last, we statistically tested the residual series (the difference between analytic and simulated results) with the non-parametric Kruskal–Wallis test followed by a post-hoc Nemenyi to identify significant differences concerning the analytic results. The mean time of one time-step simulation was computed using an Intel® Xenon® CPU-E5-1620 3.7 GHz computer with 32 GB of RAM in a Fortran-based numerical model.

b

The wave propagation over a submerged bar: This test case was an experimental model idealized by [32], and was frequently used to validate numerical models (e.g., [7,8,24,26,28,33]). This experiment was used to evaluate accuracy to represent a non-linear wave pattern due to physical changes at the bottom, by comparing free surface elevation between the FSFS approach with a different vertical resolution between simulated and experimental results. To evaluate the model’s performance was used a few metrics (RMSE, BIAS, Volume Error, KGE, NSE) and the statistical test non-parametric Kruskal–Wallis test followed by a post-hoc Nemenyi test applied at residuals series.

3. Results

3.1. 3D Standing Waves in a Closed Basin

In this test, was analyzed a flow induced by an initial wave amplitude set to $0.1$ m in a closed cubic basin with 10 m of edge (Figure 2). To discretize the spatial domain, was adopt a regular with $0.5$ m resolution, resulting in 8.000 computational cells. The simulation is carried along 30 s, with a time-step size $\mathsf{\Delta }$t = $0.01$ s. The analytic solution of free-surface water elevation is given by:

$$\eta =Acos\left(k_{x}x\right)cos\left(k_{y}y\right)cos\left(2\pi \frac{t}{T}\right),$$

(27)

where t is the time (the initial condition in the free surface may be obtained by doing $t=0$); T is the wave period equal to 3.1 s, with the wave number $kx=ky=n/L$ with the total wave number $k=\sqrt{k_x^2+k_y^2}=0.44\frac{rad}{m}$.

The analytic solution for each velocity component is described as follows:

$$u=\frac{Ag{k}_x}{\omega }\frac{cosh\left[k_x(h+z)\right]}{cosh\left(k_{x}h\right)}sin\left(k_{x}x\right)cos\left(k_{y}y\right)sin\left(\omega t\right),$$

(28)

$$u=\frac{Ag{k}_y}{\omega }\frac{cosh\left[k_y(h+z)\right]}{cosh\left(k_{y}h\right)}cos\left(k_{x}x\right)sin\left(k_{y}y\right)sin\left(\omega t\right),$$

(29)

$$w=\frac{Ag{k}_x}{\omega }\frac{sinh\left[k_x(h+z)\right]}{cosh\left(k_{x}h\right)}cos\left(k_{x}x\right)cos\left(k_{y}y\right)sin\left(\omega t\right),$$

(30)

where $\omega$ represents the wave dispersion relation, given by:

$$\omega =\sqrt{gKtanh\left(Kh\right)}.$$

(31)

The study case applied different vertical resolution scenarios (i.e., 20, 16, 13, 10, 8, and 5 vertical layers), with both $CBA$ and the proposed $FSFS$ condition (called methods). We evaluated the cumulative phase error of a 3D standing wave by comparing the model outcomes (i.e., free-surface elevation and free-surface vertical velocity component) at $x=y=0.25$ m with the analytical solution. The performance between boundary condition approaches was evaluated by comparing the free-surface elevation residuals. The performance between scenarios (considering a different number of vertical layers) was also evaluated by statistically tested the residues with a non-parametric Kruskal–Wallis test followed by a Nemenyi test since the residues series do not follow a normal distribution (Shapiro–Wilk test). The numeric experiment aims to identify a critical vertical resolution for the FSFS method, where the results are significantly different from the used best scenario for this benchmark (20 vertical layers).

In general, for both methods, the phase error (${\mathsf{\Phi }}_{\epsilon }$) increases over time step, and it becomes bigger with the reduction of the number of layers (

Table 1. Computational cost, Phases Error, and the number of wave periods between different methods and scenarios. The model was implemented with Fortran and simulated in a machine using an Intel R Xenon R CPU-E5-1620 3.7 GHz computer with 32 GB of RAM.

N°L	FSFS			CBA
	$\mathsf{\Delta }\mathit{t}\left(\mathit{s}\right)$	N° T	${\mathbf{\Phi }}_{\mathit{\epsilon }}$	$\mathsf{\Delta }\mathit{t}\left(\mathit{s}\right)$	N° T	${\mathbf{\Phi }}_{\mathit{\epsilon }}$
20-Layers	2.8	10	0.3	2.62	10	1.6
16-Layers	2.21	10	0.4	2.14	10	2
13-Layers	1.68	10	0.6	1.67	10	2.4
10-Layers	1.54	10	1	1.45	11	3.1
8-Layers	1.51	10	1.5	1.34	11	3.8
5-Layers	1.22	11	2.1	1.08	12	4.6

). For simulations CBA approach, these errors are higher; for instance, the 20-layer CBA simulation had comparable results with the 8-layers FSFS simulation (Figure 3). Furthermore, when phase error is critical, a decrease in the cumulative free surface elevation error in the low vertical resolution results occurs (less than ten) due to an increase in wave periods. Due to this, results with fewer layers appeared better than those with more layers, as can be seen by comparing FSFS-5L with FSFS-10L results (Figure 3). These effects complicate the direct comparison between methods; however, the cumulative residual free surface elevation error (Figure 4) can clarify the matter. Due to this, the calculated metrics (

Table 2. Metrics between the analytic and simulated results from the FSFS method for each scenario for the first 10 s of simulation.

Metrics	FSFS-20L	FSFS-16L	FSFS-13L	FSFS-10L	FSFS-8L	FSFS-5L
RMSE (mm)	7.40	11.52	18.29	30.53	45.30	89.20
BIAS (mm)	0.32	0.32	−0.77	−2.11	−3.54	−7.61
Error (%)	9.25	14.28	22.55	37.79	56.60	114.13
KGE	0.93	0.93	0.82	0.52	0.18	−0.88
NSE	0.99	0.97	0.93	0.81	0.57	−0.66

and

Table 3. Metrics between the analytic and simulated results from the CBA method for each scenario for the first 10 s of simulation.

Metrics	CBA-20L	CBA-16L	CBA-13L	CBA-10L	CBA-8L	CBA-5L
RMSE (mm)	47.06	57.68	68.75	83.48	95.63	107.80
BIAS (mm)	−4.35	−5.65	−6.88	−7.96	−7.57	−1.06
Error (%)	58.94	72.70	87.10	106.42	123.39	145.02
KGE	0.004	−0.30	−0.61	−0.91	−0.93	−0.20
NSE	0.54	0.31	0.017	−0.45	−0.90	−1.42

) was only comparable between 0 and 10 s of simulations (Figure 4). This effect can be better verified in CBA accumulated residual series when the graphic changes the slope, which occurs sooner as the number of layers reduced (Figure 4).

All statics metrics indicated a decrease of performance with a reduction of the number of layers. The metrics show that CBA-20L simulation had similar results as the FSFS-8L, as well as CBA-10L to FSFS-5L. The CBA method also had a decrease of performance with a reduction of vertical resolution (Table 3 and

Table 4. Nemenyi posthoc test comparing the FSFS residue series of the simulation with 20 to 5 vertical layers to identify the significative statistical difference between results.

N°-L	FSFS-20L	FSFS-16L	FSFS-13L	FSFS-10L	FSFS-8L
FSFS-16L	0.97833	-	-	-	-
FSFS-13L	0.24554	0.69401	-	-	-
FSFS-10L	<0.05	<0.05	0.05269	-	-
FSFS-8L	<0.05	<0.05	<0.05	<0.05 -
FSFS-5L	<0.05	<0.05	<0.05	0.70396	0.5515

). When CBA method was used, the phase error was 2 to 5 times higher, and with little difference in computational cost (ca. 7% smaller) (Table 1) in comparison with the $FSFS$ method. We found a phase error similar to previous works using 20-layers with $CBA$ method [8,9]. The statistical test for the FSFS approach showed significant differences between low (below 10 layers) and high (above 13 layers) vertical resolution, considering the free surface elevation residuals. Thus, simulation with 20 to 10 layers has similar (see Table 4) and acceptable results (Table 1 and Table 2).

Our findings indicated that the accuracy of the vertical velocity component is affected by the vertical resolution applied, influencing the representation of the wave phase directly (Figure 3). For the FSFS method, we have noted that the magnitude of vertical velocity components increases or decreases depending on the vertical discretization used, since the 10 layers setup a threshold value (Figure 3). For the CBA approach, the magnitude of vertical velocity components decreases above 16, leading to a less flexible critical vertical resolution. Moreover, we also observed that the variation in horizontal resolution did not affect wave phase representation, which seems to be more directly related to the wave damping (not analyzed in this paper).

3.2. Wave Propagation over a Submerged Bar

A scheme of the experiment of the wave propagation over a submerged bar with an uneven bottom may be seen in (Figure 5) [32]. At the upward slope of the bar, the shoaling wave becomes non-linear due to the generation of bound higher harmonic. At the downward slope, the depth increases rather fast, and these harmonics become free, resulting in an irregular pattern of waves [26]. The numerical reproduction of this pattern has been shown to be very demanding in terms of the accuracy of the computed dispersion frequency [7].

The computational domain has a total length of 30 m, with an initial undisturbed water level of $0.4$ m, which was discretized using a regular grid of 0.025 m resolution. The simulation is carried along 39 s, with a time-step size $\mathsf{\Delta }$t = $0.005$ s. At the left boundary, a sinusoidal wave condition, with period $T=2$ s and amplitude $A=0.01$ m, was imposed to represent the wave generator of the original experiment. At the right outflow boundary, the experimental absorbed beach was computationally represented by a $5m-sponge$ layer with a combination of a sponge layer technique [34] and a Sommerfeld-type radiation boundary condition, applied to minimize wave reflection, given by:

$${ϵ}_i=\left\{\begin{array}{c}\beta {\left(\frac{x_i-x_{io}}{l_i}\right)}^2\left(\frac{z_m-z}{z_m-z_M}\right)u_{i}if{x}_i\ge x_{io}\hfill \\ 0if{x}_i<0\hfill \end{array}\right.,$$

(32)

where ${ϵ}_i$ is the sponge layer coefficient; $x_{i0}$ is the initial point; $l_i$ the total length. This term must be added in the right side of Equations (1) and (2).

We evaluated the capability of the model in simulating the dispersion properties of the flow, comparing simulated and measured free-surface elevation between 33 and 39 seconds of simulation under different vertical discretization scenarios. We used the proposed $FSFS$ approach since the first numerical experiment showed that the $CBA$ method only might reach similar results to the $FSFS$ approach using a higher number of layers (Figure 3). The $FSFS$ approach was tested for 20, 16, 10, and 5 vertical layers comparing the free surface elevation at the six stations with the measurements of free-surface elevation obtained experimentally by Beji and Battjes [32].

The results (Figure 6) indicated that the non-linearity after the upward slope (a), at the beginning and the middle of downward slope (b–c) were well represented by vertical resolutions from 20 to 10 layers, by comparing the patterns of the experimental results of Higher Amplitude Waves (HW) (e.g., 36.5 s station b) and Low Amplitude Waves (LW) (e.g., 36 s station a). For the stations (d), (e), and (f), after the deshoaling process, the simulations using 10, 16, and 20 layers were less accurate in representing the non-linearity pattern. Specifically, the LW in the station (d) (e.g., 35 to 36 s) and the HW in stations (e) and (f) (36 s at station f). The 5-layers scenario showed lower performance in comparison with the other scenarios, with an oscillatory behavior for the HW in stations (a) and (b), substantial phase errors (station c), and low capacity to represent LW and HW amplitudes (stations e and f).

The FSFS results (

Table 5. Statistics metrics between simulated and experimental results for the six stations for each used layer scenario with the FSFS method.

	Station a: x = 13.5 m				Station d: x = 17.3 m
Metrics	20L	16L	10L	5L	20L	16L	10L	5L
RMSE (mm)	2.67	2.31	3.40	5.23	4.12	4.08	5.57	8.05
BIAS (mm)	−0.67	−0.70	−1.26	−2.70	−0.60	−0.44	−1.34	−3.14
Error (%)	26.47	25.37	34.07	53.94	59.29	57.97	79.03	116.84
KGE	−0.95	−1.06	−2.69	−6.92	−4.57	−3.06	−11.30	−27.85
NSE	0.92	0.94	0.87	0.69	0.76	0.76	0.56	0.08
	Station b: x = 14.5 m				Station e: x = 19.0 m
RMSE (mm)	3.68	4.21	4.18	6.74	4.30	4.23	5.69	7.22
BIAS (mm)	−0.80	−0.73	−1.23	−2.58	−0.65	−0.50	−1.39	−3.28
Error (%)	37.73	43.83	44.30	76.27	48.15	47.35	62.91	75.29
KGE	−0.83	−0.68	−1.80	−4.88	−54.45	−41.24	−117.00	−277.78
NSE	0.82	0.77	0.77	0.41	0.72	0.72	0.50	0.20
	Station c: x = 15.7 m				Station f: x = 21.0 m
RMSE (mm)	3.91	3.64	5.45	7.64	4.57	4.59	5.61	7.09
BIAS (mm)	−0.63	−0.48	−1.37	−3.32	−0.65	−0.48	−1.38	−3.25
Error (%)	44.46	40.98	62.63	88.06	56.63	56.12	67.62	76.43
KGE	−1.13	−0.63	−3.58	−10.09	−9.80	−7.09	−22.02	−53.11
NSE	0.79	0.82	0.60	0.21	0.69	0.68	0.53	0.24

) were capable of satisfactorily representing the phase, amplitude, and wave pattern for all stations, using a higher vertical resolution (20 to 10 layers with FSFS method) (NSE from 0.94 to 0.5). Although, the results had higher RMSE and Bias (between 13.7 and −9.9) relative to the mean maximum and minimum amplitudes for the six stations (4.1 and −6.3 mm) and had high Volume Error (between 26% and 116%), always underestimating the free-surface level, and with low performance in KGE parameter due to the increased Volume Error. As we can see in Figure 6, the high Volume Error is expected. Moreover, the Kruskal–Wallis and Nemenyi statistical tests did not indicate a significant difference between the residues series, except in station “A” for 5 layers scenario.

4. Discussion

The statement that "algorithms which define non-hydrostatic pressure at the center of free-surface computational cell need 10 to 20 vertical layers to resolve wave frequency dispersion to an acceptable accuracy" was wildly reported to justify the use or the proposal of new approaches to the free-surface boundary conditions for non-hydrostatic pressure and momentum equation discretization [7,8,9,10,11,24]. In general, this vertical resolution limitation is addressed to a classic finite difference discretization with cell center non-hydrostatic pressure without boundary treatment at the free-surface computational cell (e.g., [23]).

The improvement in non-hydrostatic hydrodynamic solutions is due to the capability to set the value of dynamic pressure close to zero at the free-surface instead of the top layer, for any vertical resolution using a FSFS condition. In a classic boundary approach, as expected, the calculated non-hydrostatic pressure in the top layer increased with the reduction of the number of layers. This occurs because, in this approach, the non-hydrostatic pressure is estimated at the center of the layer, which can be further away from the surface as the number of layers becomes smaller. When using the proposed FSFS condition, the special treatment that occurs in the upper layer minimizes these effects and provides a more physically consistent numerical solution, since this component disappears into real flows.

The results showed that the wave phase representation of short waves in deep water conditions ($\frac{H}{\lambda }>0.5$ m) is related to the contour condition for the non-hydrostatic pressure at the free-surface, but also is influenced by vertical momentum discretization applied. Our findings showed that the proposed boundary condition was capable of improving the model capability in solving wave celerity and wave frequency dispersion, in relation to the classic boundary condition approach. Also, the proposed boundary conditions have a substantial reduction in the computation effort to reach similar results (ca. 146% $FSFS$ 8 layers and $CBA$ 20 layers), and also a little additional computational effort (ca. 7%) when compared to the same benchmark setup, though with a much higher performance (comparing Table 2 and Table 3).

Besides the boundary condition treatment, the phase representation is also related to the used vertical momentum discretization less dependent on the velocities vertical profile. The number of layers may affect the vertical velocity estimative, thus affecting the non-hydrostatic pressure estimation and, hence, the free-surface elevation and velocity fields.

We also identify that there is a vertical resolution threshold where the vertical velocity component is highly affected in deep water conditions. When the number of the layer is lower than a vertical resolution threshold, the wave phase error substantially increases over simulation time, and it becomes higher with a reduction of the number of layers. The vertical resolution threshold is less restrictive using the $FSFS$ approach (10 layers) in comparison with the $CBA$ approach (20 layers). In shallow water conditions, the vertical resolution seems to have a non-significant impact on wave frequency dispersion. The $FSFS$ approach had satisfactory phase representation to 5 layers or more, and also had an adequate representation of nonlinear behavior, similar to previous work [20,23,24,35]. In summary, the $CBA$ approach may obtain similar results to the $FSFS$ approach if a feasible vertical resolution is used. For instance, the 20-layers $CBA$ simulation has similar results to the 8-layers $FSFS$ simulation, although with a substantial gain in computational cost.

When compared with the solutions of previous works [7,8,9,10,11,12], the proposed approach applied at the free-surface has (i) low cost of implementation and (ii) medium performance. The FSFS approach only required a simple local boundary treatment, and it does not change the vertical moment numeric discretization, as the existing methods do. Despite the numerical improvement and simple computational implementation, we observed a limitation imposed by the vertical resolution (number of layers) on the model accuracy, which makes a medium performance model. As the vertical discretization influences directly the accuracy of horizontal and vertical velocities (see Equations (17), (19), (22) and (25)), the proposed solution did not have good performance using a reduced number of layers (less than 10).

Previous works solution had similar satisfactory results with only 2 layers [7,8,9,12], and 4 layers [10,11] for the same benchmarks, with CPU time less than 0.4 s for these scenarios. Although, the proposed method had similar mean CPU time when compared with the same benchmark setup (10-layers simulation), 1.54 s (Intel Xenon 3.7 GHz), 1.83 s (Pentium 4 2.0 GHz) [8] and 2.02 s (i7 2.93 GHz) [11].

Moreover, the analysis of the results clarifies some questions related to the vertical resolution issue for this kind of cell-centered non-hydrostatic model. When using a treatment for non-hydrostatic pressure, the vertical resolution threshold can be more flexible, allowing to use approximately 2 times fewer layers to solve wave dispersion to acceptable accuracy, in a deep water situation. Although this paper indicated a threshold of 10 layers with the FSFS approach, more analysis is needed to establish a local relation between the number of layers and flow characteristics to solve the wave frequency dispersion properly.

5. Conclusions

The treatment of non-hydrostatic pressure boundary condition at the top layer is mandatory to the numerical model to reach satisfactory results compared to other models in the literature (e.g., [9,11,24,36]). The proposed FSFS approach is a low implementation cost method to improve the performance of non-hydrostatic models, which can reach satisfactory results with 2 times fewer layers than the CBA approach. Thus, reducing the mean computational cost of one time-step simulation by ca. 1.7 times, reaching similar results (CBA $\mathsf{\Delta }t=2.62$ s and FSFS $\mathsf{\Delta }t=1.22$ s).

Besides the improvements, the new boundary condition treatment is still limited by the vertical momentum discretization used, leading to a poor performance with low vertical resolution in a deep water situation (less than 10 vertical layers). In shallow water conditions, the vertical resolution seems to have a nonsignificant impact on wave frequency dispersion. Based on the difference between deep and shallow water conditions, more efforts are still needed to establish a local relation between the number of layers and flow characteristics to ensure that the model properly solves the wave frequency dispersion with minimum vertical resolution.