Computation of free surface waves in coastal waters with SWASH on unstructured grids

This paper aims to present the extension of the non-hydrostatic wave-ﬂow model SWASH with the co- volume method to build discretization schemes on unstructured triangular grids. Central to this method that is free of spurious pressure modes, is the use of dual pairs of meshes that are mutually orthogo- nal, such as the Delaunay-Voronoi mesh systems. The approximants sought are the components of the ﬂow velocity vector normal to the cell faces of the primal mesh. In addition to the covolume approach, a novel upwind difference scheme for the horizontal advection terms in the momentum equation is proposed. This scheme obeys the Rankine-Hugoniot jump relations and prevents the odd-even decoupling of the velocity ﬁeld accordingly. Moreover, cases with ﬂow discontinuities, such as steady bores and broken waves, are properly treated. In spite of the low-order accuracy of the proposed method, unstructured meshes easily allow for local reﬁnement in a way that retains the desired accuracy. The unstructured-grid version of SWASH is applicable to a wide range of 2DH wave-ﬂow problems to investigate the nonlinear dynamics of free surface waves over varying bathymetries. Its eﬃciency and robustness is tested on a number of these problems employing unstructured triangular meshes.


Introduction
A commonly encountered coastal engineering application involves the simulation of dispersive waves in coastal waters. Waves are usually present almost everywhere in the domain of interest, while the wave dynamics can be dominant in a small part of the entire domain, e.g. the surf zone due to wave breaking and generating currents. For many large-scale wave problems, the mesh only needs to be concentrated in particular areas with strong wave features and large current gradients. Local mesh refinement is therefore a suitable approach to offer sufficient grid resolution in otherwise under-resolved areas (e.g. breaking zones, swash zones, wetlands, coastal flood defences).
The development to increase the flexibility of structured grid methods to allow variable grid resolution and local mesh refinement has a long history. Traditional ocean and coastal models usually employ curvilinear grids with domain decomposition techniques. Such grids are typically orthogonal so that the coordinate transformation offers advantages in solution algorithm efficiency, whereas undesirable numerical artefacts arising from grid nonorthogonality, high cell aspect ratio and skewness are avoided. E-mail address: m.zijlema@tudelft.nl Structured grids greatly facilitate the use of high-order discretization methods. Such methods provide a higher spatial accuracy and hence reduce the number of grid cells required for a given level of solution accuracy. Nevertheless, we argue that loworder discretizations on locally fine grids might be more efficient. Low-order upwind-biased schemes often exhibit desirable stability properties but may suffer from numerical damping, especially at relatively coarse grids, which can easily be counteracted by grid refinements. Furthermore, they only require small discretization stencils and often deliver a good scaling capability for high performance computing (HPC). In contrast, high-order discretization methods usually lack robustness as they have the potential to destroy the stability of the underlying numerical approach. In particular, central and high-order upwind discretizations tend to display non-physical oscillations and commonly result in instabilities.
On the other hand, when local mesh refinement is applied to flow problems in water of finite depth, the accuracy is often limited by the error in topographic and forcing data rather than the discretization error [1] . Hence, grid refinement substantiates the benefit of low-order discretizations. Despite the fact that low-order (upwind) methods are typically deprecated in the literature, we believe that such methods are considered adequate for predicting wave processes including propagation, shoaling, refraction and wave breaking, provided that sufficiently fine meshes are employed. While such processes are governed by mass and momentum balances, this may not be true for the transport of contaminants in surface waters and variable density solute transport (e.g. salinity, temperature, sediment). That is to say, high resolution schemes are desired to avoid a detrimental loss of tracer concentrations [2] .
A drawback of structured meshes is the inherent difficulties associated with the generation of appropriate boundary-fitted grids and excessive grid refinement. The use of unstructured grids offers a remedy to these problems as they are more flexible and easily provide highly variable grid resolution in complex geometries. In this regard, the domain is commonly discretized into triangular or tetrahedral elements with no implied connectivity. Unstructured meshes have received great attention for the past several years in the ocean and coastal modelling community. A general review on this topic can be found in [3] . One principle feature of an unstructured grid is that advanced grid refinement techniques are easier to implement. Unstructured grid methods are therefore particularly suited for wave problems with complex geometries that require a reasonably non-uniform grid resolution. However, they do not allow for ease of implementation of high-order discretizations. Also, high-resolution schemes on unstructured grids usually do not take the full advantage of higher order accuracy that can easily be achieved on structured grids [2] . In fact, low-order methods for unstructured meshes are still in use today in established codes [1] . In spite of this lower order accuracy, we premise that a higher rate of convergence is achieved by local grid refinement rather than use of high-order (central) discretizations.
A common feature of the traditional numerical models employing unstructured meshes is the application of the finite volume method, which is a widely used approach in the CFD community. This method directly utilizes the integral formulation of conservation laws in divergence form which are subsequently discretized by means of flux approximations applied on non-overlapping control volumes in either primal or dual grids. Standard finite volume schemes thus have the benefit to discretely conserve unknowns of PDEs on arbitrary meshes. In the context of the solution of shallow water problems, such schemes typically involve the vertex-centred or cell-centred discretization of the conserved variables, i.e. the water depth and the depth-integrated velocity vector, on colocated and semi-staggered grids.
Yet, low-order vertex-centred finite volume schemes are known to be rather inaccurate or even inconsistent on irregular grids (for details see, e.g. [4][5][6] ). Such schemes require polygon-shaped duals that typically provoke highly distorted volumes. For instance, a detailed analysis of Svärd et al. [4] reveal that vertex-centred finite volume approximations yield order of accuracy less than one for non-smooth mixed grids of triangles and quadrilaterals. In particular, the non-uniformity in the size of dual mesh polygons results in a loss in discretization accuracy of one order. The cause is the differences in fluxes on opposite sides of the control volume that do not cancel in pairs when the mesh is non-uniform. It is likely that their capacity to achieve desired physical accuracy on irregular grids is rather restricted. On the other hand, the solution quality of a cell-centred finite volume scheme, operating on primal cells, can easily be greater than that of a vertex-centred scheme. However, cell-centred schemes incur larger overheads than vertexcentred schemes as they involve approximately twice as many unknowns. Therefore, we argue that low-order finite volume methods using unstructured meshes is less tractable for robust and efficient large-scale simulations of wave dynamics requiring high accuracy.
Another class of numerical methods involves the use of staggered grids. In such grids the velocity components are staggered with respect to the pressure. The pressure is placed in the centre of the grid cells, while the velocities are stored at the cell faces. Contrary to colocated and semi-staggered grid schemes, staggered grid methods prevent the well-known odd-even decoupling between velocity and pressure which is typically the key ingredient to the development of robust and efficient incompressible flow models. Staggered finite difference methods were first developed by Harlow and Welch [7] for incompressible flows on Cartesian grids 1 and has a long history of successful achievements in coastal and ocean modelling (see, e.g. [9][10][11] ).
A fruitful extension of the classical staggered grid approach to unstructured triangular and tetrahedral meshes is the covolume method [12,13] . This method is designed as a low-order method that does not encounter spurious pressure modes. The covolume method makes use of two dual orthogonal grids to approximate the flow velocity and pressure. In this respect, the Delaunay triangulation (primal grid) and the corresponding Voronoi diagram (dual grid) are the natural choice. Only the normal velocity components are defined on the faces of the primal mesh, whereas the pressure is stored at the vertices of the dual grid, i.e. the circumcentres of the primal cells. The accuracy of the covolume scheme is generally first order in space but can be second order accurate for a regular mesh [14] . However, more importantly, covolume techniques usually respect underlying physical principles due to a range of fundamental properties, such as topology, conservation and symmetry, leading to the development of highfidelity discretization methods for incompressible Navier-Stokes and Maxwell's equations [15,16] . The effect of discretization error originated from such discretization approximations is thus limited, in that the numerical solution is merely influenced by the mesh resolution and mesh quality.
Over the past two decades, the covolume method, frequently referred to as the unstructured C-grid discretization, has been successfully employed for the modelling of large-scale ocean and coastal flows. See, for instance, Perot and Nallapati [17] , Stuhne and Peltier [18] , Kramer and Stelling [19] and Kleptsova et al. [20] . Additionally, regional coastal and estuarine ocean models like UnTRIM [21] , SUNTANS [22] and D-Flow FM [23] are widely used by many researchers. Both UnTRIM and SUNTANS can also be applied to include non-hydrostatic dynamics for environmental flow problems (e.g. internal waves and oceanic fronts).
SWASH is a non-hydrostatic wave-flow model developed at Delft University and exploits the finite difference method on staggered Cartesian grids [24] . This phase-resolving model has been subjected to extensive benchmarking and has been successfully applied to a numerous wave and flow applications. Recent examples can be found in [25][26][27][28] . The present work describes the extension of SWASH to unstructured triangular grids. To this purpose, we adopt the covolume strategy because of the wish to maximize the physical accuracy, stability and efficiency of the numerical wave simulation. The focus of the present study is on tailor-made discretizations suitable for the simulation of wave dynamics. They ensure the correct handling of discontinuities and shocks as manifested in bores, hydraulic jumps and broken waves.
In what follows, the shallow water equations including nonhydrostatic pressure are introduced in Section 2 . In order to facilitate transparency on discretization schemes, the governing equations are presented in a depth-averaged form. Next, the concept of the covolume approach and its numerical implementation are treated in Section 3 . Special attention is paid to the discretization of momentum advection based on a genuinely upwind-biased approach that guarantees conservation of momentum flux locally.
The key contribution of the present study is to demonstrate the efficiency and robustness of the proposed method for wave-related simulations. In Section 4 a number of test cases are presented for validation purposes and the results are discussed. Conclusions are provided in Section 5 .

