SnapWave: fast, implicit wave transformation from offshore to nearshore

This paper presents an efficient, implicit, unstructured-grid wave propagation model, SnapWave (Roelvink, 2025), which provides a simple and fast way to predict nearshore wave conditions at specified locations, for coastline models such as ShorelineS, or wave fields and their forcing of flows, to be used in other models, such as Delft3D-FM, XBeach or SFINCS. We describe the numerical method and verify the correct implementation by comparing against analytical solutions for schematized cases. We then test the model application in four different coastal settings by propagating time series of ERA5 hourly wave conditions to observation points nearshore and through the surf zone. We conclude that the model is robust, easy to set up and fast, and can be applied on open coasts worldwide.

1 Introduction

The simulation of wave propagation and dissipation in coastal areas is important to transform wave fields from offshore areas where wave conditions are available from wave buoys or large-scale wave models to conditions nearshore. The nearshore bathymetry controls the alongshore distribution of wave heights and directions and thereby, to a large extent, the shape and orientation of coastlines. Wave energy dissipation by bottom friction is important when waves pass over large shallow areas, whereas wave breaking dominates the distribution of wave energy through nearshore areas.

One of the earliest grid-based models for nearshore wave propagation and dissipation was HISWA (Holthuijsen et al., 1989), which applied a forward-marching technique on rectangular grids and applied a parameterized frequency spectrum while resolving the directional spectrum. It was fast and reasonably accurate in nearshore areas, but had some disadvantages, such as that the wave grids had to be turned in the wave direction because of the forward-marching technique, and it lacked the full spectral description in frequency and direction. HISWA's successor, SWAN (Booij et al., 1999; Ris et al., 1999), is fully spectral and had third-generation wind growth terms similar to ocean wave models such as WAM (WAMDI-Group, 1988), WAVEWATCH (Tolman, 1991) and WAVEWATCH III (Tolman et al., 2002). The SWAN model runs on curvilinear grids, which made it relatively easy to couple with morphological models such as Delft3D (Lesser et al., 2004). Another much-applied stationary wave model, STWAVE (Smith, 2001), can be applied both in half-plane (forward-marching) or full-plane (all directions) mode; it is fully spectral, with simplified wind-wave growth formulations, and operates on a regular grid.

Although curvilinear meshes allow some flexibility in providing resolution where needed, this is limited because the meshes are structured, 2D in space. Unstructured meshes, whether triangular or a combination of quadrilaterals and triangles, do not have this constraint and can have local refinements that do not radiate out. Several models have been developed using such unstructured meshes: TOMAWAC (Benoit et al., 1997), unSWAN (Zijlema, 2010), WWM (Roland et al., 2012). While these models run much more efficiently thanks to the efficient distribution of grid resolution (see e.g. Alves et al., 2022), the fact that they fully resolve the evolution of the 2D wave spectrum, and apply complex four-wave and triad interactions, still makes them quite computationally expensive when applied over larger areas with high resolution in space and time. To overcome such constraints (O'Reilly et al., 2016) designed a characteristics-based method to link wave conditions outside the surf zone to a network of observation buoys along the California coast. This method resolved wave refraction, shoaling and sheltering of all components of the 2D spectrum, considering the travel time between the offshore buoys and the nearshore, while neglecting wave growth and dissipation, and provided reliable and fairly accurate predictions of nearshore wave conditions in most of the target locations.

For the modelling of nearshore morphological changes and for coastal inundation modelling, the wave model, inside a system or coupled with it, is usually the most time-consuming component. This is mostly due to processes, such as fully spectral modelling of non-linear interactions, that are only of secondary importance in nearshore areas. Even wave growth by wind can often be neglected on open coasts, when the waves only need to travel over tens of kilometres after having been generated over hundreds to thousands of kilometres. Therefore, we aim for a fast, stationary solver capturing only the essential physics of wave refraction and shoaling, dissipation by bottom friction and wave breaking. This solver, called SnapWave, presently serves the following purposes:

• an unstructured solver to resolve wave conditions along a nearshore depth contour, for coastline modelling in ShorelineS (Roelvink et al., 2020)

• an improved stationary wave solver for XBeach (Roelvink et al., 2009), allowing wave propagation in all directions;

• a stationary wave solver for unstructured grids consisting of triangular and quadrangular cells in Delft3D-FM (Reyns et al., 2023);

• a fast nearshore wave solver coupled with SFINCS, to resolve wave setup in inundation modelling (Leijnse et al., 2024).

In this paper, we present the first stage of this model, suitable for most open coasts. The main point of the paper is to demonstrate the use of SnapWave to efficiently transform wave conditions provided by a global model such as ERA5 (Hersbach et al., 2020) to nearshore locations anywhere in the world. After describing the numerical method, we verify the model implementation by comparing it with an analytical solution for refraction and shoaling on a straight coastline, and by analysing the iteration process for the case of a circular island and a circular reef. We then proceed with field cases of varying complexity, followed by a discussion and conclusions.

2 Model description

2.1 Coupled wave action balance and wave energy balance

We solve the wave action balance:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathrm{aa}}{\partial t}}+{\displaystyle \frac{\partial \mathrm{aa}{C}_{\mathrm{g}}\cos \mathit{\vartheta }}{\partial x}}+{\displaystyle \frac{\partial \mathrm{aa}{C}_{\mathrm{g}}\sin \mathit{\vartheta }}{\partial y}}+{\displaystyle \frac{\partial \mathrm{aa}{C}_{\mathit{\vartheta }}}{\partial \mathit{\vartheta }}}={\mathrm{ss}}_A-{\mathrm{dd}}_A\end{array}$$

With aa the frequency-integrated, directionally distributed wave action density:

$$\begin{array}{}\text{(2)}& \mathrm{aa}={\displaystyle \frac{\mathrm{ee}}{\mathit{\sigma }}}\end{array}$$

where ee is the directionally distributed wave energy density and σ is the relative angular frequency. C_g the group velocity, ϑ the wave direction, ss_A the wind source term and dd_A the directionally distributed dissipation. In case the wind-growth source term is included in the balance, the wave period cannot be assumed spatially uniform over the domain and its evolution over the interior of the domain needs to be modelled as well. This is done by simultaneously solving the action balance and the wave energy balance:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial \mathrm{ee}}{\partial t}}+{\displaystyle \frac{\partial \mathrm{ee}{C}_{\mathrm{g}}\cos \mathit{\vartheta }}{\partial x}}+{\displaystyle \frac{\partial \mathrm{ee}{C}_{\mathrm{g}}\sin \mathit{\vartheta }}{\partial y}}+{\displaystyle \frac{\partial \mathrm{ee}{C}_{\mathit{\vartheta }}}{\partial \mathit{\vartheta }}}={\mathrm{ss}}_E-\mathrm{dd}\end{array}$$

where dd is the wave dissipation density and ss_E the wind source term. The representative frequency is then calculated as:

