Numerical integration of the contravariant integral form of the Navier–Stokes equations in time-dependent curvilinear coordinate systems for three-dimensional free surface flows

We propose a three-dimensional non-hydrostatic shock-capturing numerical model for the simulation of wave propagation, transformation and breaking, which is based on an original integral formulation of the contravariant Navier–Stokes equations, devoid of Christoffel symbols, in general time-dependent curvilinear coordinates. A coordinate transformation maps the time-varying irregular physical domain that reproduces the complex geometries of coastal regions to a fixed uniform computational one. The advancing of the solution is performed by a second-order accurate strong stability preserving Runge–Kutta fractional-step method in which, at every stage of the method, a predictor velocity field is obtained by the shock-capturing scheme and a corrector velocity field is added to the previous one, to produce a non-hydrostatic divergence-free velocity field and update the water depth. The corrector velocity field is obtained by numerically solving a Poisson equation, expressed in integral contravariant form, by a multigrid technique which uses a four-colour Zebra Gauss–Seidel line-by-line method as smoother. Several test cases are used to verify the dispersion and shock-capturing properties of the proposed model in time-dependent curvilinear grids.

A different approach for the wave motion simulation is based on the numerical integration of the threedimensional Navier-Stokes equations. Some of the most recent models based on this approach use a coordinate transformation in the vertical direction (named sigma coordinate transformation [7]) by which the Cartesian vertical coordinate is expressed as a function of a moving vertical coordinate, σ , which adjusts to the free surface motions. The use of the above strategy makes it possible to obtain numerical models with excellent dispersive properties even when using a low number (less than ten) of vertical layers [8,9]. A recent development in the three-dimensional models which adopt the σ coordinate consists in the adoption of shock-capturing numerical schemes which can simulate the wave breaking in the surf zone [10][11][12].
The study of the fluid motion in three-dimensional form in domains characterised by complex geometries can be carried out by using boundary conforming curvilinear coordinate systems and by expressing the governing equations in contravariant formulation. In the context of the numerical models, the use of the contravariant formulation is not very common, mostly due to the difficulties related to the calculation of the covariant derivatives. The covariant derivatives of contravariant vectors involve the presence of the Christoffel symbols which, as is known [13], are due to the variability of the base vectors. The discretisation of these symbols introduces computational errors associated with the non-uniformities of the calculation grid that can compromise the accuracy of the numerical solution [14,15].
Some authors have realised numerical simulations of the fluid flows in domains characterised by complex geometries by using the contravariant form of the Navier-Stokes equations in a fixed curvilinear coordinate system [16][17][18][19][20]. When dealing with complex geometries which vary in time, some authors have performed a transformation of the irregular, time-varying physical domain into a regular, fixed computational domain by using time-dependent curvilinear coordinate systems. Rosenfeld and Kwak [16] used a contravariant form of the Navier-Stokes equations in moving curvilinear coordinates for the simulation of flows in cavities with variable geometry. Their contravariant formulation has been written in a discrete form, from which the corresponding differential equation is not directly derivable. Carlson et al. [21] derived a tensorial differential formulation of the contravariant Navier-Stokes equations in moving curvilinear coordinates, which is used by [22]. The contravariant differential formulation proposed by [21] does not include all the additional terms that arise when deriving the temporal derivative of a vector tensor in a time-dependent curvilinear coordinate system and, consequently, has no general validity. A complete differential contravariant formulation of the Navier-Stokes equations in time-varying, curvilinear coordinates was obtained by [23] and [24].
These equations include the covariant derivatives of contravariant vectors. Such covariant derivatives imply the presence of the Christoffel symbols which prevents the convective terms of the motion equations from being expressed in conservative form. It is known [19,25] that the numerical methods for the solution of conservation laws in which the convective terms are expressed in non-conservative form, do not guarantee the convergence to the weak solution, i.e., the solution that may contain discontinuities. To obtain a numerical model for the solution of conservation laws which can converge to the weak solution, it is necessary to express the convective terms of the differential motion equations in conservative form or express the motion equations directly in integral form [25].
In this work, in order to realise a three-dimensional numerical model which is able to simulate the discontinuities in the solution related to the wave breaking on domains that reproduce the complex geometries of the coastal regions, we propose an integral contravariant form of the Navier-Stokes equations, devoid of the Christoffel symbols, in a time-dependent curvilinear coordinate system. The proposed integral contravariant form of the momentum equation is obtained by starting from the momentum time derivative of a fluid material volume and from the Leibniz integral rule for a volume which moves with a velocity that is different from the fluid velocity. To avoid the presence of the Christoffel symbols, the proposed integral contravariant momentum equation is solved in the direction identified by a constant parallel vector field. The proposed procedure can be considered the extension of the one shown in [6] to the three-dimensional case and to time-varying curvilinear coordinate systems. In fact, contrarily to [6], which is limited to two-dimensional depth-integrated momentum equations in fixed curvilinear coordinate systems, in the present paper the resulting equation represents the general integral contravariant formulation of the three-dimensional momentum equation in a time-dependent curvilinear coordinate system. Indeed, taking the limit as the volume approaches zero, with simple passages it is easy to obtain the complete differential formulation of the contravariant Navier-Stokes equations in a time-dependent curvilinear coordinate system, which is the same as the one obtained by [24].
The proposed motion equations are numerically solved, on a time-dependent curvilinear coordinate system, by a finite volume shock-capturing scheme, which uses an approximate HLL-type Riemann solver [26]. The term representing the viscous stress tensor is replaced by a turbulence stress tensor term. Analogously to [12], the turbulent stress tensor is represented by the Smagorinsky sub-grid model and is set different from zero only in the numerical simulation of breaking waves. A coordinate transformation maps the time-varying irregular physical domain to a fixed uniform transformed computational one. The advancing in time of the numerical solution is carried out by a second-order accurate strong stability preserving Runge-Kutta (SSPRK) fractionalstep method in which, at every stage of the Runge-Kutta method, a predictor velocity field is obtained by the shock-capturing scheme and a corrector velocity field is added to the previous one, to produce a non-hydrostatic divergence-free velocity field and to update the water depth. The corrector velocity field is obtained by solving a Poisson equation, expressed in integral contravariant form. This equation is solved through a multigrid method which uses a four-colour Zebra Gauss-Seidel line-by-line method as smoother. To speed up the computation, the calculation code is parallelised by means of the Message Passing Interface (MPI).
The paper is organised as follows. In Sect. 2, the procedure is shown, whereby the proposed integral and contravariant formulation of the motion equations in a time-dependent curvilinear coordinate system is obtained. In Sect. 3, the numerical model is described. In Sect. 4, the results are shown and discussed. Conclusions are drawn in Sect. 5.