Governing equations
The governing two-dimensional, primitive variable equations for the depth-averaged, non-hydrostatic, free surface flow of an incompressible fluid over a bed topography d ( x, y ) are given by ∂ζ ∂t with the coordinate directions x, y and z aligning in the east, north, and vertical directions, respectively. The bed level d ( x, y ) is measured from the reference level z = 0 (positive downwards), whereas ζ ( x, y, t ) is the surface elevation with respect to the reference level (positive upwards). The water depth is given by is the flow velocity with the depth-averaged components u ( x, y, t ) and v ( x, y, t ) along the x and y coordinates, respectively, q = h u is the mass flux, p ( x, y, z, t ) is the non-hydrostatic pressure (normalised by the density), g is the gravitational acceleration, and c f is the dimensionless bottom friction coefficient. The shallow water equations (1) and (2) are derived from integrating the mass conservation and the momentum balance over the depth, respectively, whereas the total pressure is decomposed into its hydrostatic and non-hydrostatic components (see, e.g. [21,22,24] ). In the case of a hydrostatic pressure distribution, i.e. p ≡ 0, these equations can be reformulated as a set of nonlinear hyperbolic equations, and may thus generate discontinuous solutions featuring shock waves [29] . Such solutions can readily be understood as weak solutions in the variational context. Using Leibniz' rule, a conservative expression for the gradient of non-hydrostatic pressure is obtained [30] ζ −d with p b the non-hydrostatic pressure at the bed. This pressure is associated with the vertical motion that is governed by the following equation where w s is the velocity in the z−direction at the free surface. This equation is derived using the Keller-box method to further improve the dispersive behaviour of the waves [24,30] . The vertical velocity at the bed w b can be found by means of the following kinematic condition Finally, the system of equations is complete with the following equation which ensures conservation of local mass. The momentum equation (2) is symbolically written in vector form which is not the most natural form to describe the underlying physics in wave-dominated regimes. Instead, we consider the following momentum balance characterizing the flow below a free surface in water that moves in a nominal downstream direction defined by the unit vector n f , with u f = u · n f the local streamwise, depth-averaged flow velocity. Vector n f is assumed to be slowly varied such that pressure gradients redirect momentum through bends whereas the flow velocity u f is aligned with the direction of the pressure forces. Furthermore, the lateral variability in the free surface slope is negligible on the scale of flow induced motion, especially in coastal regions [31] . The depth-integrated momentum equation (4) essentially describes a local balance between flow acceleration, frictional deceleration and the pressure force acting on an element of fluid along streamlines, and governs momentum per unit volume normalised by the density, i.e. u f . The mass flux q is the transport velocity of momentum and the second term of Eq. (4) represents the advective acceleration of the flow field with a nonlinear spatial effect on u f . For instance, it causes the change in the momentum in an open channel due to expansion or contraction, i.e. slowing down or speeding up the fluid locally. Such flow transitions are typically rapid and frequently arise in hydraulic jumps, dam breaks, tidal bores and flows over a weir.