$$\begin{array}{}\text{(4)}& \mathit{\sigma }={\displaystyle \frac{E}{A}}\end{array}$$

where

$$\begin{array}{}\text{(5)}& A=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{aa}\mathrm{d}\mathit{\vartheta },E=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{ee}\mathrm{d}\mathit{\vartheta },\end{array}$$

Similar to HISWA and XBeach, we apply a parameterized frequency spectrum represented by a single frequency close to the peak frequency. The directional spectrum is resolved with a given, constant directional resolution and directional sector.

2.2 Simplification to wave energy balance

The main purpose of this paper is to present the numerical method of SnapWave and its application to propagating wave conditions from deep water to nearshore areas over relatively short distances. In many cases the effect of wave growth by wind and ambient currents is small, and here we focus on such cases. In a forthcoming paper we will detail the method of including wave growth by wind, which is relevant in lakes, estuaries and tidal inlets where the wave climate is dominated by local wind waves.

Under the assumption of no wave growth by wind and neglecting ambient currents, the wave frequency is constant over the domain and the wave action balance (1) reduces to the wave energy balance (3), and the wind input term ss_E goes to zero.

We assume that the directional distribution of the dissipation density dd is the same as the distribution of the wave energy density ee, so:

$$\begin{array}{}\text{(6)}& \mathrm{dd}={\displaystyle \frac{D}{E}}\mathrm{ee}\end{array}$$

Here, the total dissipation D is the sum of the dissipation by wave breaking, D_w and by friction, D_f.

$$\begin{array}{}\text{(7)}& D=D_{\mathrm{w}}+D_{\mathrm{f}}\end{array}$$

The wave breaking dissipation integrated over the directional spectrum, D_w is according to Baldock et al. (1998):

$$\begin{array}{}\text{(8)}& \begin{array}{l}{D}_{\mathrm{w}}=\mathrm{0.25}\mathit{\alpha }\mathit{\rho }g{f}_{\mathrm{p}}\exp \left(-\frac{H_{max}^{\mathrm{2}}}{H_{\mathrm{rms}}^{\mathrm{2}}}\right)\left(H_{max}^{\mathrm{2}}+H_{\mathrm{rms}}^{\mathrm{2}}\right)=\\ =\mathrm{2}\mathit{\alpha }{f}_{\mathrm{p}}\exp \left(-\frac{E_{max}}{E}\right)\left(E_{max}+E\right)\end{array}\end{array}$$

Here, α is a dissipation coefficient of order 1, ρ the density of water, f_p the peak frequency, E the wave energy integrated over the directional spectrum, $E=\mathrm{1}/\mathrm{8}\mathit{\rho }g{H}_{\mathrm{rms}}^{\mathrm{2}}$, H_rms the root-mean-square wave height, $E_{max}=\mathrm{1}/\mathrm{8}\mathit{\rho }g{H}_{max}^{\mathrm{2}}$ and H_max a depth- and frequency-dependent maximum wave height given by:

$$\begin{array}{}\text{(9)}& H_{max}={\displaystyle \frac{\mathrm{0.88}}{k}}\tanh \left({\displaystyle \frac{\mathit{\gamma }kh}{\mathrm{0.88}}}\right)\end{array}$$

The dissipation by bottom friction is given by (Collins, 1972):

$$\begin{array}{}\text{(10)}& D_{\mathrm{w}}=\mathrm{0.28}\mathit{\rho }{f}_{\mathrm{w}}{u}_{\mathrm{rms}}^{\mathrm{3}},u_{\mathrm{rms}}={\displaystyle \frac{\mathit{\omega }{H}_{\mathrm{rms}}}{\mathrm{2}\sinh \left(kh\right)}}\end{array}$$

We can write Eq. (3) in simpler form if we consider s to be the distance along each wave direction:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial \mathrm{ee}}{\partial t}}+{\displaystyle \frac{\partial \mathrm{ee}{C}_{\mathrm{g}}}{\partial s}}+{\displaystyle \frac{\partial \mathrm{ee}{C}_{\mathit{\vartheta }}}{\partial \mathit{\vartheta }}}+\mathrm{dd}=\mathrm{0}\end{array}$$

In the absence of currents, the refraction speed C_ϑ is only governed by the gradients in water depth:

$$\begin{array}{}\text{(12)}& C_{\mathit{\vartheta }}={\displaystyle \frac{\mathit{\sigma }}{\sinh kh}}\left({\displaystyle \frac{\partial h}{\partial x}}\sin \mathit{\vartheta }-{\displaystyle \frac{\partial h}{\partial y}}\cos \mathit{\vartheta }\right)\end{array}$$

2.3 Numerical grid

The unstructured numerical grid may consist of any combination of triangular or quadrangular cells (faces). The grid is defined by a list of grid points, with x, y and depth coordinates, and a list of cells, each with 3 or 4 node numbers. The x and y coordinates can be in a projected coordinate system (in m) or in WGS84 spherical coordinates (in decimal degrees). All values are defined in the nodes of the grid, so no staggering is applied.

This grid definition includes fully triangular meshes, rectangular or curvilinear meshes and stepwise refined rectangular meshes where the transitions to finer resolution are filled by triangles. No orthogonality is assumed. The grid can be specified as an ASCII file with node definitions and cell definitions, or as a NetCDF UGrid file, where only the node coordinates and the face connectivity are used.

2.4 Discretization and solution method

Equation (11) can be discretized as follows:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& {\displaystyle \frac{{\mathrm{ee}}_{k,i\mathit{\vartheta }}^{n+\mathrm{1}}-{\mathrm{ee}}_{ik,i\mathit{\vartheta }}^n}{\mathrm{\Delta }t}}+{\displaystyle \frac{c_{g,k}{\mathrm{ee}}_{k,i\mathit{\vartheta }}^{n+\mathrm{1}}-c_{g\mathrm{u},i\mathit{\vartheta }}{\mathrm{ee}}_{\mathrm{u},i\mathit{\vartheta }}^{n+\mathrm{1}}}{\mathrm{\Delta }{s}_{k,i\mathit{\vartheta }}}}\\ & +{\displaystyle \frac{c_{\mathit{\vartheta },k,i\mathit{\vartheta }}{\mathrm{ee}}_{k,i\mathit{\vartheta }+\mathrm{1}}^{n+\mathrm{1}}-c_{\mathit{\vartheta },k,i\mathit{\vartheta }-\mathrm{1}}{\mathrm{ee}}_{k,i\mathit{\vartheta }-\mathrm{1}}^{n+\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }\mathit{\vartheta }}}\\ & +{\displaystyle \frac{D_k}{E_k}}{\mathrm{ee}}_{k,i\mathit{\vartheta }}^{n+\mathrm{1}}=\mathrm{0}\end{array}\end{array}$$

