Improving predictions of swash dynamics in XBeach The role of groupiness and incident-band runup

In predicting storm impacts on sandy coasts, possibly with structures, accurate runup and overtopping simulation is an important aspect. Recent investigations (Stockdon et al., 2014; Palmsten and Splinter, 2016) show that despite accurate predictions of the morphodynamics of dissipative sandy beaches, the XBeach model (Roelvink et al., 2009) does not correctly simulate the individual contributions of set-up, and infragravity and incident-band swash to the wave run-up. In this paper we describe an improved numerical scheme and a different way of simulating the propagation of directionally-spread short wave groups in XBeach to better predict the groupiness of the short waves and the resulting infragravity waves. The new approach is tested against ﬁ eld measurements from the DELILAH campaign at Duck, NC, and against video-derived runup measurements at Praia de Faro, a relatively steep sandy beach. Compared to the empirical ﬁ t by Vousdoukas et al. (2012) the XBeach model performs much better for more extreme wave conditions, which are severely underestimated by existing empirical formulations. For relatively steep beaches incident-band swash cannot be neglected and a wave-resolving simulation mode is required. Therefore in this paper we also test the non-hydrostatic, wave-resolving model within XBeach for runup and overtopping against three datasets. Results for a high-quality ﬂ ume test show non-hydrostatic XBeach predicts the run-up height with good accuracy (maximum deviation 15%). A case with a very shallow foreshore typical for the Belgian coast at Wenduine was compared against detailed measurements. Overall the model shows correct behavior for this case. Finally, the model is tested against a large number (551) of physical model tests of overtopping from the CLASH database. For relatively high overtopping discharges the non-hydrostatic XBeach performs quite well, with increasing accuracy for increasing overtopping rates. However, for relatively low overtopping rates of less than 10 – 20 l/m/s, the model systematically underestimates measured overtopping rates.

Improving predictions of swash dynamics in XBeach: The role of groupiness and incident-band runup Dano Roelvink a, b, c, *  In predicting storm impacts on sandy coasts, possibly with structures, accurate runup and overtopping simulation is an important aspect. Recent investigations (Stockdon et al., 2014;Palmsten and Splinter, 2016) show that despite accurate predictions of the morphodynamics of dissipative sandy beaches, the XBeach model  does not correctly simulate the individual contributions of set-up, and infragravity and incident-band swash to the wave run-up. In this paper we describe an improved numerical scheme and a different way of simulating the propagation of directionally-spread short wave groups in XBeach to better predict the groupiness of the short waves and the resulting infragravity waves. The new approach is tested against field measurements from the DELILAH campaign at Duck, NC, and against video-derived runup measurements at Praia de Faro, a relatively steep sandy beach. Compared to the empirical fit by Vousdoukas et al. (2012) the XBeach model performs much better for more extreme wave conditions, which are severely underestimated by existing empirical formulations. For relatively steep beaches incident-band swash cannot be neglected and a wave-resolving simulation mode is required. Therefore in this paper we also test the nonhydrostatic, wave-resolving model within XBeach for runup and overtopping against three datasets. Results for a high-quality flume test show non-hydrostatic XBeach predicts the run-up height with good accuracy (maximum deviation 15%). A case with a very shallow foreshore typical for the Belgian coast at Wenduine was compared against detailed measurements. Overall the model shows correct behavior for this case. Finally, the model is tested against a large number (551) of physical model tests of overtopping from the CLASH database. For relatively high overtopping discharges the non-hydrostatic XBeach performs quite well, with increasing accuracy for increasing overtopping rates. However, for relatively low overtopping rates of less than 10-20 l/m/s, the model systematically underestimates measured overtopping rates.

Background
XBeach ) is an open-source, process-based morphodynamic numerical model for the nearshore and coast. XBeach was originally developed to simulate the impact of extreme storms and hurricanes on sandy barrier island systems (e.g., McCall, et al., 2010;Lindemer, 2010), including beach and dune erosion, overwash, inland flooding and barrier rollover and breaching. The model development philosophy has been to explicitly resolve dominant physical processes, where possible, and to revert to empirical relations where the physical processes are either not sufficiently well understood to be explicitly modeled, or their computation would be prohibitively expensive in typical coastal engineering applications. This philosophy has led to a relatively flexible model that allows for the seamless simulation of all levels of impact, e.g., swash, collision, overwash and inundation regimes (Sallenger, 2000), across alongshore-varying coastal geometries.
While the assumption of a saturated surf zone allows a simplification of the description of the incident-band waves into an energy balance, it also limits the application of the XBeach model on steeper coasts where considerable incident-band energy may propagate into the inner surfzone and swash. This limitation led to the development of a new branch of the XBeach model for steep gravel coasts, called XBeach-G Masselink et al., 2014), which allows a phase-resolving approach for infragravity and incident-band waves using a non-hydrostatic pressure correction term for the non-linear shallow water equations (Smit et al., 2010), in a manner similar to the SWASH model Smit et al., 2014). XBeach-G was extensively validated using field data and data collected in a large-scale physical model experiment, and the model showed great skill at predicting extreme wave run-up and barrier overtopping and overwash events.

Objective
Within the Risc-Kit project (van Dongeren et al., 2015) XBeach plays a central role in translating events generated offshore to concrete hazards for a variety of coastal types, with and without coastal defense structures. Apart from dune erosion and overwashing, runup and overtopping are essential variables to predict, as part of this model chain. Therefore the objective of this paper is a) to improve the observed deficiencies in the 2D infragravity runup, by examining, improving and validating the numerical and physical schemes for short-wave energy propagation, and b) to validate the runup and overtopping for intermediate to steep profiles with the wave-resolving, non-hydrostatic model within XBeach.