Discretization on unstructured triangular meshes
This paper is not intended to provide an outline of the overall numerical methodology of SWASH, but rather to explain the essential steps of the discretization of the governing equations on unstructured staggered grids since this topic differs from the structured-grid version of SWASH as described in Zijlema et al. [24] . Therefore, the discussion on temporal discretizations and the solution procedure is not included in the present study as they have been previously treated in depth (in, e.g. [24] ).
The starting point is that the flow velocities are segregated from the surface elevation and non-hydrostatic pressure and are assigned to different grid locations based on a staggered grid arrangement. This is further discussed in Section 3.1 . Note also that the solution of non-hydrostatic pressure is not pursued herein. The interested reader is referred to the aforementioned reference for detailed description of the pressure correction technique aiming to enforce local mass conservation. Section 3.2 outlines the concept of the covolume scheme for spatial discretization. The implementation of this method in SWASH is elaborated in detail in Section 3.3 . Next, the discretization of the advective acceleration is dealt with in Section 3.4 . Finally, a simple wet-dry scheme is treated briefly in Section 3.5 .

Staggered grid
Staggered grid methods typically employ velocity components as the primary discrete unknowns in the spatial discretization. In this respect, the flow velocity component u f is a proper one. This physical variable is located at the face midpoint in the direction n f normal to the cell face, whereas the physical scalars (e.g. water depth, pressure) are positioned at the cell circumcentre. As a result, the flow velocity normal through the face of the grid cell is influenced directly by the water depth and non-hydrostatic pressure at the circumcenters on either side, because it is the pressure gradient that drives the flow. In addition, the change in volume cell, confined by the surface elevation, is due to the mass fluxes hu f at the faces of the grid cell. This physical duality expresses a natural property of the staggered mesh approach.
Staggered schemes thus require the existence of the cell circumcentres since the line segment connecting the circumcentres of the two cells intersects with their shared face, while they are orthogonal to each other. The most common polygons with a unique circumcenter are the triangles and rectangles. Furthermore, any two non-colinear sides of such cells span a basis in R 2 . (Note that the third side of a triangle is a linear combination of the other two sides.) This implies that the velocity vector u (momentum per unit volume normalised by the density) at an arbitrary location within the cell can be described by a linear combination of the normal components u f at the cell faces; see also next section. However, in the case of rectangles, both independent velocity components, u and v , are required to characterize momentum in the two dimensional space due to their mutual orthogonality.