where k is the grid node number, iϑ the direction bin number and n the timestep/iteration number. The subscript u refers to the point upwind of grid point k, as illustrated in Fig. 1; values for cg and ee in this point are obtained from the two points p, which are upwind from point k for directional bin iϑ, and with weights w that depend linearly on the distance between the upwind point and the two adjacent points p.

Figure 1 Schematic showing the relation between point k and its upwind points p.

Wave energy is conserved since we apply the conservative form of the wave energy balance. We can write the system of equations per grid point as:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& A{\mathrm{ee}}_{k,i\mathit{\vartheta }-\mathrm{1}}^{n+\mathrm{1}}+B{\mathrm{ee}}_{k,i\mathit{\vartheta }}^{n+\mathrm{1}}+C{\mathrm{ee}}_{k,i\mathit{\vartheta }+\mathrm{1}}^{n+\mathrm{1}}\\ & =R({\mathrm{ee}}_{k,i\mathit{\vartheta }}^n,{\mathrm{ee}}_{\mathrm{prev},i\mathit{\vartheta }}^{n+\mathrm{1}})\end{array}\end{array}$$

Here, the coefficients, depending on the choice for an upwind or central scheme, are given by:

$$\begin{array}{}\text{(15)}& \begin{array}{ll}\text{Upwind},c_{\mathit{\vartheta },k}>\mathrm{0}& \text{Upwind},c_{\mathit{\vartheta },k}<\mathrm{0}\\ A=\frac{-c_{\mathit{\vartheta }k,i\mathit{\vartheta }-\mathrm{1}}}{\mathrm{\Delta }\mathit{\vartheta }}& A=\mathrm{0}\\ B=\frac{\mathrm{1}}{\mathrm{\Delta }t}+\frac{c_{g,k}}{\mathrm{\Delta }{s}_{k,i\mathit{\vartheta }}}+\frac{c_{\mathit{\vartheta },k,i\mathit{\vartheta }}}{\mathrm{\Delta }\mathit{\vartheta }}+\frac{D_k}{E_k}& B=\frac{\mathrm{1}}{\mathrm{\Delta }t}+\frac{c_{g,x}}{\mathrm{\Delta }{s}_{k,i\mathit{\vartheta }}}-\frac{c_{\mathit{\vartheta },k,i\mathit{\vartheta }}}{\mathrm{\Delta }\mathit{\vartheta }}+\frac{D_k}{E_k}\\ C=\mathrm{0}& C=\frac{c_{\mathit{\vartheta },k,i\mathit{\vartheta }+\mathrm{1}}}{\mathrm{\Delta }\mathit{\vartheta }}\\ \begin{array}{l}\text{Central scheme}\\ A=\frac{-c_{\mathit{\vartheta }k,i\mathit{\vartheta }-\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }\mathit{\vartheta }}\\ B=\frac{\mathrm{1}}{\mathrm{\Delta }t}+\frac{c_{g,k}}{\mathrm{\Delta }{s}_{k,i\mathit{\vartheta }}}+\frac{D_k}{E_k}\\ C=\frac{c_{\mathit{\vartheta },k,i\mathit{\vartheta }+\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }\mathit{\vartheta }}\end{array}& \end{array}\end{array}$$

In all cases the right-hand side is given by:

$$\begin{array}{}\text{(16)}& R={\displaystyle \frac{{\mathrm{ee}}_{k,i\mathit{\vartheta }}^n}{\mathrm{\Delta }t}}+{\displaystyle \frac{c_{g,\mathrm{prev}}{\mathrm{ee}}_{\mathrm{prev}}^{n+\mathrm{1}}}{\mathrm{\Delta }{s}_{k,i\mathit{\vartheta }}}}\end{array}$$

This is a tridiagonal system with the dimension of ntheta that can be efficiently solved for each point using a standard Thomas algorithm. The solution for each point relies on having (ideally converged) estimates of the wave energy density ee in the upwind points for each wave direction. For each point in the unstructured mesh, the spatial propagation is solved by backtracing, for each direction, to the line connecting two upwind points, in a manner similar to STWAVE and unSWAN. The combined propagation, refraction and dissipation are solved implicitly for each point. Wetting and drying is handled simply by making points inactive that have $\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\mathrm{aa}\mathrm{d}\mathit{\vartheta }$, depth less than 1.1 times hmin, set to 0.1 m by default.

To arrive at a stationary solution, we set the time step to a very large number and apply a “sweeping” technique as follows. We determine for each sweep s and each point k the distance r_k,s along the mean wave direction ϑ_m, the two orthogonal directions and the opposite direction as follows:

$$\begin{array}{}\text{(17)}& \begin{array}{l}{r}_{k,s}=x_k\cos \left({\mathit{\vartheta }}_{\mathrm{m}}+\frac{\mathit{\pi }}{\mathrm{2}}{S}_s\right)+y_k\sin \left({\mathit{\vartheta }}_{\mathrm{m}}+\frac{\mathit{\pi }}{\mathrm{2}}{S}_s\right),\\ \text{all }k,s=\left[\mathrm{1}:\mathrm{4}\right]\\ {\mathrm{S}}_s=[\mathrm{0},\mathrm{1},-\mathrm{1},\mathrm{2}]\end{array}\end{array}$$

Next, we sort the points in all four directions and store the index for all points and each sweep direction. For each sweep, this index determines the order in which we solve Eq. (14).

Generally, the first sweep already solves a major part of the wave propagation, as forward-marching techniques would do.

Secondary effects of refraction are covered by “sweeping” in all 4 directions. Since the wave dissipation is a very nonlinear function of the wave height and water depth, the whole system needs to iteratively come to a converged solution. Convergence is checked after all four sweeps in an iteration; points, where the maximum difference in energy density divided by the maximum energy density for that point is less than a user defined threshold crit, are fixed and taken out of the loop. The iteration is converged when the maximum difference in energy density, normalized by the maximum energy density, is below the same crit. As the method converges rapidly and we take out the points already converged, we can set the default crit at a comfortably low value of 10⁻⁵. The process generally converges within 4–6 iterations.

3 Verification

In comparison of model to theory or data, the error metrics Pearson's correlation coefficient (rho), scatter index (sci), relative bias (relbias), and Brier skill score (skill) are computed as shown in Table A1.

3.1 Linear shoaling and refraction

A first test of the correct implementation of refraction and shoaling is to compare the wave height and mean wave direction over a longshore uniform, double barred profile. We use an analytical representation of typical Dutch barred profiles (Bakker, 2013):

$$\begin{array}{}\text{(18)}& z_{\mathrm{b}}=z_{\mathrm{r}}-A{x}^{\mathrm{b}}-A_{\mathrm{b}}{e}^{-{\left(\frac{x-x_{\mathrm{b}}}{R_{\mathrm{b}}}\right)}^{\mathrm{2}}}\cos \left(\mathrm{2}\mathit{\pi }\left({\displaystyle \frac{x}{L_{\mathrm{b}}}}-{\displaystyle \frac{t}{T_{\mathrm{b}}}}\right)\right)\end{array}$$