Governing equations
In this section, we derive an original integral contravariant form of the Navier-Stokes equations, devoid of the Christoffel symbols, in a time-dependent curvilinear coordinate system. The proposed procedure is the extension of the one presented in [6] to the three-dimensional case and to time-depending curvilinear coordinate systems. In fact, in [6] the authors present a contravariant integral form of depth-averaged motion equations (fully nonlinear Boussinesq equations) that does not give a three-dimensional representation of the fluid flow and is limited to fixed curvilinear coordinate systems. In the model presented in [6], the fluid velocity field is represented exclusively by horizontal depth-averaged velocity components; the discretisation of the physical domain occupied by the fluid is obtained by a fixed two-dimensional system of curvilinear coordinates that subdivide the horizontal plane into grid cells bounded by fixed curvilinear lines. With the above approach, the grid cells represent fixed control volumes and, consequently, the contravariant integral form of the depthintegrated motion equations proposed in [6] is obtained by directly using the expression for the mass and momentum time derivative of a fluid material volume which holds when a fixed curvilinear coordinate system is used.
Contrarily to [6], one of the main elements on which the proposed three-dimensional model for free surface flows is based is that, at every instant, the irregular and time-varying physical domain occupied by the fluid is described by moving curvilinear coordinate lines and is divided into computational grid cells that are bounded by surfaces lying on the moving curvilinear coordinate surfaces (i.e., surfaces on which only one curvilinear coordinate is constant). The upper boundary of the computational grid coincides with the free surface and moves rigidly with it, while internal grid nodes move to preserve a given distribution along the water depth. Consequently, the grid cell faces that lie on the free surface move with a velocity equal to the fluid velocity, while the other grid cell faces move with a velocity different from the fluid velocity.
Furthermore, in the proposed model the time-varying computational grid cells are used as control volumes. These two main elements of the proposed model entail the need to write the Navier-Stokes equations in integral form on a control volume whose boundary surfaces move with a velocity different from the fluid velocity.
In this section, firstly, we introduce a time-dependent coordinate transformation from the Cartesian coordinate system to a moving curvilinear coordinate system that allows us to transform the irregular time-varying physical domain into a regular fixed computational domain. Secondly, we introduce the expression for the time rate of change of a quantity (mass and momentum) associated with a material fluid volume, which holds when an arbitrarily moving control volume is used. Finally, we derive a contravariant and integral form of the Navier-Stokes equations, devoid of Christoffel symbols, in the above time-dependent curvilinear coordinate system.
We consider a time-dependent transformation, x i = x i (ξ 1 , ξ 2 , ξ 3 , τ ), t = τ , from the Cartesian coordinate system x 1 , x 2 , x 3 , t to the curvilinear coordinate system, ξ 1 , ξ 2 , ξ 3 , τ and the inverse transformation, (l) = ∂x/∂ξ l be the covariant base vectors and g (l) = ∂ξ l /∂x the contravariant base vectors. The metric tensor and its inverse are defined, respectively, by g lm = g (l) · g (m) and g lm = g (l) · g (m) (l, m = 1, 3). The Jacobian of the transformation is given by √ g = √ |g lm |. The transformation relationships between the components of the generic vector b in the Cartesian coordinate system and its contravariant and covariant components, b l and b l , in the curvilinear coordinate system are given by Hereinafter, the symbol ∂[ ]/∂τ indicates a partial time derivative in which curvilinear coordinates, ξ i , are kept constant. Analogously, in the expressions where there is an explicit time dependence, the τ symbol indicates that such expression is represented in the curvilinear coordinate system.
The velocity vector of the moving curvilinear coordinates, expressed in the Cartesian coordinate system, is The right-hand side of Eq. (2) represents, in the Cartesian coordinate system, the temporal derivative of the position vector (with coordinates x 1 (τ ), x 2 (τ ) and x 3 (τ )) of the points that move rigidly with the curved coordinate lines and are associated with a set of three fixed values, ξ 1 , ξ 2 and ξ 3 (ξ = cost). Let us indicate by v l the l-th contravariant component of v G . To find an expression for v l , let us indicate by A a generic scalar, by ∂ A/∂x the vector gradient of A and by the point "·" the scalar product between vectors defined in the Cartesian coordinate system. It is well known [27] that the relation between the partial derivative with respect to time of the scalar A in the two coordinate systems is By applying Eq. (3) to the generic curvilinear coordinate ξ l (x 1 , Since ∂ξ l ∂τ = 0 (because the coordinate ξ l is not dependent on τ ) and ∂ξ l /∂x = g (l) (by definition), Eq. (4) reads The term ∂x ∂τ on the left-hand side of Eq. (5) represents the velocity vector (in the Cartesian coordinate system) of the point which moves rigidly with the moving coordinate lines and which is identified by a fixed value of ξ . Such velocity vector is the above-mentioned velocity vector of the moving coordinates v G , defined in Eq. (2). By recalling that v l = g (l) ·v G , Eq.
To express the integral formulation of the contravariant motion equations in a time-dependent curvilinear coordinate system, let us start from the contravariant expression of the three-dimensional Leibniz integral rule. Let ρ and u l be, respectively, the density and the lth (l = 1, 3) contravariant component of the fluid velocity vector. Let V (τ ) and A (τ ) be, respectively, the volume and surface area of a fluid material volume, i.e., a time-varying volume which moves with the fluid and always encloses the same fluid particles. Let V 1 (τ ) be a time-varying control volume that at instant τ coincides with the fluid material volume and is bounded by a surface with an area of A 1 (τ ), every point of which moves with a velocity that is different from the fluid velocity.
By applying the three-dimensional Leibniz rule both to the time derivative of the integral of ρ over the fluid material volume V (τ ) and over the arbitrarily moving control volume V 1 (τ ) and by combining the resulting equations, the time rate of change of the mass of the above fluid material volume, d dτ V (τ ) ρdV , takes the form [28] d dτ where n m (m = 1, 3) is the outward unit vector normal to the surface with area of A 1 (τ ) and v m is the mth (m = 1, 3) contravariant component of the velocity vector with which the points belonging to this surface move.
By equating to zero the right-hand side of Eq. (7), the following integral contravariant form of the continuity equation is obtained d dτ From a general point of view, to express the momentum conservation law in integral form, the rate of change of the momentum of a material volume and the total net force must be projected in a physical direction. The direction in space of a given curvilinear coordinate line changes, in contrast with the Cartesian case. Thus, the volume integral of the projection of the momentum equation onto a curvilinear coordinate line has no physical meaning, since it does not represent the volume integral of the projection of the mentioned equation in a physical direction [13].
We identify a physical direction with the direction of a constant and parallel vector field λ. This vector field is represented in the Cartesian coordinate system by constant and uniform components and in the curvilinear coordinate system by constant and space-varying (not uniform) covariant components, λ l . Indeed, since the base vectors of the curvilinear coordinate system vary from point to point, the values of vector components λ l , that are relative to those base vectors, must vary from point to point to represent the same physical direction at every point.
By using Eq. (7), in which the integrand ρ is substituted by ρu l λ l (which represents, in curvilinear coordinates, the scalar product between the momentum per unit volume vector, ρu, and the vector field λ), the projection of the rate of change of the momentum of the above material volume onto the direction of We equate the rate of change of the momentum of the material volume V (τ ) , expressed by the right-hand side of Eq. (9), to the net force in this direction d dτ where f l (l = 1, 3) represents the external body forces per unit mass vector and T lm are the contravariant components of the viscous stress tensor, that can be expressed in the form T lm = −pg lm + 2μS lm , where μ is the dynamic viscosity and S lm are the contravariant components of the rate of strain tensor [13]. As a parallel vector field, we choose the one which is normal to the coordinate line on which the ξ l coordinate is constant at point P 0 ∈ V . We indicate by ξ 1 0 , ξ 2 0 and ξ 3 0 the coordinates of P 0 . The contravariant base vector at point P 0 , indicated by g (l) ξ 1 0 , ξ 2 0 , ξ 3 0 , is, by definition, normal to the coordinate line on which ξ l is constant and is used in this work to identify the constant parallel vector field. Let λ k ξ 1 , ξ 2 , ξ 3 be the covariant component of g (l) ξ 1 0 , ξ 2 0 , ξ 3 0 , given by For the sake of brevity, we indicateg (l) = g (l) ξ 1 0 , ξ 2 0 , ξ 3 0 and g (k) = g (k) ξ 1 , ξ 2 , ξ 3 . By introducing Eq. (11) into Eq. (10) we obtain d dτ Let us introduce a restrictive condition on the control volume V 1 (τ ): in the following, V 1 (τ ) must be considered as the volume of a physical space that is bounded by surfaces lying on the curvilinear coordinate surfaces. In the curvilinear coordinate system, the volume is the corresponding volume in the transformed space, which is defined as V 0 = ξ 1 ξ 2 ξ 3 . Analogously, in the curvilinear coordinate system, the area of a surface of the physical space that lies on the coordinate surface in which ξ α is constant is indicates the corresponding area in the transformed space which is defined as A α 0 = ξ β ξ γ . It must be noted that the volume V 1 (τ ) and the surfaces A α 1 are functions of time, because they are expressed as functions of the base vectors, g (l) , and the Jacobian of the transformation, √ g, whose values change over time as the curvilinear coordinates follow the displacements of the free surface. Conversely, the volume V 0 and the areas A α 0 are not time dependent. By adopting the volume V 1 (τ ) (defined above) as control volume, in the transformed space, the integral Eq. (12) reads where A α+ and A α− indicate the contour surfaces of the volume V 0 on which ξ α is constant and which are located at the larger and the smaller value of ξ α , respectively. Here, the indices α, β and γ are cyclic. By adopting the same control volume, It must be emphasised that Eq. (13) and Eq. (14), in which the Christoffel symbols are absent, represent the general integral form of the contravariant Navier-Stokes equations in a time-dependent curvilinear coordinate system. Indeed, taking the limit as the volume approaches zero, with simple passages it is easy to demonstrate that Eq. (13) and Eq. (14) reduce to the complete differential formulation of the contravariant Navier-Stokes equations in a time-dependent curvilinear coordinate system. The equations obtained are the same as those obtained by [24].
In order to simulate the fully dispersive wave processes, Eq. (13) can be transformed in the following way.
t be the total water depth, where h is the undisturbed water depth and η is the free surface elevation with respect to the undisturbed level. We indicate the acceleration due to gravity by G, and we split the pressure p into a hydrostatic part, ρG η − x 3 , and a dynamic one, q. To accurately represent the bottom and surface geometry and to follow the wave induced free surface evolution, we consider the following time-dependent transformation from the Cartesian system of coordinates, x 1 , x 2 , x 3 , to the curvilinear system of coordinates, ξ 1 , ξ 2 , ξ 3 where the curvilinear coordinates ξ 1 and ξ 2 conform to the boundaries of the physical domain and the vertical coordinate ξ 3 varies over time to adapt to the free surface movements. From Eq. (15), we get the following relation Consequently, from Eq. (6) it can be deduced that the velocity vector of the moving coordinates, expressed in contravariant components, in the curvilinear coordinate system is The proposed coordinate transformation basically maps the irregular, varying domain in the physical space to a regular, fixed domain in the transformed space, where ξ 3 spans from 0 to 1.

Let
√ g 0 = k· g (1) g (2) , where k indicates the vertical unit vector and indicates the vector product.
It is not difficult to verify that in the specific case of this transformation the Jacobian of the transformation can be written in the form √ g = H √ g 0 . This makes it possible to write the integral contravariant momentum equation expressed by Eq. (13) in an original three-dimensional conservative form, in which the conserved variables are given by the cell-averaged product between the water depth H and the three contravariant components of the punctual velocity u l with l = 1, 3. To this end, let us define the following cell-averaged values in the transformed space Since H and √ g 0 do not depend on ξ 3 , the cell-averaged value of the water depth H represents a twodimensional quantity given by the averaged value of the water depth H on the base area of a water column, Since the velocity components u l are dependent on the spatial coordinates ξ 1 , ξ 2 , ξ 3 , the cell-averaged value Hu l represents a three-dimensional quantity that varies along the three spatial coordinates. By adopting the above defined volume V 1 (τ ) as control volume and by using the definition of cell average given in Eq. (18), in the transformed space, the integral Eq. (13) reads where the gradient of the hydrostatic pressure is split into two parts by using η = H − h, the last integral on the right-hand side of Eq. (19) is related to the gradient of the dynamic pressure and the viscous stress tensor term T kα /ρ is omitted (because it is negligible) and is replaced by the turbulent stress tensor, R kα = ν t S kα , in which ν t is the eddy viscosity. Following [12], ν t is estimated by the Smagorinsky sub-grid model and is set different from zero only for the numerical simulations involving wave breaking. As demonstrated by several authors [11,12,29], the inclusion of a turbulent stress term in the numerical model has a negligible effect on the predicted free surface elevation but improves the representation of the vertical turbulent diffusion, thus producing better predictions of vertical velocity gradients. Let us integrate the continuity Eq. (14) over a vertical water column (between the bottom and the free surface) which is bounded by coordinate surfaces. In the transformed space, this integral reads o on which ξ α is constant and which are located at the larger and at the smaller value of ξ α , respectively, and the last two terms on the right-hand side are calculated, respectively, at the free surface (ξ 3 = 1) and at the bottom (ξ 3 = 0). By applying the impenetrability condition, u 3 = 0, to the bottom (where v 3 = 0, because the bottom is fixed), and by applying the kinematic condition, v 3 = u 3 , to the free surface (that moves with the same velocity as the fluid), the last two integrals on the right-hand side of Eq. (20) vanish and Eq. (20) reads Equation (21) represents the governing equation for the free surface movements. Equation (19) and Eq. (21) represent the expressions of the three-dimensional motion equations as a function of the new conserved variables, given by the two-dimensional cell-averaged values of the water depth, H , and by the three-dimensional cell-averaged values, Hu l , in the time-dependent coordinate system, ξ 1 , ξ 2 , ξ 3 , τ .

Numerical model
The numerical integration of Eq. (19) and Eq. (21) is carried out by a finite volume scheme which uses a staggering strategy, proposed by [30] and used by [12], consisting in defining the velocities at the cell centroids and the pressures at the centres of the upper and lower cell faces (Fig. 1).

Time integration
The state of the system is known at the centre of the calculation cells and is defined by the cell-averaged values of the conserved variables, Hu l and H . The solution procedure uses a second-order accurate strong stability preserving Runge-Kutta (SSPRK) fractional-step method [31] in which, at every stage of the Runge-Kutta method, a predictor velocity field is obtained by means of a shock-capturing scheme (which uses an HLL Riemann solver) and a corrector velocity field is added to the predictor one, to produce a non-hydrostatic divergence-free velocity field and update the water depth.
Starting from Hu l (n) and H (n) , related to the generic time level τ (n) , the updated values Hu l (n+1) and H (n+1) , related to the time level τ (n) + τ , are calculated by a two-stage second-order Runge-Kutta method described below in which superscripts (0RK) , (1RK) and (2RK), refer to values obtained at the end of every stage of the Runge-Kutta method, while superscripts (1) and (2), refer to the intermediate values obtained after a predictor-corrector procedure described hereafter.
For each stage s (with s = 1, 2) of the Runge-Kutta method, the predictor fields of Hu l (s) * and H (s) * are calculated by means of the following expressions where Hu l (s RK −1) and H The predictor velocity field u l (s) * (obtained by dividing the conserved variable Hu l (s) * by H (s) * ) does not consider the dynamic pressure and is not divergence free. To obtain a divergence-free velocity field, u l (s) , that considers the dynamic pressure, at each stage s of the Runge-Kutta method, a corrector field, u l (s) c , is added to the predictor field u l (s) * where u l (s) c = ∂Ψ (s) ∂ξ m g ml represents (in contravariant form) the gradient of a scalar field Ψ (s) . The values of Ψ (s) are obtained by solving the following integral contravariant formulation of the well-known Poisson equation that arises by imposing that the divergence of the corrector field is equal in magnitude and opposite in sign to that of the predictor velocity field Equation (25) is discretised by means of a second-order accurate finite volume scheme and is solved by an iterative method. The solution of Eq. (25) provides the scalar field Ψ (s) by which (Eq. 24) the predictor field is corrected both at the cell centroids and the lateral cell faces. The corrected velocities at the lateral cell faces, u l (s) , are then used to obtain the correct value of the water depth at the generic intermediate stage s of the Runge-Kutta method, by the following equation

Ortho and de-orthonormalisation
At each stage s of the Runge-Kutta method, the advancing in time of the water depth and fluid velocity is performed by solving a local Cartesian-based Riemann problem at the cell face centres. To this end, the fluid velocity vector is expressed with respect to a local system of orthonormal base vectors by extending the procedure proposed by [32] and used in [33,34] to the three-dimensional case. By adopting this strategy, Cartesian-based Riemann problems are solved by an HLL solver commonly used in literature [12], to obtain the discontinuity solution at the cell faces. Let us define with B (1) α , B (2) α and B (3) α the orthonormal base vectors defined at the centre of the face on which the generic ξ α is constant. At the cell faces, these local orthonormal base vectors are easily identified starting from the covariant and contravariant base vectors associated with the curvilinear coordinate lines. For example, at the centre of the face on which the ξ 1 coordinate is constant, the contravariant base vector g (1) is aligned, by definition, with the direction normal to the face. Therefore, B (1) α=1 is obtained by dividing g (1) by its magnitude g 11 . On the same face, the covariant base vector g (3) is aligned, by definition, with the direction normal to g (1) ; thus, B (2) α=1 is obtained by dividing g (3) by its magnitude √ g 33 . B (3) α=1 is found by performing the vector product between B (1) α=1 and B (2) α=1 . By generalising this procedure, the set of three orthonormal base vectors B (1) α , B (2) α and B (3) α is given by where α, β and γ are cyclic. The scalar product between the vectors expressed by Eq. (27) and the point values of the velocity vector reconstructed at the cell face centres, indicated by ( v L ) α and ( v R ) α , provides the orthonormal velocity components which are used, as initial values, in the HLL solver The solution of the local Riemann problem, indicated by u RS α , v RS α and w RS α , gives the updated point values of the velocity vector expressed with respect to the orthonormal base vectors. The scalar product between the above velocity vector and the contravariant base vectors, g (1) α , g (2) α and g (3) α , defined at  (s) * , does not consider the dynamic pressure and is not divergence free. This velocity field is corrected by means of the gradient of a scalar field Ψ (s) that has zero curl and a divergence which is equal in magnitude and opposite in sign to the divergence of the predictor field. The scalar field Ψ (s) is calculated by solving an algebraic system obtained by applying the integral contravariant form of the Poisson equation (Eq. 25) at a control volume V i, j,k+ 1 2 (Fig. 2), that is centred at the upper face cell centres (where the pressure is defined). As is seen in Fig. 2, the V i, j,k+ 1 2 volume lateral faces lie on the same vertical plane as the V i, j,k volume lateral faces, while the V i, j,k+ 1 2 volume lower and upper faces lie over the cell centres (i, j, k) and (i, j, k+1), respectively. By applying Eq. (25) to V i, j,k+ 1 2 and using a second-order accurate finite volume scheme, the following expression is obtained ∂Ψ (s) ∂ξ m g m1 The sum of the terms on the right-hand side of Eq. (30) represents the sum of the fluxes by which the integral form of the divergence predictor velocity field (changed in sign) is discretised. It is possible to notice that, in this sum, the V i, j,k+ defined at the grid nodes i, j + 1 2 , k + 1 2 and i, j − 1 2 , k + 1 2 . It is known from literature [35] that the calculation of such velocity components as a simple interpolation from the velocities defined at the cell centroids leads to pressure-velocity decoupling problems, typical of collocated-grid schemes, and arise in the form of spurious oscillations with high frequency, known as checkboard solution [36]. To solve this problem, in literature the Rhie-Chow velocity interpolations, which introduce numerical diffusivity in the solution and smooth the pressure gradient, are generally adopted. In this work, to avoid the pressure-velocity decoupling problems and preserve the ability of the shock-capturing scheme to represent high gradients of water elevation and fluid velocity at the wave breaking, a different strategy is employed. The contravariant velocity components defined at the V i, j,k lateral faces, are computed starting from the values obtained by solving the local Riemann problems. As an example, at the i + 1 2 , j, k + 1 2 grid node, the values of the contravariant velocity component 1 that are obtained by solving the local Riemann problems at the grid nodes i + 1 2 , j, k and i + 1 2 , j, k + 1 are used An analogous expression is adopted for the contravariant velocity component 2 defined at i, j This strategy makes it possible to avoid the pressure-velocity decoupling problems and take advantage of the shock-capturing abilities of the proposed scheme to represent steep fronts at the wave breaking. The calculated contravariant predictor velocity components (Eqs. 31 and 32) are corrected by means of the gradient of the Ψ (s) potential, calculated in the same grid nodes. Indeed, the ∂Ψ (s) ∂ξ m g m1 By discretising the corrector velocities, ∂Ψ (s) ∂ξ m g m1 where the R i, j,k+ 1 2 term on the right-hand side represents the divergence (changed in sign) of the predictor velocity field (right-hand side of Eq. 30) and the 19 a l coefficients are shown in "Appendix". The matrix associated with the algebraic set of equations expressed by Eq. (34) has 19 non-null diagonal coefficients. The system solution is obtained by an iterative multigrid method [37] that uses a combination of the four-colour Zebra Gauss-Seidel line-by-line as a smoother [16].

Boundary conditions
The upper boundary of the computational domain on which the momentum equation expressed by Eq. (19) is integrated corresponds to the water free surface. The above boundary corresponds to the upper face of the V i, j,n3 computational cells, where n3 indicates the number of cells in the vertical direction. The curvilinear coordinate transformation defined by Eq. (15) allows the grid nodes identified by i, j, n3 + 1 2 (which lie on the upper surface of cell V i, j,n3 ) to represent, at each instant, the water free surface (on which ξ 3 is constant). Moreover, it is to be noticed that, by definition, the contravariant component 3 of the velocity vector is orthogonal to the coordinate surface on which ξ 3 is constant. Consequently, in our work, the application of the kinematic condition at the free surface is achieved by imposing that the difference between the contravariant component 3 of the fluid and grid velocity at the cell nodes i, j, n3 + 1 2 is null On the i, j, n3 + 1 2 grid nodes, zero-value condition of the Ψ potential is imposed The lower boundary of the domain coincides with the cell centroids of the first calculation cell along the vertical direction, V i, j,1 . At these nodes, the no-slip condition can be assigned by imposing null velocity contravariant components at the i, j, 1 cell centroids. The slip condition can also be assigned by imposing a null normal derivative for the fluid velocity contravariant components 1 and 2, and a null value for the contravariant component 3.
The lower boundary of the computational domain on which the Poisson equation is solved corresponds to the lower face of V i, j,1 . Hence, such lower boundary is placed on the i, j, 1 2 grid nodes, i.e., half-cell down the i, j, 1 nodes of the impermeable bottom of the domain. In such a way the nodes corresponding to the bottom are placed between the nodes at which the Ψ (s) potential is defined. The lower boundary condition of the Poisson equation is obtained by imposing a zero-value condition of the Ψ (s) potential derivative in the direction normal to the bottom. This condition reads Such condition is equivalent to imposing the contravariant corrector velocity component 3, u 3(s) c , (normal to the bottom) to be null at the bottom. Consequently, the lower boundary conditions for the bottom-normal velocities are those assigned in the predictor field calculation At the lateral boundaries of the domain, a ghost cell is used. On the closed boundaries, a null derivative in the normal direction is assigned for the free surface elevation and a null flux is imposed through the boundary. On the open boundaries, a null derivative in the normal direction is imposed for the free surface elevation and velocity components. At the open boundary on which the wave motion is generated, the free surface elevation and velocity components, obtained by the Airy theory, are assigned at the ghost cell centroids. The

Results
The proposed three-dimensional non-hydrostatic shock-capturing numerical model is validated by reproducing experimental test cases on time-dependent curvilinear grids. The ability of the model to simulate wave propagation over constant and varying depths and the nonlinear hydrodynamic wave phenomena related to it, is verified.

Standing wave in closed basin
In this subsection, the dispersion properties of the proposed three-dimensional non-hydrostatic shock-capturing numerical model are verified by reproducing a test case proposed by [38] and used by several authors [3,4,12,39]. The main objective of simulating this test case is to verify the accuracy of the proposed model when modelling linear dispersive waves [39].
The test case simulates a standing wave (of amplitude a = 0.1 m and wave number k = π/H ) in a closed basin with length L = 20 m and depth H = 10 m. At instant t = 0, the water velocity is set to zero everywhere and the initial free surface elevation is given by the Airy linear theory, η = acos(kx). Starting from this initial condition, the free surface elevation is left free to oscillate. For this standing wave, the dimensionless parameter ε = a/H which measures the nonlinearity is equal to 0.01. This value results much smaller than 1, and, consequently, the considered wave can be analytically represented by using the Airy linear theory. The turbulent stress terms are set to zero, and a slip condition is imposed at the bottom (null horizontal velocities gradients). In this way, the numerical solution is characterised by such small values of the water velocity components and free surface elevation gradient as to make the nonlinear effects in the wave celerity negligible. Thus, the obtained results are comparable with the Airy wave theory analytical solution. Such analytical solution consists of a standing wave whose length is equal to that of the basin and wave period T = 3.5857 s. The numerical test is run for 30 s by using a time step of t = 0.002 s and a horizontal grid spacing of x = y = 0.1 m. Figure 3a     To verify the dispersive properties of the proposed model with respect to the number of vertical layers and the dispersive wave parameter, the test case has been repeated by increasing the number of vertical layers. For the three different wave cases, the depth of the basin is discretised with 3, 5, 8 and 10 uniform layers. The numerical errors are calculated by using the normalised root mean square expression proposed by [12] numerical error = 1 H in which H is the wave height at the considered abscissa x, N is the number of the compared data, η i is the value obtained by the numerical simulation, and η α is the value of the analytical solution. Figure 6 shows the numerical errors at x = 17.5 m as a function of the number of vertical layers and wave dispersion parameter. It can be noticed that the numerical errors decrease as the number of vertical layers increases. These results show that the proposed model has good dispersion properties. An accurate solution (error less than 2 %) is obtained even when using only three vertical layers (unlike [38,39] that employed 20 vertical layers). Such result is in accordance with the use of a Keller-box scheme (in which the pressure is located at the cell faces) that allows the computation of the non-hydrostatic pressure also in the surface cells [30] and with the results obtained by [12,30]. The proposed model guarantees good dispersion properties with few vertical layers, thus containing the computational costs.

Wave train propagating over a bar
In this subsection, the proposed model performances are tested by numerically reproducing the experiment carried out by [40,41] and used by several authors [3,8,10,12,38,39,42,43]. In this experiment, steep but nonbreaking waves (with period T = 2.02 s and height H = 0.01 m) propagate in a wave flume (of length L = 37.7 m and initial water depth equal to 0.4 m) with a submerged trapezoidal bar consisting of an offshore slope of 1/20, a 2 m horizontal crest (upon which the water depth is reduced to 0.1 m) and an onshore slope of 1/10. At the end of the flume, a plane beach with a 1/25 slope is present. The isolated submerged obstacle poses a challenging case due to the phenomenon of wave decomposition taking place in the deepening region behind the obstacle. In fact, it has been found that the shoaling would occur when the wave train passes over the upward slope [43]. Due to the interaction with the bar, the shoaling wave becomes nonlinear through the generation of bound higher harmonics which travel phase-locked to the primary wave [3,43]. This wave decomposition past the obstacle gives rise to rapidly varying wave forms, which can be simulated accurately only by a sufficiently nonlinear model with good dispersion characteristics [41].
To measure the free surface displacements, seven gauges were placed along the wave flume: the first one (G1) and the last one (G7) were placed in the constant-depth sections, respectively, before and after the submerged trapezoid; gauge 2 (G2) was placed at the offshore slope of the trapezoid, while gauge 5 (G5) and gauge 6 (G6) were placed at the onshore slope of the trapezoid; gauge 3 (G3) and gauge 4 (G4) were placed on the horizontal crest of the trapezoid. Bottom topography and gauge positions are sketched in Fig. 7. Figure 8 shows the free surface elevation on the submerged trapezoidal bar bottom topography and contours of p at t = 40 s. Figure 9 shows the spatial evolution of steep regular waves propagating over the submerged trapezoidal bar by means of a time domain comparison between the numerical and experimental surface displacement at the location gauges 2 through 7.
At the first two gauge locations, the wave train remains almost sinusoidal with nearly perfect agreement between the numerical results and the experimental data (Fig. 9a). From gauge 2, where the wave starts to climb the slope, the gradually increasing wave steepening due to shoaling effects is observed. Here, the   [41] shoaling produced by the bar increases the wave height and reduces the wave length, increasing in such a way the ratio H/L and therefore the nonlinear effects. At gauge 3 and gauge 4 (Fig. 9c, d), where the wave rides over the top of the breakwater, the growth of a secondary wave becomes apparent, resulting in a profile distortion. Behind the breakwater, the secondary wave mode gains energy from the main wave mode and the wave length becomes shorter. This can be seen at the last two gauge locations (Fig. 9e, f). The numerical model simulates the evolution and decomposition of the wave field very well, with the initial steepening on the seaward slope and the enhancement of the wave profile distortion. A good agreement of the numerical results with the experimental data can be observed in all the gauges considered. The shoaling process on the offshore slope of the bar is well described by the proposed model by using only three vertical layers, while on the onshore slope some discrepancies between the numerical results and the experimental data are to be noticed, where the wave deformation is rapid. The proposed model can reproduce the nonlinear wave transformations over and behind a submerged bar, thus well reproducing the wave-structure interaction.

Wave train propagating on a varying depth
Wave breaking on a gently sloping plane beach is a key point for any testing of a breaking model [44]. In this subsection, the ability of the proposed model to represent the nonlinear wave phenomena (shoaling and breaking) is verified by reproducing the laboratory test for breaking waves on a sloping beach performed by [45] and used by many authors [46]. The experiment simulates an incoming wave train of height H 0 = 0.156 m and period T = 1.79 s which propagates in a 55 m long wave flume characterised by an initial constant depth of 0.85 m, followed by a plane sloping beach of 1:40 (see Fig. 10). The breaking behaviour at the breaking point of this test is characterised by a surf similarity parameter equal to ξ 0 = tan α/ √ H 0 /L 0 = 0.125 (where α is the bed slope angle and the ratio H 0 /L 0 = 0.032) which falls into the category of "spilling breaking" [45].
To demonstrate that the results obtained with the proposed model are not influenced by the grid distortion, the simulation is performed by using a Cartesian computational grid and a highly distorted curvilinear computational grid. With respect to the highly distorted curvilinear computational grid, Fig. 11a, b shows the grid distortion and a detailed area of the computational domain.
The results obtained with the two different grids are compared ( Table 1) by means of the root mean square error, σ v y , of the difference between the values numerically obtained and the null values of the v y velocity  Fig. 12 Wave train propagating on a varying depth. Instantaneous wave field obtained with the curvilinear highly distorted computational grid component (because the propagation direction is parallel to the x-axis). The σ v y error calculated for the highly distorted grid case differs by less than 1 % from the corresponding Cartesian grid error. Figure 12 shows an instantaneous wave field obtained with the curvilinear highly distorted computational grid. Figure 13 shows a three-dimensional and plane view detail of such instantaneous field in the area where the grid distortion (solid lines) is maximum. From these figures, it can be seen that, despite the grid distortion, as the wave propagates from deep water up to the shoreline, the wave train maintains even wave fronts and does not show spurious oscillations. From such results, it can be deduced that the grid distortion does not affect the ability of the proposed model to simulate the nonlinear wave phenomena. Moreover, these figures show the ability of the model to represent the wave shoaling and the steep front of the breaking wave without significant disturbances in the solution. Figure 14 shows the detail of the instantaneous water depth, eddy viscosity and fluid velocity along the central cross section, obtained with the proposed model. In this Figure, it is possible to see the steepening of the wave front, as the waves approach the shore. This occurs up to the breaking point, after which the numerical solution can be considered as a moving discontinuity in the flow variables and is characterised by a sharp reduction in the wave height. From this figure, the wave breaking turbulent production can be deduced by the contours of the eddy viscosity: as the wave train propagates towards the shoreline, the turbulence production reaches a peak at the breaking point and gradually drops after breaking.
One of the main concerns is the validation of the shock-capturing features of the numerical model, and for this reason the attention is focused on the ability of the model to catch the breaking position and to predict wave height decay after the breaking point. Figure 15a, b shows, respectively, the time averaged wave height and the mean water level comparisons between the numerical results obtained with the highly distorted curvilinear computational grid and the experimental data. In the deepwater region as well as in the wave breaking region, the overall agreement in terms of wave height and water level is good. In particular, by observing Fig. 15a, b it can be seen, respectively, that the position of the breaking point and the wave set-up are well predicted by the numerical model. Table 1 shows the root mean square error, σ wh , of the difference between the time averaged numerically computed wave height and the corresponding experimental data, and the root mean square error, σ mwl , of the difference between the numerically computed mean water level and the experimental data. From Table 1 it can be deduced that, even for these two quantities, the numerical error is not influenced by the grid distortion.

Rip current test
In this subsection, the laboratory experiment carried out by [47] has been reproduced. This test aims to verify the ability of the proposed model to reproduce wave propagation from deep water up to the shore, the wave breaking in a spatially varying bottom and the induced nearshore current. The laboratory experiments were carried out in a 30 × 30 m basin characterised by a plane sloping beach of 1:30 with a rip channel excavated along the centreline. The bottom variation is given by cos 10 π (15−y) 30 7 < x < 25 x ≤ 25 Figure 16a shows the bottom variation of the calculation domain. Because of the presence of an axis of symmetry perpendicular to the wave propagation direction, only half of the experimental domain has been reproduced. At the lateral boundary that coincides with the symmetry axis of the experimental domain, reflective boundary conditions are imposed in order to simulate the effects of the specular half portion of the experimental domain. Figure 16b, c shows the beach profile at two significant cross sections, one inside the rip channel (y S1 = 14.9625 m) and one at the plane beach (y S2 = 1.9875 m), where the experimental data reported by [47] are available. The alongshore non-uniformity of the bottom causes the wave shoaling and breaking to occur more offshore at the plane sloping beach than at the rip channel. It produces a significant gradient of the water surface elevation in the alongshore direction. This gradient drives the flow from the plane sloping beach towards the channel. For the mass conservation, at the centreline of the channel the water flowing from the plane sloping beach is then driven seaward. It gives rise to a nearshore circulation, which is expected to be numerically reproduced by the proposed model.
The experiments considered several different incident wave conditions. Here, the monochromatic, regular, incident waves are considered with a period of T = 1.25 s and wave height of H = 0.07 m. The simulation is carried out on a curvilinear computational grid (Fig. 17) and is run for 520 s with a time step of 0.0025 s. Figure 18 shows a three-dimensional and a plane view detail of an instantaneous wave field (contour) carried out with the curvilinear computational grid (solid lines), at the time when the breaking-induced circulation is fully developed. From this figure, it is possible to notice that, when approaching the beach, the wave fronts rotate and the wave field shows an increase in the wave height in correspondence with the channel location due to the occurrence of a pronounced rip current along the channel. Figure 19a, b shows the instantaneous water depth and velocity field along the cross sections S1 and S2 at the time when the breaking-induced circulation is fully developed. In both figures, it is possible to notice the steepening of the wave front and the sharp reduction in the wave height subsequent to the wave breaking. In section S1, due to the rip channel, incipient breaking occurs closer to the shore than at section S2. Moreover, from the comparison between the two figures it is possible to see that at the centre of the rip channel, before breaking, the wave height reaches a maximum that is significantly higher than in the plane sloping beach. This effect is due to the presence of a seaward current (rip current) that produces a further increase in the water surface elevation.
The occurrence of the rip current is due to the fact that the differences in wave breaking position (shown in Fig. 19) induce different water surface elevation in the alongshore direction. In fact, such gradient of the water level drives an alongshore current that turns offshore producing the rip current at the rip channel position, as can be seen in Fig. 20. This figure shows a plane view detail of the time averaged velocity field near the bottom, from which it is possible to focus on the circulation which has developed. It should be noted that vectors are shown only at every four grid points.
The mean current velocity variation along the rip channel is derived from the time averaged velocity field near the bottom. In Fig. 21, such mean current velocity variation is compared with the experimental data obtained by [47] for unidirectional and multidirectional random waves with the same characteristics as the numerically reproduced monochromatic wave. It can be noticed that the numerical results are in good agreement with the experimental data. The model makes it possible to capture both the general trend of the mean current velocity, and its peak value.
In Fig. 22, the wave heights computed with the proposed model are compared with the wave heights measured by [47] along the two above-mentioned cross sections and reported by [48]. It can be noticed that Fig. 16 Rip current test. Bathymetry and beach profiles S1 and S2 located, respectively, along the rip channel (a) and a plane beach section (b) Fig. 17 Rip current test. Curvilinear grid of the computational domain. Only one out of three coordinate lines is shown the numerical results in terms of wave height are in good agreement with the laboratory measurement for the incised-channel on a plane beach. The shock-capturing ability of the model is proved as the breaking point is well predicted both in the rip channel and the plane beach.

Conclusions
A three-dimensional non-hydrostatic shock-capturing numerical model for the simulation of wave propagation, transformation and breaking has been presented. The proposed model is based on an original integral formulation of the contravariant Navier-Stokes equations, devoid of Christoffel symbols, in general timedependent curvilinear coordinates. The contravariant Navier-Stokes equations are numerically solved, on a time-dependent curvilinear coordinate system, by a finite volume shock-capturing scheme, which uses an approximate HLL-type Riemann solver. A coordinate transformation maps the time-varying irregular physical domain to a fixed uniform transformed computational one. The advancing in time of the numerical solution is carried out by a second-order accurate strong stability preserving Runge-Kutta (SSPRK) fractional-step method in which, at every stage of the Runge-Kutta method, a predictor velocity field is obtained by the Fig. 18 Rip current test. Three-dimensional and plane view detail of an instantaneous wave field (contour) at the time when the breaking-induced circulation is fully developed and the curvilinear computational grid (solid lines) Fig. 19 Rip current test cross sections. Instantaneous water depth (contours) and velocities (arrows) along the rip channel (a) and the plane beach (b), obtained with by the proposed model at the time when the breaking-induced circulation is fully developed shock-capturing scheme and a corrector velocity field is added to the previous one, to produce a non-hydrostatic divergence-free velocity field and to update the water depth. The corrector velocity field is obtained by solving a Poisson equation, expressed in integral contravariant form. This equation is solved through a multigrid method which uses a four-colour Zebra Gauss-Seidel line-by-line method as smoother. The proposed model has been validated by reproducing experimental test cases on time-dependent curvilinear grids. It has been shown that the proposed model has good dispersive and shock-capturing properties in time-dependent curvilinear grids.   [47] for unidirectional (times) and multidirectional (square) random waves Fig. 22 Rip current test. Wave height along the rip channel (a) and along a plane beach section (b). Comparison between the numerical results (solid line) obtained with the proposed model and the experimental data (symbols) from [47]  Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.