The covolume method
The main merit of the above staggered treatment of the primitive variables is the proven robustness and efficiency as these variables are located in optimal positions on the mesh with two unique consequences. First, it prevents the odd-even decoupling between the pressure and the normal velocity component. Second, it yields excellent conservation properties. For example, both primary (mass and momentum) and secondary (e.g. kinetic energy and enstrophy) quantities are conserved by Cartesian staggered grid methods for the incompressible Navier-Stokes equations (see, e.g. [14,32,33] ). The conservative behaviour of staggered schemes is due to the discrete divergence and gradient operators that mimic the physical properties of their continuous counterparts. As a consequence, the physical accuracy is preserved at the discrete level. In this context, the staggered discretization is classified as the mimetic approach [15] .
The covolume method [12,13] is one of the many mimetic discretizations that have appeared in literature, of which an excellent overview is given by Lipnikov et al. [15] . The mimetic character of discrete covolume operators provides conservation properties in the resulting numerical methods. For instance, the local and global conservation properties of the covolume scheme applied to the incompressible Navier-Stokes equations have been extensively examined by Perot [14] and Zhang et al. [34] .
The covolume method usually requires two orthogonal meshes to compute the change in mass and momentum of shallow water flows. The Delaunay mesh and its dual, the Voronoi tessellation, are the obvious choice. The vertices of the Voronoi mesh are the circumcentres of the Delaunay cells. Like the Cartesian staggered approach, only the normal vector component is employed at the face midpoint of the primal grid and the pressure is consequently defined at the cell circumcentres, i.e. the vertices of the dual grid. In contrast, for instance, the method implemented in EL-CIRC [35] is not based on the covolume approach as it requires the solution of the velocity component tangential to the cell face. Such a semi-staggered scheme can give rise to erroneous pressurevelocity oscillations [36] . In general, preservation of the key physical variables, i.e. the pressure at the cell centre and the normal velocity component u f at the cell face rather than the full momentum vector u , in distinct discrete structures such as primal and dual meshes is crucial to prevent non-physical artefacts. (See Tonti [16] for details on the central role of the relationship between physical variables and mesh topologies, i.e. points, lines, surfaces and volumes, in computational physics.) In addition, vector fields can be reconstructed by means of an interpolation method. Examples related to triangular meshes can be found in [37][38][39][40] .
It is important to note that the use of covolume or C-grid discretizations might be problematic for the simulation of geophysical flows [41] . For instance, spurious geostrophic modes may occur on both structured and unstructured meshes in the presence of the Coriolis force and varying bathymetry. A proper reconstruction of the tangential velocity is the key solution to this problem [39] . Furthermore, three-dimensional hydrostatic models for resolving large-scale baroclinic currents tend to exhibit checkerboard pattern in the divergence of the horizontal velocity field, and hence in the vertical velocity (cf. Eq. (3) ), due to the combined use of triangular geometry, notably equilateral triangles, and variable staggering [41,42] . Some remedies have been put forward in alleviating this problem such as divergence averaging, increased diffusion, filtering and mimetic discretizations (see, e.g. [42][43][44] ), though the unwanted modes are usually not entirely absent. Nevertheless, such checkerboard modes have not been encountered in the numerical wave simulations as presented in this study. Therefore, no further attention will be paid to this issue in this paper.
Although the covolume method is known to be first order accurate for general unstructured meshes and second order for a class of regular grids 2 , it yields physically consistent solution to the governing PDEs. Furthermore, discretizations should reflect accurate changes in the response of the water wave motion to advective acceleration and pressure gradient and, in turn, bathymetric forcing. From this perspective, we postulate that low-order mimetic discretizations are a viable alternative to high-order finite-volume discretization schemes, including colocated and semi-staggered ones employing the full momentum vector as the primary unknown. Despite their higher order accuracy, numerical methods that do not meet certain physical and mathematical properties can produce unacceptably large errors due, principally, to the presence of nonphysical parasitic modes and the nonlinear amplification of such modes and discretization errors [32] . The above reasoning motivates the implementation of the covolume scheme in SWASH which is presented in the next section.