With z_b the bed level, z_r a reference level of 6 m, A_b the bar amplitude of 1 m, R_b the bar scale of 200 m, x_b the location of maximum bar amplitude (300 m), L_b the bar wavelength (200 m) and T_b the bar migration period (10 years). The expression describes a bar system that grows, migrates seaward and damps in a periodic fashion; the time t was taken arbitrarily as 0 years, which means that the bar crest is at the location of maximum amplitude. The water level was set at 0 m, and as purely refraction and shoaling, no breaking, were tested, the depth was cut off at 1 m and the breaker parameter gamma set to a high value.

Three grid configurations were applied: one with a uniform grid size of 20 m by 20 m (denoted “uniform_20”), one with uniform resolution of 10 m by 10 m (denoted “uniform_10”) and one where the resolution varied from 40 to 10 m through two uniform refinements (“variable_40_10”). The domain was 2000 m cross-shore by 10 000 m cross-shore; the coarse uniform grid had 50 000 nodes, the fine uniform grid had 200 000 nodes and the non-uniform grid approximately 31 000 nodes.

Uniform boundary conditions were specified on the offshore boundary and Neumann boundary conditions (no longshore gradient) at the lateral boundaries. The boundary conditions and model settings are specified in Table 1.

Table 1 Parameters shoaling and refraction test.

Figure 2 Refraction and shoaling test; comparison between analytical solution for 0, 45 and 60° angle of incidence (drawn black lines) and SnapWave results for uniform grid (blue dotted line), fine uniform grid( red dotted line) and unstructured grid (green dotted line).

The results show a good agreement between the analytical and numerical wave direction and wave height, as is shown in Fig. 2 and in the statistics Tables A2 and A3, with a scatter index of less than 1 % for the wave height and wave direction, and a bias in wave height of less than 1 % and less than 1.2 % in wave direction. As expected, the model that is refined in the barred area has slightly better skill than the coarse uniform grid, comparable to the fine uniform grid, at less than a sixth of the number of nodes; no deviations are found at the transitions in grid resolution.

3.2 Circular island

The circular island testcase is included to illustrate the capability of SnapWave to compute the wave refraction and shoaling all around an island, and to show how the solution scheme progresses. The conditions are taken from the case of a sandy circular island (Kamphuis and Nairn, 1985), with a radius of 350 m, and a 1:12 slope until a depth of 20 m. A circular curvilinear grid was applied with uniform cross-shore resolution of 5 m; the directional resolution was 5° and the sector was 360°. Wave conditions were imposed using an Hm0 wave height of 2 m, a peak period of 15 s and directional spreading of 20°. Various angles of incidence were tried, all uniformly applied on the outer boundary, resulting in symmetrical patterns. In Fig. 3 the computational grid and bathymetry are shown, as well as the wave height distribution for incident wave angle of 270° N. The resulting focusing of the waves in front of the island, and the reduced wave height on the sides and in the back of the island agree well with earlier results shown by (Roelvink et al., 2013) for XBeach, both in stationary and nonhydrostatic mode.

Figure 3 Circular island test. Left panel: bathymetry and computational grid. Right panel: wave height Hm0 distribution. Points 1–4 correspond to locations of directional distributions.

The sweeping process converges rapidly and resolves the wave pattern all around the island. In Fig. 4 the directional distributions of the directional energy density ee are shown for the first 5 sweeps, at four points surrounding the island. As the first sweep is plotted last, a purely green line indicates that all subsequent sweeps are hidden behind it and therefor have not changed much. This is clearly the case for point 1 on the windward side, where the first sweep going from East to West is almost fully converged. In point 4, sweep 2 going from South to North almost fully resolves the distribution. It modifies the peak in point 2 but not completely, and it adds the purple peak in point 3 at the leeward side. Sweep 3 proceeding from North to South produces the second peak at point 3, and brings the peak in point 2 to the same level as in point 4. In point 3 at the lee of the island, sweep 4 brings the peak at around 40° at the final level. Subsequent sweeps and iterations have very little impact and quickly converge to high accuracy, as indicated in Table 2.

Table 2 Convergence characteristics of circular island test.

An interesting aspect of the solution is that at the leeward side we have waves from almost opposing directions. In the nonhydrostatic solution in Roelvink et al. (2013) this was also observed.

Figure 4 Circular island test. Change in directional distribution after first 5 sweeps, first one plotted last, for the four indicated points.

We can conclude that the wave patterns are realistic and that the method quickly converges for waves incident from any direction.

3.3 Circular reef

The circular reef case was inspired by the work of Mandlier and Kench (2012) who considered analytical solutions to the refraction problem using ray tracing. The case we present has a flat circular reef with a radius of 350 m, a depth of 1.5 m on the reef and deep water (taken as 100 m) all around it.

To be able to compare our model with the analytical solution in terms of wave height distribution, we reproduced the wave ray refraction pattern as described in Mandlier and Kench (2012) and added the computation of wave heights, by counting the number of wave rays passing within a certain distance, taken as 4.5 times the initial ray distance, from each grid point in a regular 5 m by 5 m grid, and computing the wave height as the rms value of the wave heights associated with each refracted wave ray within this distance. The incident wave height was 0.1 m, ensuring wave breaking did not play a role. The resulting refraction pattern and wave height distribution are shown in Fig. 5, showing a highly concentrated wave height region around 90 m East of the centre of the reef, and two areas of very low or undetermined wave heights where the wave rays cannot reach.

Figure 5 Reproduction of Mandlier and Kench analytical solution for a flat circular reef; left panel: refraction pattern of wave rays; right panel: derived wave height.

For SnapWave, we constructed a circular grid with 5 m radial resolution and 1° angular resolution, except for the part within 50 m of the centre, which was filled in with triangles with sides around 2.5 m, resulting in the grid shown in the left panel of Fig. 6. The grid was rotated to check whether the implementation was sensitive to the grid orientation, which it was not. The wave conditions were a mean direction of 270° N, a peak period of 12 s and a small directional spreading of 5°. The wave angle resolution was 1°. A small incident wave height of 0.1 m was applied to enable a comparison with the analytical solution of Mandlier and Kench (2012). In the right-hand panel the wave height distribution is shown, where we see a narrow area of concentrated wave height at around 60–150 m from the centre of the reef. This corresponds reasonably well with the area of concentrated wave energy in the analytical solution, which centred around 90 m from the centre. It must be noted that our model provides seemingly reasonable results on the leeward side of this caustic and does not blow up; for higher incident waves the wave breaking mechanism kicks in and limits the growth of the wave height near the caustic.

Figure 6 Reef refraction test. Left panel: bathymetry and computational grid. Right panel: wave height Hm0 distribution. Black circle indicates maximum wave height; points 1–4 correspond to locations of directional distributions.

The building up of the solution is almost complete in the first sweep, for points 1, 2 and 3. For points 2 and 4 the fourth, East to West sweep brings some additional energy peak from almost easterly direction, due to refractive trapping along the edge of the reef. In any case it is symmetrical and relatively small.