Outline
In the following section we will describe in general terms the different physical representations of short waves within XBeach and indicate their application areas. In section 3 we describe improvements in the numerical scheme for the time-varying wave energy balance, which become important when propagating wave groups over large distances. We identify the causes for a strong damping of the wave groupiness in the existing physical representation for directionally-spread waves and propose a solution that not only resolves this adequately but is also two to three times faster. This is then validated against field datasets at Duck, NC and Praia de Faro, Portugal. In section 4 we compare the runup simulated by the non-hydrostatic model with a series of physical model tests for a dike behind a complex bar profile. The overtopping predicted by the model is first evaluated in detail for a dike behind a shallow foreshore. Finally, the model is tested against a large dataset of 551 physical model tests of overtopping over dike profiles of varying toe depths, slopes and crest heights. Conclusions and recommendations are given in the last section.

Model descriptionstationary, surf-beat and non-hydrostatic mode
XBeach was originally developed as a short-wave averaged but wavegroup resolving model, allowing resolving the short wave variations on the wave group scale and the long waves associated with them. Since the original paper by Roelvink et al. (2009) a number of additional model options have been implemented, thereby allowing users to choose which time-scales to resolve: Stationary wave model (efficiently solving wave-averaged equations but neglecting infragravity waves; Surf-beat mode (instationary), where the short wave variations on the wave group scale (short wave envelope) and the long waves associated with them are resolved; Non-hydrostatic mode (wave-resolving), where a combination of the non-linear shallow water equations with a pressure correction term is applied, allowing to model the propagation and decay of individual waves.
In the following these options are discussed in more detail (see Fig. 1).

Stationary mode
In stationary mode the wave-group variations and thereby all infragravity motions are neglected. This is useful for conditions where the incident waves are relatively small and/or short, and infragravity motions would be small anyway. The model equations are similar to HISWA (Holthuijsen et al., 1989) but do not include wave growth or wave period variations. Processes that are resolved are wave propagation, directional spreading, shoaling, refraction, bottom dissipation and wave breaking, and a roller model (Nairn et al., 1990;Roelvink and Reniers, 2011); these processes are usually dominant in nearshore areas of limited extent. For the breaking dissipation we use the Baldock et al. (1998) model, which is valid for wave-averaged modeling. The radiation stress gradients from the wave and roller model force the shallow water equations, drive currents and lead to wave setdown and setup.
Further details and possible applications are discussed in the XBeach User Manual (online at xbeach.readthedocs.io).

Surf beat mode (instationary)
The short-wave motion is solved using the wave action equation with time-dependent forcing. This equation solves the variation of shortwaves envelope (wave height) on the scale of wave groups. It employs a dissipation model for use with wave groups (Roelvink, 1993;Daly et al., 2012) and a roller model (Svendsen, 1984;Nairn et al., 1990;Roelvink and Reniers, 2011) to represent momentum stored at the surface after breaking. These variations, through radiation stress gradients Stewart, 1962, 1964) exert a force on the water column and drive longer period waves (infragravity waves) and unsteady currents, which are solved by the nonlinear shallow water equations (e.g. Phillips, 1977). Thus, wave-driven currents (longshore current, rip currents), long (infragravity) waves, and runup and rundown of long waves (swash) are included.
Using the surf-beat mode is necessary when the focus is on swash zone processes rather than time-averaged currents and setup. It is fully valid in the swash on dissipative beaches, where the short waves are mostly dissipated by the time they are near the shoreline. On intermediate beaches and during extreme events the swash motions are still predominantly in the infragravity band and so is the wave runup.

Non-hydrostatic mode (wave resolving)
For non-hydrostatic XBeach calculations depth-averaged flow due to waves and currents are computed using the non-linear shallow water equations, including a non-hydrostatic pressure correction. The depthaveraged normalized dynamic pressure (q) is derived in a method similar to a one-layer version of the SWASH model . The depth averaged dynamic pressure is computed from the mean of the dynamic pressure at the surface and at the bed by assuming the dynamic pressure at the surface to be zero and a linear change over depth. Full derivation and validation results are given in Smit et al. (2010).
Under these formulations dispersive behavior is added to the long wave equations and the model can be used as a short-wave resolving model in intermediate to shallow water depths (kh 2.5, where k is the wave number and h the water depth). Although wave breaking is implicitly solved in the model equations, we apply the hydrostatic front approximation of Smit et al. (2013) to improve the computed location and magnitude of breaking.
In case the non-hydrostatic mode is used, the short wave action balance is no longer required. However, in the wave-resolving mode we need much higher spatial resolution and associated smaller time steps, making this mode generally much more computationally expensive than the surf-beat mode.
The main advantages of the non-hydrostatic mode are that the incident-band (short wave) runup and overtopping are included, which is especially important on steep slopes such as gravel beaches. Another advantage is that the wave asymmetry and skewness are resolved by the model and no approximate local model or empirical formulation is required for these terms. Finally, in cases where diffraction is a dominant process, wave-resolving modeling is needed as it is neglected in the short wave averaged mode. An application of the non-hydrostatic mode is XBeach-G, which is a branch of the main XBeach source code but is specifically developed to simulate storm impacts on gravel beaches Masselink et al., 2014). The formulations for gravel beaches were developed and extensively tested for the non-hydrostatic mode. Sandy morphology can be simulated using the wave-resolving mode but has not yet been validated as extensively, though promising results are presented in literature, e.g. Daly et al. (2017).
3. Improving IG wave motions in surf-beat: role of wave groupiness 3.1. Introduction Stockdon et al. (2014) found that 2D XBeach significantly underpredicted infragravity runup for a large number of field observations at the Field Research Facility at Duck, NC, while 1D simulations tended to overpredict for the higher wave conditions. The latter may be expected as wave directionality tends to reduce the infragravity wave forcing (e.g. Herbers et al., 1994); however, the underprediction for the 2D case was unexpected and has led us to further analyze the causes.
In reviewing some of the simulations at the base of Stockdon et al. (2014) we noticed that the variation of the short wave height, which can be visualized as moving 'blobs', rapidly reduced towards the shore in the standard two-dimensional horizontal (2DH) propagation scheme. This variation is responsible for forcing bound infragravity waves (Longuet Higgins and Stewart, 1962) and is also at the base of the breakpoint mechanism (Symonds et al., 1982). The variability of the short wave energy can be expressed in the so-called Groupiness Factor GF (Funke and Mansard, 1980), which we define in our case as where H is the slowly-varying wave height, the overbar indicates timeaveraging and σ H the standard deviation of H. In idealized cases the bound infragravity waves are linearly proportional to this groupiness. Recent improvements of the code are given in Roelvink et al. (2015) and the XBeach manual. Here we limit ourselves to two improvements relevant to the cases in this paper: a new numerical scheme to solve wave group energy propagation that removes the potential for inaccurate/ unphysical numerical growth of wave groupiness, and a new method for the derivation of directional wave group energy propagation (single-dir option) that improves wave groupiness coherence.