Covolume discretization
We consider a two-dimensional grid consisting of Delaunay triangles. Note that the faces in the dual mesh are orthogonal to the cell faces in the primal mesh. Let indices c and f enumerate triangular cells and faces, respectively. In the discretization approach we denote S c the set of faces forming the boundary of cell c and vector n c,f the outward-pointing normal to face f of cell c . Rather than employing the outward normal, the unit normal to face f coincided with the streamwise direction n f is adopted. Since this normal can either be outward or inward with respect to cell c , we select one of its direction to be fixed at all faces in the grid, so that u f · n f = u f with u f the velocity vector at face f . The mutual orientation of the outward normal n c,f and the unique normal n f at face f of cell c is given by (cf. Fig. 1 ) n c, f = α c, f n f , α c, f = n c, f · n f = ±1 The continuity equation (1) Next, the momentum equation (4) is solved at the faces of each interior cell. The approach to follow is essentially a finite difference method as the spatial discretization will be carried out directly for the equation instead of its integral form expressing the conservation balance of momentum. Fig. 1 shows cell face f and its two adjacent cells c L and c R , with the subscripts indicating the position of the cells with respect to the face.
The distance between the associated circumcentres is given by with s c,f the distance from the circumcentre of cell c to face f . We compute the rate of change of the streamwise momentum per unit face length between two circumcentres adjacent to face f , i.e. s f h f u f , with u f the velocity component normal to face f and h f the water depth at face f . This is determined by, amongst others, the pressure forces acting on the area s f (per unit face length) and the advective flow acceleration altering spatially the momentum in the upstream part of the area, i.e. s c,f h f u f (per unit face length), either with c = c L (if u f > 0) or c = c R (if u f < 0). This upwind bias ensures the stability of the present method since it inherits certain physical conditions, as we will see in the next section.
By virtue of the mesh orthogonality property, we obtain the following finite difference form where a f is the component of momentum advection in the streamwise direction. This will be treated in Section 3.4 . Furthermore, p c , h c and d c denote the pressure, the water depth and the bed level located at the circumcentre of cell c , respectively. Finally, ˜ h f , h f , ˜ p f and ˜ u f are the interpolated water depths, non-hydrostatic pressure and velocity vector at face f , respectively; see below. Note that subscript b has been dropped from the pressure variables. It is worth noting that, owing to the mesh orthogonality, the discrete pressure gradients of Eq. (6) and the discrete divergence of the mass flux, Eq. (5) , are mutually adjoint. This duality contributes to the robustness of the covolume method (see, e.g. [15] for details).
We now provide appropriate interpolations for h f , p f and u f , respectively. Let us consider a flow across an abrupt bed change at cell face f . A condition that must be met is the hydrostatic balance between the (hydrostatic) pressure flux and the topography term where the operator [ · ] f denotes the jump across face f . This condition preserves the quiescent water over varying bathymetry, i.e. u = 0 and ζ = constant. Hence, Eq. (6) reduces to By defining so that which is the hydrostatic balance at the discrete level.
Next, ˜ h f in the time-derivative term of Eq. (6) represents the volume per unit area at face f . Hence, for consistency, we define Likewise, the non-hydrostatic pressure at face f is approximated by taking the volume-weighted average of the pressures of the two cells adjacent to the face Since the covolume scheme employs the face normal velocity components u f , a velocity vector u = (u, v ) needs to be reconstructed from these components. One reconstruction method is due to Perot [14] and computes the cell-based velocity vector u c by means of the divergence theorem. The result is given by which expresses the cell velocity vector in terms of the normal velocity components at the faces of the cell. Perot [14] argued that computing the cell-based velocity vector by means of interpolation (8) implies discrete conservation of momentum (see also Appendix A ). This interpolation is first order accurate and second order if the grid is regular. An alternative is the method proposed by Vidovic [38] using a least squares technique which provides better accuracy, but is rather involved and computationally demanded. We therefore prefer the use of Perot's interpolation scheme (see also [44] ). Finally, the velocity vector at the face can be found as the volume-weighted average of the cell-based vectors of the two cells adjacent to the face, as follows