We may conclude that, although there is not a perfect match, the SnapWave model produces a very similar wave height pattern at the windward side and an area of highly focused wave height over an area similar to the analytical solution. Interestingly enough, the SnapWave method is considerably faster than the analytical solution.

Figure 7 Change in directional distribution after first 5 sweeps, first one plotted last, for the four indicated points.

4 Field validation

The main objective of the field validation cases is to demonstrate a methodology to hindcast or predict nearshore wave conditions based on ERA5 data at locations ∼ 0.5° offshore, global or local bathymetry and SnapWave to transform the wave conditions to specified nearshore points. We consider four testcases, spanning a range of conditions and geographical locations.

• Coast3D campaign at Egmond, the Netherlands, situated on an open, barred coast

• Ameland inlet, the Netherlands, under the influence of a large ebb tidal delta

• St Croix, US Virgin Islands, with operational buoys on either side of the island

• Ningaloo Reef, Australia, with an array of pressure sensors across a wide, shallow reef.

For all cases we use a similar setup starting from ERA5 model output points at 0.5° resolution. For two of the cases, Coast3D and Ningaloo, we compare these results with those of a local model driven by locally measured wave conditions, in order to distinguish between errors in SnapWave and those inherent in the ERA5 model.

4.1 Coast3D

4.1.1 Local model vs dcsm-fine

This testcase concerns the hindcasting of wave conditions at the 15 m depth contour and the subsequent propagation and dissipation of the waves throughout the surf zone at Egmond, the Netherlands. These wave measurements were part of a large EU project COAST3D (Soulsby, 2001); Egmond was one of the test sites and the main campaign at this location, in November 1998, is described extensively in Ruessink (1999); the wave measurements are also detailed in Reyns et al. (2023).

In Fig. 8 the extent of the large-scale model is shown along with the ERA5 output points used as boundary locations. In Fig. 9 the details of both the large-scale model and the local model are shown in the area of the field campaign. The large-scale model has a resolution ranging from approximately 800 to 100 m near the entire coast, with three subsequent local refinements to approximately 14 m in the measurement area. It must be noted that for providing boundary conditions to coastline models typically a grid size of 100 m would be sufficient, but the finer resolution is needed to resolve the breaker bars in the surf zone. The local model has a curvilinear setup with cross-shore resolution from 70 to 13 m and a longshore resolution from 125 to 25 m; in other words, the resolution in the measurement area is similar.

Figure 8 Overview of Holland coast, with bathymetry of large-scale model domain; black rectangle in North Holland: location of Egmond field campaign. Green dots: locations of ERA5 boundary points. Map data: © OpenStreetMap contributors 2025. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

Figure 9 Details of the computational grids for the large-scale model (left panel) and the local model (middle panel) and measurement locations (right panel). Map data: © OpenStreetMap contributors 2025. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

The measurement period considered here covered the period of 1 November, 00:00 UTC until 12 November 12:00, 1998. ERA5 data were downloaded and time series were extracted for the indicated boundary points (Fig. 8) at hourly intervals. Data for point 8, at 15 m depth, were used as validation data for the large-scale model, and as boundary conditions for the local model. Points 1a–1d covered the transformation over the barred profile, as indicated in Ruessink (1999). For the water levels, the measured timeseries, interpolated from 2 tidal stations as in Ruessink et al. (2001) was applied uniformly over both models.

Sensitivity tests indicated that the results were little sensitive to the directional resolution, so a directional step of 10° was applied. Breaker parameter gamma values of 0.75 (default) and 0.70 were applied.

4.1.2 Results local model

First, we discuss the results for the local model, see Fig. 10. At the outermost station 2 we see modest events on 2, 3 and 11 November, and a major event on 6 November. At this location, the wave height variation is mostly due to the variation in offshore wave conditions. Unfortunately, only a short period is available in the observations. As we move through the surf zone in points 1a thru 1d the effect of depth-induced breaking becomes more obvious, leading to a strong tidal modulation in the wave height time series. These results are very similar to those of Ruessink et al. (2001) applying a profile model.

Figure 10 Time series of Hm0 wave height at 30 min intervals, measured (red dots) vs computed (drawn blue line). Local model driven by observed wave conditions at 15 m depth, γ = 0.7.

In Fig. 11 we show scatterplots of the computed vs. observed wave heights, for a gamma value of 0.7, which showed slightly less bias and higher skill than the default value of 0.75. Error metrics for this test are given in Table A4. The results for station 2 show the highest bias and scatter, but it must be noted that these points only cover a short period. The surf zone points 1a thru 1d show a modest scatter index in the order of 15 % and a bias of around 10 %. Reducing gamma further reduces the bias but results in poorer performance for the higher wave conditions. Note that we did not take the substantial bed level changes in the inner surf zone over the course of the measurement campaign into account. Propagating the wave heights through the surf zone is performed with a skill of over 96 %.

Figure 11 Scatterplots (heat maps) of computed vs observed Hm0. Local model driven by observed wave conditions at 15 m depth, gamma = 0.7.

4.1.3 Results large-scale model

The time series of the large-scale model simulation are shown in Fig. 12. First, as an indication of the quality of the ERA5 hindcast, the observations at the point 8 at 15 m water depth are generally reproduced quite well, except for a small event at 3 November, which is completely missed by ERA5; during that period, both wave directions and wind directions are offshore in ERA5 so there is no possibility of getting such nearshore wave heights in the order of 2 m. The other peaks are generally predicted well, sometimes with a phase shift in the order of a few hours.

The results through the surf zone, though less accurate than for the local model, generally reproduce the observed time series quite well, particularly around the main event between 5 and 7 November.

Figure 12 Time series of Hm0 wave height at 30 min intervals, measured (red dots) vs computed (drawn blue line). Large-scale model driven by ERA5 boundary conditions, gamma = 0.7.

The scatterplots shown in Fig. 13 confirm this narrative, as do the error metrics in Tables A5 and A6. Point 8 is mostly indicative of the skill of the ERA5 hindcast and has a low bias of 7 % and a scatter index of 22 %. In point 2 the scatter index is quite high since the short time series includes the event that was missed by ERA5. The points through the surf zone have a bias of around 10 % and higher scatter indices than the local model, mostly for missing the event on 3 November.

Figure 13 Scatterplots (heat maps) of computed vs observed Hm0. Large-scale model driven by ERA5 boundary conditions, gamma = 0.7. For color legend refer to Fig. 11.

4.1.4 Conclusions

From this case study we may conclude that SnapWave has an adequate skill of more than 95 % in propagating waves through the surf zone; for given wave conditions at 15 m depth, the bias is in the order of 10 % and the scatter index in the order of 15 %. When using it to propagate waves from the nearest ERA5 output points, the bias at the 15 m depth point was a low 4 % and the scatter index 21 %. For the subsequent propagation through the surf zone the bias remains low at around 10 % and the scatter index is less than 30 %; the higher values for the scatter index compared to the locally driven model are mostly due to the phase shifts between observations and model at the 15 m depth point, and due to one event being missed by ERA5. According to the used definition, the model skill is around 0.9 and higher.