Improved propagation scheme for wave action balance
The second-order upwind scheme implemented in XBeach has very little dissipation and is therefore very suitable for propagating the wave groupiness over large distances; however, in some cases it was reported by users that the scheme could lead to too sharp variations in wave energy, which could lead to wiggles in energy and water level. We analyzed this behavior by propagating long-crested waves over a uniform depth of 10 m and comparing the time series of wave height at the boundary with that after 7 wave lengths L, shifted in time by 7L/c g , where c g is the group velocity. Ideally these time series should be on top of each other; what we saw is that for simulations with fine grid resolutions the higher-frequency variations in wave height tended to grow, leading to steeper wave height variations and frequent occurrences of zero wave height. The results can be expressed in terms of the mean wave height and GF.
In Fig. 2 we show the results for the second-order upwind scheme, which indicate that the mean wave height is propagated without error but the variation, expressed in GF, tends to grow for small grid sizes and damp for larger grid sizes; the optimum lies around grid sizes typically recommended, in the order of 15-20 m in 10 m water depth. Note that the spectral significant wave height H m0 indicated in the figure is a time- To overcome the undesired effects of steepening of wave groups we implemented a correction to the second-order upwind scheme according to Warming and Beam (1976), which implies a small additional diffusion term which is a function of time step and group velocity. This gives results as indicated in Fig. 3. The behavior is now always slightly damped, with an acceptable dissipation of wave groupiness using typical grid sizes, which converges to zero for decreasing grid sizes. This behavior is much preferred over that where the amplitude error increases for decreasing grid sizes, and the improved scheme will be used in the following sections.

Multi-dir vs. single-dir
The standard mode of solving the time-varying wave energy balance in XBeach is to compute the propagation of wave energy (or action) in x, y and θ space simultaneously, by solving the 3D advection equation: where D w is the dissipation by wave breaking, D f by bottom friction and D v by vegetation. Here the wave action A is calculated as: Aðx; y; t; θÞ ¼ S w ðx; y; t; θÞ σðx; y; tÞ (3) where θ represents the angle of incidence with respect to the computational x-axis, S w represents the wave energy density in each directional bin and σ the intrinsic wave frequency. The intrinsic frequency σ and group velocity is obtained from the linear dispersion relation. For each directional bin i the horizontal propagation speeds are equal to: The refraction of the waves is produced by the 'refraction speed' c θ . The way the directional wave groups are propagated leads to excessive smoothing of the longshore wave groupiness: the wave energy from different directional bins is simply added up, without considering the interference of the different wave components.
An alternative to this approach, which we named the 'single-dir' option, is to first calculate the mean wave directions and to propagate the short wave energy along these directions. This can be achieved by alternating a stationary run to obtain the mean wave direction with instationary runs where the following reduced equation is solved: Fig. 2. Hm0 wave height and Groupiness Factor (GF) after 7 short wave lengths, relative to their values at the boundary, as a function of relative grid size; second order upwind scheme. Fig. 3. Hm0 wave height and Groupiness Factor (GF) after 7 short wave lengths, relative to their values at the boundary, as a function of relative grid size; corrected second order upwind scheme according to Warming and Beam (1976).
D. Roelvink et al. Coastal Engineering 134 (2018) 103-123 This effectively reduces the problem from a 3D to a 2D problem, where the occasional (say, every 10 min) stationary run described in 2.1 to obtain the direction takes relatively little time. It is comparable to the approach taken in the 'roller model' in Delft3D (e.g., Reniers et al., 2004), which used intermittent SWAN runs to predict the mean wave direction, which was then used to solve the time-varying wave action and roller energy balances.