Discretization of momentum advection
In [45] , the incorporation of the Rankine-Hugoniot jump relations in a numerical scheme is devised to preclude the odd-even decoupling in the velocity field. In addition, fulfillment of these jump conditions is desired to accurately capture flows with discontinuities. In this regard, the advective flow acceleration is projected on the direction n f at face f representing a local shock, a f = a c · n f with a c the upstream cell-based advection vector. The manner in which this nonvolumetric vector is to be determined is closely linked to the fulfillment of the jump conditions at the discontinuity. More specifically, consider steady, frictionless flow across a jump at face f , under the assumption of hydrostatic pressure, then Eq. (6) becomes and a c is constructed such that Eq. (9) reduces to the following Rankine-Hugoniot relation with q the mass flux per unit width. This will improve physical accuracy at the discrete level since the balance between the advective acceleration and the pressure gradient is respected. Additionally, such a scheme remains accurate on coarse meshes, whereas any consistent scheme that does not adhere to jump conditions explicitly converges relatively slower to a weak solution. Examples where such a latter scheme has been employed are [14,22,23] . Thus, fulfilling the Rankine-Hugoniot condition (10) has the potential to further enhance the robustness and efficiency of the numerical method [45] . To achieve the desired jump condition at face f , we use onesided discretizations for momentum advection [45] . First, we consider a cell c upwind of face f whereby velocity u f is pointing out of the cell. Let choose this cell c L , i.e. u f > 0 (cf. Fig. 1 ). Vector a c L is the cell-based advection vector of the upwind cell and is computed using Gauss' theorem The transported flow velocity vector û k at the faces k of cell c L is obtained by means of a one-sided interpolation (see, e.g. [19,20] ) with u c L ,k the cell-based velocity vector of the upwind cell c L of face k (see Fig. 1 ; keep in mind that u f > 0). In turn, this cellbased velocity vector is interpolated from face normal components using Eq. (8) . Finally, we arrive at the following expression for the face normal component of momentum advection Approximation in the case of negative flow, i.e. u f < 0, is done similarly. This completes our discretization of the momentum advection.
We now demonstrate that jump condition (10) is indeed fulfilled. Substituting Eq. (11) in Eq. (9) and using Eq. (7) gives Let k 1 , k 2 and f denote the faces of the upwind cell c L , see Fig. 1 .
the rate of volume flow per unit width of the upwind cell. Hence, we have which is the Rankine-Hugoniot jump condition at the discrete level. Note that this condition holds on face f representing the discontinuity, whereas the first two terms express the jump upstream of that face. For reasons of efficiency, Eq. (6) is not solved as it stands, but the solution of an equation for flow velocity u f is considered instead. For the purpose of deriving this equation, the right hand side of Eq. (6) is set to zero. Furthermore, we consider the case of positive flow at face f , i.e. u f > 0. We obtain as the outflow at face f only alters the water depth upstream . Next, substituting Eq. (5) yields and subsequently, the component of momentum advection a f . We then have The outgoing flux at face f can be omitted without changing the amount of momentum because u c L , f · n f = u f . Hence, including the right hand side of Eq. (6) , the final semi-discretized momentum equation is given by (cf. Fig. 1 ) (12) A similar semi-discrete equation can be derived for the case with negative flow ( u f < 0). Note that the current discretization differs from the one proposed by Kramer and Stelling [19] in that it is able to resolve flow discontinuities more accurately. This is demonstrated in Section 4.2 .
An attractive characteristic of the resulting finite difference approximation, apart from preserving physical accuracy, is that it ensures total momentum is conserved at the discrete level. Such a mathematical property is desired to keep the overall discretization error low by eliminating the conservation errors. A proof is given in Appendix A .

Wetting and drying
The mass flux integrated along a cell face in Eq. (5) is evaluated with where ˆ h f is the water depth approximated at face f . Since the water depth is stored at the circumcentres an upwind-biased interpolation is applied, as follows so that the water depth in the outgoing mass flux of any cell is that of the cell itself. As a consequence, a sufficient condition to guarantee a non-negative water depth in cell c can be derived. Details can be found in Kramer and Stelling [19] . This condition is given by and so it forces at most one cell per time step to be flooded. If a wet cell becomes dry, the water depth in the cell can be arbitrary small, while there is no outgoing mass flux. The flow velocity at face f is then set to zero if ˆ h f < δ with δ the threshold depth ( δ = 10 −10 m in the present study), while momentum equation (12) is not solved. This completes our description of the wetdry algorithm.