4.2 Ameland

Ameland Inlet is typical for tidal basins with a large tidal inlet, and the wave refraction, shoaling and dissipation over such ebb delta is important for the coastline evolution along the seaward-facing coast. In the framework of the project SBW (Strength and Loads on Coastal Defences) a number of directional wave riders was installed from 2007 onwards; for an overview see (Elias, 2017). Several studies have focused on the wave penetration into the Wadden Sea and on the effect of wave growth by wind and current refraction, but for the purpose of this study we focus on the wave distribution around the ebb delta. We extracted two-month time series from the MATROOS system used by the Dutch government and knowledge institutes, for the period of 1 November thru 31 December 2008. As there was uncertainty over the location of one of the buoys in early November we used the period of 5 November until 31 December 2008.

We created an unstructured grid covering all the Dutch Wadden islands and extending to the nearest reliable (i.e. not affected by land) ERA5 points, as indicated in Fig. 14. The resolution ranged from 800 m offshore to approximately 100 m in the nearshore. The bathymetry in the area of interest was updated with area soundings (“Vaklodingen”) from 2008. Six observation points were selected, as shown in Fig. 15. The two points AZB11 and AZB12 are outside the ebb delta and are indicators of the quality of the ERA5 hindcast. The other four are spread out over the ebb delta and should give an impression of the quality of the wave propagation model in a complex area with shoaling, refraction and wave breaking.

The tide level was imposed uniformly based on a nearby output location from the Global Tide and Surge Model (GTSM, Muis et al., 2016).

From ERA5 the data for Hm0 wave height, peak period, mean wave direction and directional spreading were extracted. The data for the observation points contained Hm0 wave height and Tm10 wave period, which based on our experience we converted to peak period by multiplying by 1.1.

Default parameter settings were chosen with gamma of 0.75, a directional resolution of 10° and a directional sector of 360°.

Figure 14 Overview of model domain and bathymetry for Ameland hindcast. Green dots indicate ERA5 boundary points; red dots observation points. Map data: © OpenStreetMap contributors 2025. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

Figure 15 Detail of bathymetry and observation points Ameland inlet. Map data: © OpenStreetMap contributors 2025. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

Time series for Hm0 wave height are shown in Fig. 16. Clearly, the points AZB11 and AZB12 are in depths where the tidal modulation does not play a role yet, and generally the computed wave heights follow the measurements closely, except for the event around 22/11 where the wave height is clearly underestimated. This is likely due to an underestimation of wave heights by ERA5 for this event, though some additional wave growth due to wind (not included in this SnapWave model) could play a role as well.

The results for the four other points are clearly modulated by the tide, as wave breaking plays an important role. A change of gamma value to 0.8 improved the error statistics somewhat. Relative bias and scatter index are in the same order of magnitude as for the Coast3D case, around −10 % and 25 %–30 % respectively.

Figure 16 Time series of Hm0 wave height for 6 observation points, Ameland Inlet.

Though SnapWave without wind growth terms assumes a uniform distribution of the peak period, it is useful to test this assumption against the wave data. Figure 17 shows that this is not a very bad assumption; from the error statistics in Tables A7 and A8 we see that the relative bias is around zero and the scatter index in the order of 25 %.

Figure 17 Uniform Tp vs observed 1.2× Tm01, Ameland Inlet.

The scatterplots in Fig. 18 confirm that the systematic underestimation of the higher wave heights originates with the ERA5 data and propagates through the nearshore area. The scatter in the nearshore points is rather consistent at around 30 %. According the used definition, the model skill is consistently over 0.9, indicating an adequate performance for such cases.

Figure 18 Scatterplots SnapWave Hm0 vs observations, coloured by wave direction.

The performance for this model is similar to that of the Coast3d large-scale model; the number of nodes is around 250 000 and the run time per wave condition is around 2.5 s, or approximately 10 µs per node and wave condition.

4.3 St Croix