Testcase propagation over uniform depth
To illustrate the effect of the single-dir option on wave groupiness, we carried out tests for a model domain of 1200 m longshore by 1000 m cross-shore and a constant depth of 10 m with a JONSWAP wave spectrum with mean direction 20 from shore normal, H m0 wave height of 2 m, peak period T p of 10 s and directional spreading coefficient s of 10 (equivalent to directional spreading of 24 ; typical sea conditions) and 64 (10 ; swell), respectively. We compare non-hydrostatic, wave-resolving mode ('nonh'), standard surf-beat mode ('multi-dir') and the new singledirection mode ('single-dir'). In all cases cyclic boundary conditions were used to allow obliquely incident waves to leave the southern boundary and to continue on the northern boundary. This avoids any disturbances at the lateral boundaries. Since kh, the product of wave number k and water depth h, is 0.68 for this case, well within the validity range of the non-hydrostatic model (kh < 2), we consider the non-hydrostatic model as the benchmark. In order to compare the wave-resolving with the surfbeat mode, we computed the slowly-varying wave height in the nonhydrostatic mode by low-pass filtering the squared surface elevation with a cutoff at 0.5 times the peak frequency as is usual; to obtain the infragravity surface elevation we applied the same filter to the surface elevation itself. Only short test runs of 300 s were carried out to illustrate the behavior; the first 3 min were omitted from the analysis to allow the waves to pass through the domain.
In Fig. 4 we compare the wave height fields and long wave fields for three situations: 1. Non-hydrostatic mode; we have filtered the results to obtain the short wave height field and the water level fields; 2. Single-dir mode, the new implementation; 3. Multi-dir, with a directional bin size of 10 deg.
From the figure we can see that the single-dir wave height pattern retains much more variability of the wave height than the multi-dir standard version, in comparison with the non-hydrostatic result. The patterns between non-hydrostatic and single-dir are very comparable. The effect on the infragravity waves is less pronounced but still clearly present. In Fig. 5 we present the same results for a directional spreading coefficient s ¼ 64. Clearly, more longshore groupiness is preserved in this case, but again the 'single-dir' option shows patterns much more similar to the benchmark non-hydrostatic result.
In Figs. 4 and 5, bottom panels, the longshore-averaged GF is shown for the three simulation modes, for directional spreading s ¼ 10 and s ¼ 64, respectively. Over approximately 6-7 wavelengths, the GF in the Fig. 4. Snapshot of slowly-varying wave height (top panels) and infragravity surface elevation (middle panels) for wave-resolving, non-hydrostatic mode (left panels), single-directional mode (middle panels) and multi-directional mode (right panels); bottom panels: longshore-averaged groupiness factor (GF) as a function of distance from sea boundary, for uniform depth of 10 m and directional spreading coefficient s ¼ 10. non-hydrostatic model remains more or less constant; the 'multi-dir' simulation shows a reduction of the GF to a third in the case of s ¼ 10, and approx. two thirds in the case of s ¼ 64; the 'single-dir' simulation shows a decrease in the order of 15% over the first 10 wave lengths for the case of s ¼ 10 and 5% for the case of s ¼ 64.

Area of interest
In order to verify the 2DH hydrodynamics of XBeach when forced by directionally-spread short waves, a simulation was set up to compare model results to field measurements. In this case the DELILAH field experiment at Duck, North Carolina was selected as a suitable test location. The period that is modeled is October 13th, 1990, which was a stormy day, between 16:00 and 17:00 h. This period was chosen earlier in Roelvink et al. (2009) andVan Dongeren (2003) as a case representative of energetic conditions; though we could have run the model throughout the measurement campaign we believe the single experiment gives a good insight into the ability of the model to describe the crossshore transformation of wave and flow parameters and spectra. The significant wave height at 8 m water depth was 1.81 m, with a peak period of 10.8 s and a mean angle of incidence of À16 relative to the shoreward normal. This period was selected because the wave conditions are energetic enough to generate a significant infragravity wave component and the incident wave spectrum is sufficiently narrowbanded (comparable to a JONSWAP spectrum with s ¼ 6) to justify the model assumption of a narrow-banded frequency spectrum. The model is forced with the actual wave spectrum measured at 8 m water depth (Birkemeier et al., 1997). A measured tidal signal is imposed on the model boundaries of which the mean level is 0.69 m above datum (see Fig. 6).

Model set-up
The model was directly taken from the XBeach model skill bed (Deltares, 2012) and has been presented in Roelvink et al. (2009). A uniform grid size of 5 m cross-shore by 10 m longshore was applied, with a horizontal extent of 850 m cross-shore by 700 m longshore. The Warming and Beam propagation scheme was used for all simulations. To avoid disturbances at the lateral boundaries, cyclic boundary conditions were applied. Though testing of breaker formulations in itself was not an objective of this paper, we ran simulations with two different breaker formulations to assess the robustness of our conclusions: The default Roelvink (1993) The Roelvink-Daly formulation (Daly et al., 2012), with its default gamma ¼ 0.52 ('RoelvinkDaly') For both of these we compared multi-dir with 15 deg directional bins with single-dir simulations. The simulation time was 3800 s of which the first 200 were ignored in the analysis. Time series with 1 s intervals were created at the measurement locations.

Analysis
The time series from the model were treated in the following manner to allow accurate comparison with the observations provided by Birkemeier et al. (1997). The long wave water level variance in XBeach, though starting at the boundary with no energy in frequency components Fig. 5. Snapshot of slowly-varying wave height (top panels) and infragravity surface elevation (bottom panels) for wave-resolving, non-hydrostatic mode (left panels), single-directional mode (middle panels) and multi-directional mode (right panels); bottom panels: longshore-averaged groupiness factor (GF) as a function of distance from sea boundary, for uniform depth of 10 m and directional spreading coefficient s ¼ 64.  above 0.5 times the peak frequency f p , develops higher components due to shoaling and steepening of the infragravity waves. Therefore we filtered the water surface elevation time series with a frequency filter into high-frequency components above 0.5 f p (z s,hi ) and low-frequency components below 0.5 f p (z s,lo ). We then computed high-frequency and lowfrequency root-mean-square wave height H rms according to: The longshore velocity v was time-averaged. Simulated time series (3600 s at 1 Hz) were detrended and converted to power spectra, applying a Hanning filter. Measured spectra were available averaged over 0.0078 Hz bins. We therefore averaged the simulated highresolution power spectra over the same bins for direct comparison.