Introduction
This section presents results for four test cases with an increasing degree of difficulty in the unstructured mesh configuration. With the exception of the first case, the selection of these cases was motivated by the wish to investigate the main features of the nearshore wave dynamics, such as wave transformation due to shoaling, refraction, diffraction and runup, and nonlinear processes, in particular the resonant three-wave interaction and depth-induced breaking. The objective is to assess the ability of the present unstructured staggered mesh approach to reproduce these wave processes with an emphasis on predictive accuracy. The predictive performance of SWASH for simulating these wave-related cases on rectangular grids is well documented in, e.g. [24,46,47] . The first academic test case is added to this study to verify the numerical accuracy with which the momentum advection is approximated.
Some brief remarks are given here in relation to marching in time and imposing boundary conditions. The semi-discrete Eqs. (5) and (12) are integrated in time using a second order explicit leapfrog scheme in conjunction with the explicit Euler method for the advection term and the implicit Euler method for the non-hydrostatic pressure term [24] . For stability reasons the time step is restricted to comply with the CFL condition that depends on the wave celerity. This condition is given by with t the time step and C f the Courant number evaluated at a face midpoint. A variable time step is employed so that the Courant number is less than a prescribed number ( < 1) in all wet faces. The balance of experience derived from our various waverelated studies suggests that this route has emerged as providing the best compromise in terms of accuracy and stability. For further details on the time integration see Zijlema et al. [24] . Furthermore, boundary conditions considered here pertain mainly to the generation of incident free short waves and bound infragravity waves. They are imposed as weakly reflective at the offshore boundary. Details on their implementation can be found in [24,47] .
As a final note, the staggered mesh method described in this paper will soon be released in a future version of SWASH ( http: //swash.sourceforge.net ), though the source code and a detailed description of the implementation can be obtained from the author upon request.

Dam break over wet bed
The aim is to demonstrate the accuracy with which flow discontinuities are reproduced using the proposed upwind scheme outlined in Section 3.4 and the scheme of Kramer and Stelling [19] for the momentum advection. We consider an unsteady flow with propagation of a dam break wave in a horizontal frictionless channel. The pressure is assumed to be hydrostatic. The length and the width of the basin was 100 m and 10 m, respectively, whereas the dam was initially located in the centre of the basin. The initial upstream water depth was 1 m, and the downstream water depth was set to 0.1 m. The velocity was initially zero. The rather uniform triangular mesh employed has an averaged grid size of 0.4 m, whereas the time step was set to 0.01 s.
At time t = 7 s after dam failure, the model results of the two schemes are depicted in Fig. 2 . Also shown is the analytical solution.
The results obtained with the proposed scheme are in better agreement with the analytical solution than the scheme of Kramer and Stelling [19] . In particular, the simulation shows that a correct speed and height of the bore was virtually rendered by the proposed scheme that satisfies the Rankine-Hugoniot jump conditions.

Wave deformation by a submerged shoal
In this test case, the UnSWASH model is applied to study wave propagation over a submerged shoal on a sloped bottom. The experimental arrangement conducted by Berkhoff et al. [48] has served as a classical test case for this study. The simulation is considered in a basin of 20 m wide and 35 m long with a plane slope of 1/50 on which an elliptic shoal is rested. A monochromatic, unidirectional wave with a period of 1 s and a height of 4.64 cm enters the area and allows to propagate throughout the domain. Detailed information on the bathymetry and the wave conditions can be found elsewhere (e.g. [30,48] ).
The domain was meshed at a resolution in between 0.03 m and 0.07 m, with a total of approximately 370,0 0 0 triangular elements. The mesh is quite uniform as the dominant waves are present everywhere in the domain. Furthermore, it is relatively coarse with 20 to 50 grid cells per incident wave length, while the wave length becomes smaller on the slope and shorter waves are generated over the shoal. The time step was taken initially as 0.005 s. The maximum Courant number was set to 0.8 and the simulation time equaled 30 s.
First, a few words concerning the computational efficiency. The performance of the unstructured-mesh version of SWASH is compared to the structured-grid version. The serial mode simulations were run on an Ubuntu 16.04 desktop with a 64-bit Intel Xeon processor (2.3 GHz). The execution time of UnSWASH needed for the current case appeared to be about 108 CPU minutes, with the highest level of resource required by the solution of the pressure Poisson equation. Furthermore, the total CPU time per grid cell per time step required for UnSWASH was approximately 2.7 μs, whereas for the regular-grid SWASH was about 1.8 μs. This increase in the computation time is probably due to the overhead related to unstructured mesh data structure access, causing low cache efficiency. It is clear that high-resolution UnSWASH must utilize HPC resources in order to be scalable and efficient. A parallel code will be build using message-passing paradigms and will be tested on a commodity computer cluster. This will be reported in a future paper.
Attention is turned next to the comparison between the computed and observed normalised wave heights shown in Fig. 3 .
Wave heights were measured along eight transects at regular intervals (see [48] for details), of which four interesting ones are considered in the present study. The wave focussing in the wake of the shoal (transects 2 and 4), the wave shoaling and refraction on the slope (transects 6 and 7), and the interference pattern caused by diffraction (transect 4) are well captured by the model, despite the use of low-order discretizations. In addition, these results are in good agreement with previously computed results discussed by Stelling and Zijlema [30] .

Breaking waves over a barred topography
The laboratory flume test of Boers [49] is considered with random waves propagating over a barred sandy beach (see Fig. 4 ).
The case 1C is discussed with waves generated at the wavemaker based on a target JONSWAP spectrum with a significant wave height of 0.103 m and a peak period of 3.33 s. This case has a relatively low wave steepness, while waves shoal over the sloping bottom until they break in the surf zone. For detailed analysis on this test case including the infragravity wave dynamics, see Rijnsdorp et al. [47] .
The computational domain is 32 m long and 1 m wide. The non-uniform unstructured grid employed in the present test case consists of approximately 12,0 0 0 triangles with the cell size varying in between 0.2 m offshore and 0.01 m near the coast. This mesh resolution should provide an economical representation of the bathymetric variability in the considered area, while the numerical approach must remain accurate in the presence of relatively large depth gradients. An initial time step of 0.004 s was taken with a maximum Courant number of 0.5. This time step turned out to be unchanged throughout the simulation. Since the present simulation is depth averaged, the so-called hydrostatic   front approximation (HFA) of Smit et al. [46] is employed. Implementation details can be found in Appendix B . Furthermore, a friction coefficient c f = 0 . 001 was adopted. At the offshore boundary, a weakly reflective condition is imposed based on the recorded time series of surface elevation at the most seaward located wave gauge. The moments m n = f n E( f ) df are computed from the energy density spectra E ( f ) (with f frequency) of the surface elevation. Clearly, the trend of the integral wave parameters throughout the domain is well reproduced by the model. Next, we compare the computed wave spectra with the observed ones in Fig. 6 .
The location of the selected wave gauges is indicated in Fig. 4 . Both the generation of higher harmonics and the transfer of energy towards low frequencies are properly predicted. Also, the energy loss at the mid-to high-frequency range due to wave breaking is resolved quite accurately. In general, the model is able to predict the primary characteristics of the spectral evolution, both in the shoaling region and the surf zone.
Though the potential of the present low-order unstructured staggered mesh method is clearly demonstrated, the results can be further improved by using two vertical layers [46] . Extension of the Fig. 6. Observed (thick line) and predicted (thin line) energy density spectra at selected gauges of shoreward propagating waves for case 1C of Boers [49] . unstructured-mesh version of SWASH to include vertical resolution will be done in the future work.

Runup of waves on a conical island
In this last test case the runup of a solitary wave around a conical island is discussed. The experimental configuration of Briggs et al. [50] is selected and is depicted in Fig. 7 . Surface elevations were recorded using wave gauges at eight different locations as in- The grid used to obtain accurate wave runup and inundation on the island contains approximately 171,0 0 0 triangular cells and is locally refined near the island, see Fig. 8 .  Fig. 7 ). The corresponding cases were simulated with UnSWASH using the following settings. The initial time step was 0.001 s with the maximum CFL of 0.8 for the first two cases, whereas for the last severe case they were 0.0 0 05 s and 0.5, respectively. The simulation period was set to 25 s for all cases. Since wave breaking was observed on the lee side of the island for the case H = 0 . 181 d, the HFA is applied accordingly. Fig. 9 provides a comparison between the model results and the measured data at the selected gauges. The observed surface elevations are very well predicted by the UnSWASH model.
In particular, the computed arrival time and height of the incoming wave agrees with the measured data. There is, however, a phase shift observed in the peak at gauge 22, especially for the case H = 0 . 181 d. This known issue can be resolved by including two vertical layers in the model. Finally, it is worth noting that the model results are consistent with previously obtained results for the structured mesh case [24] . The associated mesh is uniform with a grid size of 0.05 m and the number of grid cells is approximately twice as large as that of the unstructured grid.