The island St Croix (US Virgin Islands) is used as a case study where open boundaries are applied at all sides, and ERA5 data are specified all along these boundaries. There are two operational CDIP buoys (https://cdip.ucsd.edu/m/about/, last access: 13 September 2017) at the edge of the shelf, called “Fareham” on the southern end and “Christiansted” on the northern side. In terms of processes needed, the case is not too challenging. We test mainly if the ERA5 hindcast is accurate and if the shielding and for some wave directions the refraction on the shallow reef areas are properly accounted for. The model setup is shown in Fig. 19; the cell sizes range from 800 m offshore to 200 m near the coast; higher resolution was not needed here as the observation points were on the edge of the shelf, still in relatively deep water. We obtained wave height and period records from the CDIP buoys for the period of 1 June until 1 November 2010, and downloaded ERA5 wave data for Hm0 wave height, peak period, mean wave direction and directional spreading for the same period, for the locations indicated by the green dots.

Figure 19 Grid layout and bathymetry, St Croix case. Green dots indicate ERA5 boundary locations; red dots indicate the observation points at the site of the CDIP buoys. Map data: © OpenStreetMap contributors 2025. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

In Fig. 20 the time series comparison is shown for the two observation points. In general, the model follows the observations closely, except for the event on 18 September which is severely underestimated at the Christiansted buoy, and an event on 6 October, which the model underestimates at the Fareham buoy. Such behaviour where ERA5 misses some of the extreme peaks due to a lack of resolution is well documented (e.g. Fanti et al., 2023).

Figure 20 Time series of Hm0 wave height for stations Christiansted and Fareham near the island of StCroix, computed (drawn black line) vs. observed (blue dots).

The scatterplots (heat maps) confirm the fact that for most conditions the agreement is quite good, and only for some individual events the ERA5 model misses the peaks. Overall, as is also apparent from the error statistics, the bias is less than 10 % and the scatter index in the order of 20 %, and skill over 95 %.

Figure 21 Scatter plots (heat maps) of computed vs. observed Hm0 wave heights (left panels) and Tp wave period (right panels), for stations Christiansted and Fareham, St Croix. For color legend refer to Fig. 11.

The computation for the 5 months took 51 min, on average 0.8 s per wave condition or 20 µs per node per wave condition (TBD: clean performance check).

4.4 Ningaloo Reef

Ningaloo Reef is a wide and extensive, pristine coral reef in NW Australia. The reef has been the subject of a number of studies on hydrodynamics and sediment transport, and data collected there (Pomeroy et al., 2012) has been used to validate other wave models such as XBeach in Van Dongeren et al. (2013). That study focused on the generation of infragravity waves but also considered the propagation of the swell waves as we do here. One important finding in these studies was that the roughness of the reef was very high, and could be mimicked by using a high friction factor fw of 0.6 on the reef. Here we used this value, making use of the option to impose space-varying roughness fields as random samples. We focus in a cross-shore transect with pressure sensors C1 on the forereef and C3 through C6 on the reef flat.

We used two model setups: one local model (square cells, resolution 16 m by 16 m) driven entirely by locally measured wave conditions, and one unstructured grid with square cells, refined 5 times, with resolution from 500 to 16 m. The overall model grid is shown in Fig. 22 and details at the measurement site are shown in Fig. 23. For the water level in the large-scale model we extracted time series for a nearby location from, the GTSM (Muis et al., 2016) for the month of June 2009 and imposed this uniformly.

Figure 22 Overview of large-scale Ningaloo reef model with bathymetry, boundary points (in green dots) and the location of the local model. Map data: © OpenStreetMap contributors 2025. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

Figure 23 Detail of large-scale model at measurement site (left panel) and local model (right panel). Observation points in rod dots. Map data: © OpenStreetMap contributors 2025. Distributed under the Open Data Commons Open Database License (ODbL) v1.0.

First, we compare time series of the Hm0 wave height for the local model. As noted in the literature, the swell heights rapidly decay behind the reef edge, a process that is dominated by the bed friction. Also, as in the Coast3D and Ameland cases, the wave heights over the reef flat are strongly modulated by the tidal water level elevation. The model results follow the observations reasonably well, given that the wave heights decrease by an order of magnitude. As SnapWave by itself does not consider the wave setup, water depths on the reef flat are underestimated, which is apparent particularly in the most shoreward points. As shown in the statistics, the bias is in the order of centimetres; the relative bias is in the range of 0 %–25 % as the mean Hm0 on the reef flat is very low. The same holds for the rms error, which is a few centimetres, whereas the scatter index is in the order of 30 %.

Figure 24 Time series of HM0 wave height across the reef at Ningaloo; observations (black dots) against model simulation (blue drawn lines). Local model.

Figure 25 Scatterplots (heat maps) of computed vs. observed Hm0 wave heights. Local model. For color legend refer to Fig. 11.

For the large-scale model the results are shown in Figs. 26 and 27. First, we see that the ERA5 model predicts the general trend in the wave height time series at the outer reef location, but underestimates the Hm0 around yearday 164, and overestimates for much of the remainder of the period, particularly around yearday 174. Still, the relative bias of 10 % at this location and the scatter index of 29 % are in line with the other case studies.

For the reef locations the relative bias is less than 10 % for most locations except C4, but the scatter index is rather high, at 40 %–65 %. This is mostly due to a small phase shift between the GTSM hindcast water level and the observed water level in situ, as can be seen in Fig. 28. When we apply this shift to the simulated model results and compare them with the observations, as shown in Figs. 29 and 30, the skill over the reef flat improves considerably, from 0.72 to 0.84.

Figure 26 Time series of HM0 wave height across the reef at Ningaloo; observations (black dots) against model simulation (blue drawn lines). Large-scale model.

Figure 27 Scatterplots (heat maps) of computed vs. observed Hm0 wave heights. Large-scale model. For color legend refer to Fig. 11.

Figure 28 GTSM hindcast of water level (blue) vs. observed (red). Lower panel: uncorrected except for 8 h shift from GMT to local time; top panel: GTSM model shifted by one hour to GMT + 7 h.

Figure 29 Time series of Hm0 wave height across the reef at Ningaloo; observations (black dots) against model simulation (blue drawn lines). Large-scale model, simulation results shifted by one hour.

Figure 30 Scatterplots (heat maps) of computed vs. observed Hm0 wave heights. Large-scale model. Simulation results shifted by one hour. For color legend refer to Fig. 11.

4.4.1 Effect of wave setup on wave propagation and decay.

For a given water level, the local wave setup can have an important influence on the wave decay. In such cases, a coupled hydrodynamic model is needed to provide this non-uniform water level. We tested the effect of wave setup by using Delft3D-FM with the in-built SnapWave solver to check how important this effect is. The model was set up for the local grid and fed with constant water level boundary conditions, over the same period as for the other simulations. In Fig. 31 the wave heights for the simulation without wave setup are compared with those including the effect of wave setup. The effect is significant, in the order of 2–3 cm or about 10 %–20 % of the local wave height. Still, in light of the uncertainty in the local water level, this is still a relatively minor effect.

Figure 31 Comparison of coupled Delft3D-FM – SnapWave model with and without setup, for observation points 1, 3 and 4.

4.5 Haringvliet intercomparison with unSWAN

The Haringvliet mouth is a shallow area seaward of a closed-off estuary, with sand banks, shoals and channels and has long been used as a test case for the models HISWA, SWAN and recently unSWAN. The model can be readily downloaded from the SWAN website https://swanmodel.sourceforge.io/download/download.htm (last access: 21 October 2025) and was used here to check the difference in results and performance between SnapWave and unSWAN, on the exact same triangular grid. Two modifications were made to the unSWAN input:

• instead of a 1D spectrum, a similar parametric JONSWAP spectrum was imposed, with peak period Tp = 8 s, wave height Hm0 = 3.2 m and mean direction of 270° N; peak enhancement factor was 3.3.

• the bed friction was set to “COLLins” with value of 0.02, to facilitate comparison with SnapWave, fw = 0.02.

• wind was turned off in both cases.

For SnapWave, the unSWAN grid files were converted to a NetCDF UGrid file. In both cases, the same convergence criterion of 0.02 was used. SnapWave converged for 100 % of grid nodes in 3 iterations, which took 0.03 s, whereas unSWAN took 5 iterations to converge for 99.7 % of the grid nodes, in approximately 6 s, a factor of 200 slower.

The results are shown in Fig. 32; in most of the area, the pattern of Hm0 wave height is very similar. We present the point-by-point comparison in Fig. 33, with a low relative bias of 0 %, a scatter index of 9 % and a skill of 0.98.

We may conclude that in this kind of nearshore application dominated by refraction, shoaling, bed friction and wave breaking dissipation, differences between SnapWave and unSWAN are small in the light of other uncertainties, and that SnapWave is two orders of magnitude faster.

Figure 32 Comparison of unSWAN and SnapWave for the Haringvliet case. Top left and top right panels: Hm0 computed by unSWAN resp. SnapWave; lower left panel: bed level; lower right panel: difference between the models.

Figure 33 Scatterplot of Hm0 SnapWave against Hm0 unSWAN, for rectangular area indicated in Fig. 32.

5 Discussion

5.1 Model features

The SnapWave model uses unstructured grids based on the NetCDF ugrid convention (https://ugrid-conventions.github.io/ugrid-conventions/, last access: 21 October 2025), which can consist of a combination of triangular or quadrangular cells, for which no particular restrictions apply. The numerical method converges quickly to high accuracy (typically a relative error of $\mathrm{10}\times {\mathrm{10}}^{-\mathrm{5}}$ within 10 iterations or less), and in closed-coast cases the first sweep of the first iteration resolves most of the final solution. For the linear shoaling and refraction case we show that the model results between rectilinear and unstructured meshes are quite comparable and that the having the right resolution where the steepest gradients are governs the accuracy. In more general cases, such as determining wave fields around an island, the model allows omnidirectional propagation and refraction. The boundary conditions can be a combination of Dirichlet or Neumann conditions. The model provides convenient NetCDF output using CF conventions.

5.2 Model behaviour on open coasts and islands

In the schematic verification cases SnapWave accurately reproduces linear theory for longshore uniform coasts. For the case of a flat circular reef the model produces qualitatively similar results as the analytical solution by Mandlier and Kench (2012), with a fair match in wave height patterns, showing focusing towards the same area leeward of the shoal centre. For the case of a circular island we illustrate the rapid convergence around it and the smooth and realistic wave height pattern.

5.3 Model efficiency

All simulations were carried out on a HP ZBook Studio 16 inch G10 Mobile Workstation PC, with 13th Gen Intel(R) Core(TM) i7-13800H, 2500 Mhz, 14 Core(s), 20 Logical Processor(s) and 64 GB RAM.

The large-scale COAST3D model has approximately 340 000 net nodes, as it covers the Dutch coastal zone excluding the Wadden Sea, at a resolution down to 100 m. The extra refinements near the Coast3D site only added relatively few extra points, and the higher resolution does not influence the implicit solution in any way. The computation of one wave field took 1.9 s on average; for the 12 d of the simulation at hourly intervals this took 9 min.

The local model has approximately 5000 nodes; computation of one wave field took 16 ms, and the 12 d period at half-hour intervals took 9 s.

Per node, directional bin and condition, the large-scale model took 0.31 µm microseconds, whereas the local model only needed 0.17 µs. This difference can be attributed to the fact that the more complex large-scale model typically took 10 iterations to fully converge, where the local model typically took 6.

In Table 3 the run times and model characteristics are shown for all field validation cases. In general we can conclude that the model takes around 0.15 µs per node, directional bin and wave condition for the simplest rectangular grids, and around 0.3 µs for more complex, unstructured grids. The St Croix model is an outlier with 0.6 µs, which cannot be explained by its convergence characteristics, which are very similar to e.g. the Ningaloo large-scale grid.

The Haringvliet case compares the performance between SnapWave and unSWAN, and shows that on the exact same horizontal grid, the SnapWave model is two orders of magnitude faster while giving, for this shallow nearshore case with complex bathymetry, very comparable results.

Table 3 Overview of run time characteristics for all field validation models.

5.4 Method to transform wave conditions from ERA5 to nearshore

ERA5 performed well in all cases, with absolute bias typically less than 10 %, scatter index 20 %–25 %. Extreme events may be underestimated where ERA5 cannot resolve the atmospheric scale of the depression, as in the case of the US Virgin Islands. Results for nearshore locations have similar relative bias (∼ 10 %) and somewhat higher scatter index (∼ 30 %). The case of Ningaloo Reef poses a severe challenge because of the high friction losses, represented by a uniform friction coefficient of 0.6, and because of its sensitivity to the water level, where even a small phase error leads to large deviations in water level and hence shallow water wave heights. In this case the averaged scatter index for points on the reef is around 40 % for the large model forced by ERA5, against around 30 % for the purely local model. The effect of neglecting wave setup in the case of Ningaloo Reef was significant but small compared to this scatter. As this is the case study that is most sensitive to wave setup, we may conclude in general that neglecting wave setup in predicting nearshore wave conditions is acceptable in view of other uncertainties.

5.5 Limitations

The SnapWave model considers directionally spread waves with a single representative frequency, which introduces errors for multi-peaked spectra. Particularly the use of the peak frequency as the characteristic frequency can introduce large fluctuations in the modelled wave celerity and group speed, if swell and local wind waves compete for dominance of the spectrum. The use of a characteristic period $T_{\mathrm{m}-\mathrm{1},\mathrm{0}}$ that best represents the mean group velocity is then recommended:

$$\begin{array}{}\text{(19)}& T_{\mathrm{m}-\mathrm{1},\mathrm{0}}={\displaystyle \frac{\int S{f}^{-\mathrm{1}}\mathrm{d}f}{\int S\mathrm{d}f}}\end{array}$$

where S is the spectral density (m² Hz⁻¹) and f is the frequency (Hz). There is a wide range of bi- or multi-model wave conditions that could occur and the current model does not represent this well; this is likely one of the causes for the scatter in the nearshore comparisons against data. However, we see possible improvements in future versions to better deal with this.

First, we could improve the input functionality by reading in 2D spectra and converting these to 1D directional spectra and the Tm10 wave period. This would already improve the prediction for systems with spectral partitions from different directions but similar frequencies. Alternatively, such 1D directional spectra could be generated from integral parameters of the sea and swell bands. Both are relatively minor implementation issues that would not affect the overall solution method.

Second, we could apply different characteristic frequencies per directional bin. This would allow us to represent swell and sea from different directions and with clearly different frequencies. This is worth investigating but less straightforward, particularly in combination with wind growth.

The functionality described here does not include wave growth by wind, although this process has been implemented and is currently being tested. The model is stationary and is therefore suited for swell propagation and wave propagation over limited distances, as is typically the dominant situation in coastal areas. It can provide a fast alternative to more complex models such as SWAN when the dominant processes are wave shoaling, refraction and dissipation by friction and depth-limited wave breaking.

6 Conclusions

The SnapWave model presented here provides an efficient way to propagate wave conditions from the ERA5 hindcast, or similar global wave hind- or forecasts, to the nearshore. We have shown that the model correctly simulates nearshore wave propagation and dissipation for directionally spread waves specified at points typically 50–100 km offshore. For the cases we tested the ERA5 hindcast provides adequate boundary conditions and the combination of ERA5 and SnapWave is able to reproduce time series of wave heights at nearshore locations with significant skill. Although we have only tested the method in a few locations, we believe this approach can be used on open coasts anywhere and has the potential to be used as part of large-scale to global assessments that rely on nearshore wave conditions.

Appendix A: Error metrics

Table A1 Definitions of error metrics.

Table A2 Error metrics shoaling and refraction test, Hm0.

Table A3 Error metrics shoaling and refraction test, wave direction.

Table A4 Error statistics Coast3D local model, gamma = 0.70.

Table A5 Error statistics Coast3D large-scale model, gamma = 0.7.

Table A6 Error statistics Coast3D large-scale model, gamma = 0.75.

Table A7 Error statistics Ameland Inlet, Hm0.

Table A8 Error statistics Ameland Inlet, Tp.

Table A9 Error statistics St Croix, Hm0.

Table A10 Error statistics Ningaloo, local model.

Table A11 Error statistics Ningaloo, large-scale model, uncorrected water levels.

Table A12 Error statistics Ningaloo, large-scale model, corrected water levels.