Results
In Fig. 7 the resulting wave height distributions over the profile are shown for the default wave breaking formulations, and for both multi-dir and single-dir options. As could be expected, the mean Hrms HF distributions do not differ significantly, since it is mainly the groupiness that can be expected to be different. Nevertheless there is an effect of the increased groupiness in that the (non-linear) dissipation is higher during the passage of maxima in the wave height, which leads to slightly lower mean wave height for single-dir. The LF wave heights are very well reproduced for the single-dir simulation, whereas they are slightly underestimated for the multi-dir option; the difference is easily explained by the reduced groupiness of the HF waves.
The longshore velocity shown in the middle panel misses importantly the observed peak in current velocity in the trough in both cases, though the single-dir option comes closer.
Mainly to address the discrepancy in the longshore velocity we examined the effect of using the breaker formulation according to Daly et al. (2012), which is similar to that by Dally (1992) in that the breaking of waves is tracked along with the propagation of the wave energy, and waves that have started breaking only stop doing so after their height gets below a second threshold (typically 0.3 times the water depth). This leads to more dissipation and hence forcing of the longshore current beyond the bar crest. Fig. 8 show the same results as before for this formulation. The HF wave height, for the default settings with this formulation, is now slightly higher in both cases, and the single-dir results match the shoaling region better, while the multi-dir results slightly overestimate it. For the LF wave height, the single-dir results are slightly on the high side whereas the multi-dir Hrms LF are a little underestimated. As expected, the longshore velocity is reproduced much better, especially in the singledir option.
The LF spectra are compared with the observations in Fig. 9 and Fig. 10. In general there is a good agreement in location of the peaks and in the quantitative reproduction of the power spectrum (values presented in absolute spectral variance density). For the case of the Roelvink93 formulations the single-dir option clearly outperforms the multi-dir simulation, with the exception of the most offshore point. For the Roel-vinkDaly formulation the results are more mixed, as a result of the light overprediction of the LF energy in the single-dir simulation, but especially in the nearshore area the agreement is still good.
The performance statistics for these tests are shown in Table 1. We present correlation coefficient R 2 , Scatter Index (SCI) and Relative Bias. The following conclusions can be drawn: The RoelvinkDaly formulation performs better for the Hrms HF , both in single-dir and multi-dir.
For the LF wave height single-dir in combination with the Roelvink93 scores best and much better than multi-dir for that formulation. The performance for multi-dir somewhat surprisingly is quite good in combination with the RoelvinkDaly formulation. Possibly a somewhat higher friction coefficient for the long wave motion could lead to quantitatively better results for single-dir; this has not been investigated further here.
For the longshore velocity RoelvinkDaly and single-dir perform clearly better; their combination gives quite acceptable agreement, though the model still underestimates the magnitude of the velocity in the trough. The spectra are represented very well by the Roelvink93, single-dir combination, with the lowest SCI and bias, though the combination RoelvinkDaly, single-dir gives the highest correlation.  The run-up was measured using a video system from which time stacks of water lines were extracted; with the help of 40 accurate beach surveys the water line horizontal positions could be converted to vertical levels and the R 2% run-up level could be determined. Simultaneous wave height, period and direction information was obtained from a nearby wave buoy and water levels from a nearby tidal station. The details of the measurement procedures are given in Vousdoukas et al. (2012). We obtained from Dr Vousdoukas a relevant selection of 301 data points where all this information was synchronized and quality checked. An important parameter included in this series was the slope of the intertidal beach. Vousdoukas et al. (2012) investigated several possible empirical parameterizations starting from Stockdon et al. (2006); we have used his optimum formulation as a reference to evaluate the XBeach results.

Setup of XBeach simulations
An automated Matlab procedure was set up for this validation. All input data, observations and Matlab scripts to generate the results are stored on the XBeach repository, under folder testcases/Vousdouka-s2012_Praia_de_Faro. We created grids with cell sizes varying from 10 m offshore to 1 m in the swash zone. Based on the measured intertidal beach slopes the profiles were adjusted in this region, as illustrated in Fig. 11.
All cases were run in XBeach 1D (both non-hydrostatic and surf-beat mode) and in 2DH (surf-beat mode).
The beaches at Praia de Faro have a relatively steep upper beach slope, following an approximately 1:30 slope from À15 m to À3 m. Under such circumstances the assumption that most of the swash zone energy is in the infragravity band is questionable, as there will be some energy in the incident band. Therefore for the 1D runs we compared both with the non-hydrostatic and surf-beat mode.
The 1D runs were carried out without directional spreading and assuming the waves to be perpendicular to the beach. For the 2DH runs we could only carry out the simulations in surf-beat mode, as these runs would have required a very high resolution and hence computation time.