Conclusions
An extension of the wave-flow model SWASH to unstructured triangular meshes has been discussed. The main motivation for the application of such meshes for the simulation of wave dynamics is the ease of local grid refinement. The covolume method has been adopted for the spatial discretization of the depth-integrated shallow water equations with the primitive variables. The velocity components normal to the cell faces are employed as the primary unknowns in the discretization. To account for wave dispersion the non-hydrostatic pressure is included.
The covolume method has the great benefit of allowing to construct finite difference discretizations that mimic desirable properties of PDEs (e.g. topology, conservation, jump relation) while keeping the computational stencil compact. This significantly improves both the robustness and the efficiency of SWASH. In addition, mimetic discretizations routinely enable physically meaningful results to be obtained on relatively coarse meshes.
Still, the application of the covolume scheme is practically limited to Delaunay-Voronoi meshes, which may impede the user flexibility to generate adequate grids comprising the necessary refinements. However, the associated orthogonality requirement is found not to be a limiting factor in wave-related applications. This is explained by the fact that the required mesh resolution is reasonably non-uniform for the scale of wave dynamics across the entire domain.
A robust and efficient upwind-biased scheme has been proposed. The scheme complies with the Rankine-Hugoniot jump relations and is specifically designed with a view to preserving the local momentum flux. This is crucial to the simulation of breaking waves and unsteady bores. This has been demonstrated by the dam break test case.
The proposed method is generally first order accurate on nonuniform meshes. Despite this fact, the model is capable of reproducing essential wave processes, including shoaling, refraction, diffraction, nonlinear interactions, depth-induced breaking and wave runup, owing to the mimetic nature of applied discretizations and the use of locally fine grids. This has been verified by the three test applications presented in this paper.