Results
In Fig. 12 the wave conditions, water levels, beach slopes and resulting R 2% run-up height computed with the non-hydrostatic mode simulations are shown and compared with the observed R 2% and those predicted by Vousdoukas et al. (2012), eq. (7): Here H 0 is the deep water H m0 wave height, L 0 the deep water wave length and ξ the Iribarren parameter.
The comparison of the series of measurements, the empirical predictions and the non-hydrostatic model results shows a significant overestimation of the run-up heights for this case, very likely due to the absence of directional spreading and because the waves were modeled as perpendicular to the beach, whereas the average deep water angle relative to the coast orientation of 218 w.r.t. N is 28 .
The scatter plots for this situation, shown in Fig. 15, show a trend line for the simulations against the predictions with a slope that is 28% too steep if the trend line is forced through zero; if the zero intercept is free, the slope of the trend line is only overestimated by 6%.
By comparison, the empirical relationship slightly under predicts the trend in case the regression line is forced through zero; when the zero intercept is let free, we see that the regression line has a slope of 0.51, thus severely underestimated; this is particularly important for extreme conditions.
For the 1D surf-beat mode simulations the results are shown in Fig. 13  and Fig. 16. A similar pattern can be seen for this case, where the trend line for the XBeach results against the observations now has a slope that is 15% too high, both when the line is forced through zero and for the case of a free zero intercept.
The 2DH runs were carried out in single-dir mode. The main purpose of running in 2DH mode was twofold: 1. To take into account effects of directional spreading 2. To take into account the effect of wave refraction Various experiments were carried out to find an optimum alongshore grid resolution and extent. The result of this was that a longshore grid size of 50 m and alongshore extent of 1000 m sufficiently captured the typical size and shape of the wave groups so further refinement or extension in longshore direction was not necessary.
The results for these simulations are shown in Fig. 14 and Fig. 17. Clearly, the simulated R 2% values much more closely match the observed ones. As the scatter plots show, though considerable scatter remains, the  Fig. 11. Base, mildest and steepest profile for simulations.
D. Roelvink et al. Coastal Engineering 134 (2018) 103-123 trend is quite accurate, both in the regression line forced through zero and the one with a free zero intercept.

Introduction
In the following we compare the runup simulated by the nonhydrostatic XBeach model with a series of physical model tests on a dike behind a complex bar profile. The overtopping predicted by the model is first evaluated in detail for a dike behind a shallow foreshore. Finally, the model is tested against a large dataset of 551 physical model tests of overtopping over dike profiles of varying toe depths, slopes and crest heights.

Introduction
Within the framework of the European MAST-OPTICREST project, prototype measurements were performed on the Petten sea defense in the Netherlands (van Gent, 1999(van Gent, , 2001. The main characteristic of this dike is a complex shallow foreshore, which makes it rather difficult to apply standard runup formulas; the non-hydrostatic XBeach model should be able to resolve the wave transformation over the shallow foreshore and the runup on the steep dike profile, but this had not been validated yet. There are two types of measurements for these tests, of which wellcontrolled physical model tests have been used to see the response of the XBeach for different wave conditions. The physical model tests were performed in the Scheldt Flume of Deltares in Delft.    D. Roelvink et al. Coastal Engineering 134 (2018) Fig. 18 shows the profile of the foreshore in the last kilometer, which could be modeled in the flume, and the detailed profile of the dike. Wave run-up levels were measured relative to SWL at the upper slope by sensors, acting as a step-gauge within the smooth slope. Also in the physical model tests, the measurements were performed on this upper slope. It should be noted that the minimum water layer thickness on the upper slope is considered as 0.10 m for these prototype tests. The crest elevation is NAPþ12.9 m and all slopes are considered to be smooth.
JONSWAP type wave conditions were specified using parametric spectra defined case by case. Physical models of Petten cases were analyzed based on specified peak period and significant wave height. The directional spreading coefficient s is set to 1000 (unidirectional waves). All data are given and all models were run at prototype scale.  D. Roelvink et al. Coastal Engineering 134 (2018) 103-123 A uniform grid with resolution of 1 m was applied. The simulation was run over a period of 2 h to obtain a good statistical representation of the R 2% run-up height.

Results
The conditions for these tests are listed in Table 2. The R 2% run-up height as measured in the physical model and predicted by XBeach is shown in Fig. 19. The results show that for a typical 1D (cross-shore) application with a complex shallow foreshore and dike with berm, under controlled conditions with second-order steering and reflection compensation, non-hydrostatic XBeach predicts the run-up height with good accuracy, in these cases with a maximum deviation of 15%.

Introduction
In this section, a detailed case was modeled with XBeach to investigate its capacity to simulate wave overtopping and run-up for the Flemish coastal town of Wenduine in Belgium. This town has been highlighted as one of the weak links in the Flemish master plan for coastal safety.
The model in this section is based on a study that has been done to test the capability of the SWASH model  for wave transformation and overtopping discharge for four selected impermeable dikes with different characteristics (Suzuki et al., 2011). The SWASH model was validated based on the comparison to the physical model measurements in terms of significant wave height, spectral analysis, maximum wave height, wave set-up, average period and peak period for one of the four cases. Also, wave overtopping discharge was obtained  from SWASH and physical models for all the tests. Fig. 20 shows the geometry of this physical model including the positions of six wave gauges starting from offshore to the toe of the structure. The topography and the dike, consisting of a shallow foreshore 1:35 and a relatively steep slope dike (1:2) represents the typical configuration of the coastline at Wenduine, Belgium. Table 3 indicates the detailed characteristics of the first physical model test. This case is a very specific case and falls outside standard coastal structure design rules. The water level at the toe of the dike is only a few decimeters, whereas the initial significant wave height in deep water is almost 5 m. The main difference with "conventional" coastal structure design is this very low water depth. A rule of thumb gives that the depth limited wave height is roughly half the water depth. That would mean one or 2 dm. Reality is much different, due to the transformation of a short wave spectrum to a spectrum which has predominantly long waves. Note however that the long wave phenomena in this situation will likely be overestimated compared to reality since the short-crestedness of waves (2D-effect) will reduce the long wave energy.

Case description
The physical model topography and sea dike were constructed at a Froude scale of 1:25. The topography was simplified into a 1/35 foreshore slope starting 13.3 m from the wave paddle up to the toe of the sea-dike, and constructed from smooth concrete. The dike has a wide crest with a seaward 1/100 slope. In this section one test was conducted with irregular waves (Jonswap γ ¼ 3.3) with one storm condition and dike condition shown in Table 3.
In the physical model, the wave flume length is 70 m, width is 4 m and height is 1.45 m. Fig. 20 shows the model geometry, with the vertical axis at prototype scale while the horizontal axis is at model scale.
Numerical modeling of the mentioned case has been carried out by XBeach model using a grid size of 0.04 m in the horizontal. The geometry mentioned in Fig. 20 was reproduced in the numerical domain at physical model scale. Given a selected grid size of 4 cm, the length of the numerical flumes was 52 m long with 1300 grid cells. The time duration of the numerical simulation was 40 min, the same as the physical model experiment. Non-hydrostatic 1D boundary conditions and absorbinggenerating (weakly-reflective) back boundary condition in 1D were applied for the front and back side of the numerical model, respectively.

Results
In Fig. 21 through Fig. 23 the measured wave spectra are compared with those simulated in XBeach. Obviously, there is a difference in the exact spectral shape at the offshore point, where XBeach shows a typical JONSWAP shape and the model is less peaked. This could have been overcome by imposing the measured time series or spectra, but apparently, the model is not very sensitive to the exact spectral shape since in shallow water (points 5 and 6) the spectra are dominated by the infragravity waves and the shape and energy in these low-frequency spectra are reproduced well by the model.
The H m0 evolution is also reproduced nicely as is shown in Fig. 24. The most sensitive parameter here is maxbrsteep: the criterion determining when the non-hydrostatic correction is turned off and wave dissipation inherent in the nonlinear shallow water equations takes over for breaking waves (cf. HFA of Smit et al., 2013). As shown in Fig. 24, the optimum value for this case is 0.4, but the deviation for the default value of 0.6 is relatively small.
In the same Fig. 24, as could also be seen in the spectra, the peak period jumps to values in the range 20-40 s, indicating the dominance of infragravity waves in the very shallow area.

Overtopping rate
The instantaneous discharge qx was monitored at the end of the crest    Fig. 25 a snapshot is shown of the water level during one of the overtopping events.
In Table 4 the mean overtopping discharge is given for different values of the roughness and maxbrsteep parameter. Simulated overtopping rates using a maxbrsteep value of 0.4 are close to the observed value of 0.58 l/s/m; for the default setting of 0.6 the values are somewhat overestimated. The results are insensitive to the roughness setting.
Overall however, the model shows correct behavior for this case and even when default settings are applied the overtopping discharge is overestimated by a factor 1.5, which can be considered acceptable in view of the large scatter usually found in the measurements.

Introduction
In many of the RISC-KIT case studies the problem area has a shallow foreshore; however, this is not always the case, especially during storm surge conditions. Therefore it is useful to test the applicability of XBeach for representing overtopping over dikes and breakwaters without a particularly shallow foreshore.
To this end, a collection of 551 relatively simple cases was taken from the CLASH database (Steendam et al., 2004), of which 366 points from very recent tests by Victor and Troch (2012) with accurate second-order wave generation and active reflection compensation. Comparison of XBeach with this large number of data points provides an insight in the predictive capacity of the model and into the bias and scatter.

CLASH database
The CLASH database consists of an excel sheet with over 10,000 data points described by a name, 8 wave parameters (H m0 ,T p , mean wave period T m and spectral wave period T m-1,0 at deep water and at the toe of dike or breakwater), 18 parameters defining the structure, two factors defining the reliability and the complexity, and the measured mean overtopping rate. With the help of Prof. J.W. van der Meer 511 reliable and simple (smooth surface, single slope, deep foreshore as defined in Van Gent (1999): H m0,deep /toe depth < 0.4) cases were selected for this validation.    D. Roelvink et al. Coastal Engineering 134 (2018) 103-123 4.4.3. Setup of the simulations An automated Matlab procedure was set up for this validation; details of the procedure can be found in Roelvink et al. (2015). An important aspect here is that the CLASH database contains the incoming H m0 wave height only; to ensure that the model simulation had the same incoming wave height, the incoming and outgoing waves were separated in the model results and where necessary the model was rerun with adjusted input wave heights.
With the help of the parameter depthscale, an input parameter in XBeach that allows to scale all criteria such as the drying/flooding criterion eps, we attempted to avoid all numerical scale effects. We set the depthscale at 20 m divided by the water depth near the wave maker.
The grid size varies: in deep water 50 cells per wave length are used and in shallow water and on the slope the grid size equals the offshore depth divided by 50; for small-scale tests in 50 cm this means the horizontal resolution on the slope is 1 cm; for Delta Flume conditions of approx. 5 m depth this is still a fine 10 cm resolution. As for other settings, all were taken at default values except for the roughness, where a Manning's n of 0.01 was applied and the viscosity, which was turned off completely.
Separating the incoming and reflected waves in the numerical results was done using time series of water level and velocity, the short wave celerity and the local water depth following Guza et al. (1984): where z s is the water level, u the velocity, c the wave celerity and h the water depth. The incoming wave height produced by XBeach was typically in the order of 10% below the target value; adjusting for this led to improved results. An example of the results of such an analysis is given in Fig. 26. The incident wave height from this analysis at half a wavelength from the boundary was made to match the observed H m0,deep from the CLASH database.

Results
The mean overtopping discharges for both measurements and simulations were converted to the same prototype scale, taking a 20 m offshore depth as prototype scale. This allows an easy comparison of the different tests related to the prototype overtopping discharge, regardless of the scale the tests were carried out at.
In Fig. 27 the simulated discharge rates are plotted against the measured ones for all tests. Clearly, there is a good correspondence for the higher overtopping rates, while the model performance suffered below values of 10-20 l/m/s. Part of the scatter may be explained by the fact that older tests are included, which for instance do not have reflection compensation. However, as the results for the cases from Victor and Troch (2012) show in Fig. 28, there is still a systematic underestimation for low overtopping rates, though, importantly, there are hardly any 'false negatives' in these cases and the overall performance is much better.
By the colouring of the data points with the ratio crest height A c over H m0 wave height we clearly see both that the overtopping strongly increases with decreasing A c /H m0 ratio but also that scatter and bias increase with increasing A c /H m0 . We can use this to obtain meaningful statistics of the error as a function of A c /H m0 ; we do this by computing the SCI and relative bias per bin of A c /H m0 of 0.5. The results are shown in Fig. 29. Apart from an unexpected peak in SCI for values for Ac/Hm0 between 1 and 1.5, there is a clear trend for the scatter to increase with this ratio and for the bias to become more negative; at A c /H m0 around 2, for instance, the overtopping is systematically underpredicted by a factor of two. Since Ac/Hm0 is firmly known one could think of correcting the predicted values for the bias; this is a subject of further study.
These results are qualitatively comparable with the uncorrected results presented in Suzuki et al. (2017), in that they also showed, in their Fig. 7, a similar pattern of increasing scatter and underestimation for lower overtopping discharges. Their cases concerned dikes with shallow or very shallow foreshores, and include the Wenduine case of which we present one case in 4.3, with good results. It is quite possible that for numerical model predictions, the 'difficult' cases with complex shallow foreshores are in fact well suited, since both SWASH and XBeach-nonhydrostatic have been shown to predict wave transformation over shallow profiles and reefs accurately; the overtopping discharges in these cases are dominated by low-frequency waves, the generation of which is well predicted by the models; the breaker type of the low-frequency waves themselves will generally be surging or non-breaking, which makes it relatively easy to simulate overtopping. For the deep foreshore cases presented here the short-wave breaking processes on the dike, which are often of collapsing or plunging type, are not captured as accurately given the inherent problems in representing such processes in a one-layer model. Obvious next steps will be a) to compare wave propagation, breaking and individual overtopping events for deep foreshore cases in detail, for which the CLASH database does not provide the data, and b) to run through the entire set of CLASH datasets in order to cover the full spectrum from deep to very shallow foreshores.
We may conclude that for relatively high overtopping discharges the non-hydrostatic XBeach performs quite well; recent, high-quality data are reproduced for 91% within a factor 10 and 85% within a factor 2. These rates improve for increasing overtopping rates.
However, for relatively low overtopping rates of less than 10-20%, the model systematically underestimates the overtopping rates. This will be the subject of further studies. There is a consistent trend in the relative bias and slope with A c /H m0 , which could possibly be employed to correct the simulated overtopping rates.

Conclusions
For short-wave averaged, 'surf-beat' simulations, an improved numerical scheme and a different way of simulating the propagation of directionally-spread short wave groups results in better prediction of the groupiness of the short waves and the resulting infragravity waves. The new approach consists of intermittently computing the mean wave direction using the stationary solver within XBeach, and then propagating the wave energy along these directions. Apart from being more accurate the approach is also 2-3 times faster in typical applications.
The new approach was tested against field measurements from the DELILAH campaign at Duck, NC. We found that the combination of an improved numerical scheme, the breaking model according to Daly et al. (2012) and the single-dir option gives the best overall performance, with accurate reproduction of HF and LF wave heights, longshore currents and LF spectra. This test was mainly aimed at validating the correct short-and long wave transformation through the surf zone, for which we focused on one particular case with storm waves.
In a test against over 300 video-derived runup measurements at Praia de Faro, the 2D single-dir model gave the best prediction for runup, Fig. 27. Comparison of measured and simulated mean overtopping discharges for all selected cases; results scaled up to prototype scale. The colors indicate the ratio crest height A c over H m0 wave height. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 28. Comparison of measured and simulated mean overtopping discharges for cases in Victor and Troch (2012); results scaled up to prototype scale. The colors indicate the ratio crest height A c over H m0 wave height. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) compared to 1D simulations in non-hydrostatic and surf-beat mode which overestimated the runup by 30% and 13% respectively. The slope in the trend was approx. 0.9, indicating a slight underestimation due to the omission of incident-band swash. Compared to the empirical fit by Vousdoukas et al. (2012) the XBeach model performs much better for more extreme wave conditions, which are severely underestimated by existing empirical formulations. The tests show that for cases with only partially saturated surf zones the surf-beat approach, with the new singledir option, gives only a slight underprediction of the runup by approx. 10%, and effects of directional spreading and oblique incidence have to be accounted for to avoid overestimating runup.
The non-hydrostatic, wave-resolving model within XBeach was tested for runup and overtopping against three datasets. Results for the Petten, the Netherlands case show that for a typical 1D (cross-shore) application with a complex shallow foreshore and dike with berm, under controlled conditions with second-order steering and reflection compensation, nonhydrostatic XBeach predicts the run-up height with good accuracy, in these cases with a maximum deviation of 15%.
A case with a very shallow foreshore typical for the Belgian coast at Wenduine was compared against detailed measurements. Overall the model shows correct behavior for this case; the transformation of wave spectra towards very shallow water is predicted correctly and the overtopping rate is predicted accurately for calibrated parameter settings that reflect the correct wave height decay; when default settings are applied the overtopping discharge is overestimated by a factor 1.5, which can be considered acceptable in view of the large scatter usually found in the measurements. The bed roughness only has a limited effect on the overtopping rate.
Finally, the model was tested against a large number (551) of physical model tests of overtopping from the CLASH database. For relatively high overtopping discharges the non-hydrostatic XBeach performs quite well; recent, high-quality data are reproduced for 91% within a factor 10 and 85% within a factor 2. These rates improve for increasing overtopping rates. However, for relatively low overtopping rates of less than 10-20%, the model systematically underestimates the overtopping rates. This will be the subject of further studies. There is a consistent trend in the relative bias and slope with the freeboard to wave height ratio, which could possibly be employed to correct the simulated overtopping rates.