Wave Shape Evolution from a Phase-Averaged Spectral Model

: In spectral wave models, the nonlinear triad source term accounts for the transfer of energy to the bound higher harmonics. This paper presents an extension to commonly used spectral models that resolves the evolution of the bound wave energy by keeping track of the energy that has been bound by the triad interactions. This extension is referred to as the bound wave evolution (BWE) model. From this, the spatial evolution of the bound wave height is obtained, which serves as a proxy for the nonlinear wave shape. The accuracy of these bound wave heights, and thus wave shape predictions, is highly dependent on the accuracy of the triad source term. Therefore, in this study, the capability of the LTA and SPB triad formulations to capture the growth of the bound wave height is evaluated. For both of these formulations, it is found that slope dependent calibration parameters are required. Overall, despite being computationally more expensive, the SPB method proves to be significantly more accurate in predicting the bound wave evolution. In the shoaling zone, where the bound wave energy is dominated by triads, the BWE model is well capable of predicting the nonlinear wave’s shape. In the surf zone, however, where a combination of triads and wave breaking control the spectral evolution, the BWE model over-predicts the bound wave height. Nevertheless, this paper shows the promising capabilities of spectral models to predict the nonlinear wave shape.


Introduction
In deep oceanic waters, sea-swell surface waves can generally be described by Gaussian statistics where the wave field consists of independent harmonic waves, provided that the waves are not too steep [1].When the waves propagate towards the coast, nonlinear triad wave interactions occur under the influence of decreasing water depth and variable ambient currents [2,3].This changes the initially harmonic wave shape into a nonlinear wave shape due to the presence of bound waves accompanying the freely propagating primary waves [4,5].The nonlinear wave shape ranges from skewed waves with steeper crests and flatter troughs to asymmetric waves where the wave front has pitched forward, creating a saw-tooth wave shape at breaking.This deformation of the wave shape is also present in the near-bed velocity signal, generally resulting in a net-onshore transport due to both skewness [6][7][8][9] and asymmetry [10,11].This wave-shape induced sediment transport can be the dominant transport mechanism under certain conditions, and is instrumental in properly predicting the evolution of bars [12,13] and tidal deltas [14,15].
Morphodynamic process-based models typically rely on phase-averaged spectral models for prediction of the wave field.In such models, the wave shape is not explicitly resolved and thus needs to be parameterized before calculating wave-induced sediment transport.For instance, in the morphodynamic model Delft3D [16], the approach by Isobe and Horikawa [17] was initially adopted, which estimates that the bound wave and associated non-linear wave shape using a combination of fifth order Stokes and third order cnoidal wave theory.More recently, the empirically derived wave shape parametrization of Ruessink et al. [18] was implemented, leading to improved long-term prediction of coastal morphodynamics [19].Ruessink et al. [18]'s parametrization predicts the wave shape as a function of the Ursell number, which is a nonlinearity parameter based on the local wave height, wave period and depth.This means, in particular, that the effect of the beach slope [20,21] and the spectral shape [21] are not accounted for, while they are known to influence the development of higher harmonics and, thus, the wave shape.More generally, such a parametrization intrinsically assumes that, at any point of the domain, the waves have reached some equilibrium state, such that the wave shape can be accurately described as a function of local hydrodynamic variables only.This was proven to be a limitation when considering the wave transformation over a rapidly changing bathymetry [22,23].
In this paper, we present an alternative approach to estimate the wave shape from phase-averaged spectral models.To do so, we use the fact that the sea-swell wave shape can be estimated when the bound sea-swell wave energy is known [23], and propose a method to estimate the bound sea-swell wave energy from standard spectral model outputs.More specifically, a bound wave evolution (BWE) equation is added to an existing spectral wave model SWAN [24][25][26].The extended model (SWAN+BWE) not only resolves the evolution of the total wave spectrum, but also keeps track of the evolution of the bound super-harmonic wave spectrum.This is performed through the introduction of an energy balance equation for the bound super-harmonic spectrum that relies on existing source term formulations for triad wave-wave interactions and breaking-induced dissipation.From this, the spatial evolution of the bound wave height is computed which, in turn, is used to estimate the wave shape using an analytical expression, thereby removing the need for a local parameterization.
Section 2 introduces the spectrum and bispectrum, and elaborates on how the nonlinear wave shape, the bound wave spectrum, and resulting bound wave height can be estimated from those.The analytical expression linking the bound wave height to the nonlinear wave shape, key to the proposed approach, is furthermore presented.The bound wave height evolution equation and its practical implementation is discussed in Section 3.This includes a detailed description of two existing source term parametrizations for the triad interactions, for which the key calibration parameters are identified.The model is then evaluated by comparing the predicted bound wave height and wave shape to those obtained from a wave-resolving model, SWASH [27], for a range of offshore peak periods and beach slopes.The model set-up and test simulations are explained in Section 4, while the model performance for different triad formulations is presented in Section 5. Section 6 provides a discussion on the uncertainties associated with the modeling, and Section 7 gives the conclusions.

The Spectrum, Bispectrum and Wave Shape
Spectral changes in variance, the triad source term and the wave shape are all related to each other via the bispectrum.This section describes how the spectrum and bispectrum are defined, and how the wave shape and bound wave height can be obtained from those.

The Spectrum
The single-sided variance density spectrum is calculated as: in which E[...] denotes the expected value and C m and C * m the complex amplitude and its complex conjugate at the discrete frequency f m = m∆ f , with ∆ f being the frequency resolution, obtained by applying a Fourier transform to the surface elevation.N is the number of discrete positive frequencies from the Fourier transform.
From this, the sea-swell spectral moments can be calculated as: in which m j is the j-th order spectral moment, and in which i min and i max are the indices defining the lower and upper bound of the sea-swell frequency range.The spectral moments are used to obtain the significant wave height and spectral period estimates:

The Bispectrum
The bispectral density is here presented in double-sided form as: in which p = m + n.From the bispectrum, the biphase can be obtained as: in which ℜ and ℑ denote the real and the imaginary parts.

Wave Shape and Bound Wave Height
From the spectrum and bispectrum, the normalized sea-swell wave shape [18], which represents the combination of skewness Sk and asymmetry As, is given by: in which i b,min and i b,max are the indices referring to the lower and upper frequency limits of the bound super-harmonic frequencies.
Following de Wit et al. [23], based on the work of Herbers et al. [28], the total amount of bound variance associated with super harmonics can be obtained using: in which the factor 4 originates from the single-sided variance densities used in this equation, instead of double-sided variance densities.Analogous to the wave height, the bound wave height can be obtained from the bound variance density: The normalized wave shape can also be expressed as a function of the normalized bound wave height [23]: in which proportionality factor Ψ is obtained after substitution of Equations (3a), ( 6) and (8) in Equation ( 9): The fraction on the right hand side of Equation ( 10) equals 1 (and, thus, Ψ = 3) if the sums are evaluated over the full frequency range.However, if only part of the integrals are included (for instance, in our case, by only including sea-swell frequencies), this fraction is smaller than 1.In shallow water, where relatively more energy is present in the infragravity frequency band and in the high-frequency tail, this fraction decreases.The formulation for Ψ in Equation ( 10) is, as such, consistent with the empirically derived values (Ψ between 2.70 and 2.80) found by de Wit et al. [23] using data collected from three different field sites.

Spectral Evolution Equations for Total and Bound Variance Density
This section describes how information extracted from SWAN is used to solve the evolution equation for the bound wave variance density, from which the bound wave height and, ultimately, the wave shape is calculated (Equations ( 8)-( 10)).As SWAN does not distinguish between free and bound waves, the default evolution equation is referred to as the total wave evolution equation.The evolution equation for the bound wave variance density is the novel part introduced in this paper, which allows us to keep track of the evolution of the bound wave height and wave shape.For reading purposes, the equations are presented here in a 1D form, consistent with unidirectional normally incident waves propagating over an alongshore uniform profile.

Total and Bound Wave Energy Balance
The 1D evolution equation for the total variance density in stationary conditions, as computed by SWAN, is: where E( f , x) represents the variance density as a function of frequency f and cross-shore distance x, and in which c g is the group celerity, which is computed using the linear dispersion relationship.This implicitly assumes that, in such models, all energy is freely propagating.S nl3 and S break are source terms accounting for the effect of nonlinear triad wave interactions and wave breaking, respectively.Other commonly used source terms (wind, quadruplet, white-capping and friction) are ignored in this study to simplify the analysis.Furthermore, their effects are expected to be negligible compared to S nl3 and S break in the intermediate water depth and over the limited propagation distance investigated in this study.
In the following, a separate balance is introduced for the bound energy.The S nl3 source term accounts for energy transfers due to nonlinear triad wave interactions, which are responsible for the growth of bound super-harmonics.It is positive at frequencies towards which the (bound) energy is transferred and negative at the frequencies that this energy originates from.For the bound energy equation, only the positive part of S nl3 is considered as a source term.Furthermore, a wave breaking source term S b,break is added to account for the dissipation of the bound wave energy due to depth-induced breaking.Subsequently, the evolution equation for the propagation of bound variance density in steady conditions reads: in which E b ( f , x) is the bound variance density.For consistency with the evolution equation for the total variance (Equation ( 11)), it is also assumed here that the bound variance propagates at the free wave group celerity c g .By numerically solving Equation (12), the spatial evolution of E b is obtained (see Section 3.4 for more detail), from which H b and, ultimately, S are estimated.The wave shape estimates are thus directly dependent on the accuracy of the source terms representing the triad wave-wave interactions and wave breaking.These source terms and their parametrizations are described in Sections 3.2 and 3.3.

Nonlinear Triad Energy Transfers S nl3
Different formulations are available for the S nl3 source term in spectral models [29][30][31][32], which are all derived from the following equation relating energy transfers due to triad interactions to the imaginary part of the bispectrum: in which W m,n is the nonlinear interaction coefficient for two interacting components with frequencies f m and f n .The first term on the right hand side of Equation ( 13) accounts for the sum interactions (energy exchanges between components f m , f n and f p = f m + f n ), whereas the second term accounts for the difference interactions (between f m , f n and In this section, we discuss two commonly used formulations that are both implemented in SWAN: the LTA method [29] and the SPB method [30,33].Both methods use the same nonlinear interaction coefficient derived by Madsen and Sørensen [34], based on the Boussinesq wave theory [35]: in which g is the gravitational acceleration, c and k are the wave celerity and wave number according to linear wave theory, with the subscripts indicating the frequency for which these are calculated, and d the local water depth.These two formulations differ in the number of nonlinear interactions that are accounted for, and in the way the imaginary part of the bispectrum is estimated.More specifically, the bispectral estimate needed in Equation ( 13) is derived from the evolution equation for the bispectrum, which itself depends on the trispectrum.Estimating the bispectrum, therefore, requires a closure approximation in which the trispectrum is expressed in terms of lower order spectra, hence the spectrum and bispectrum.

LTA Method
The Lumped Triad Approximation (LTA) is a computationally efficient method which assumes that the energy transfer to a certain frequency as a result of many interactions can be represented by the self-self interaction only.LTA further relies on the quasi-normal closure approximation [36], in which the trispectrum is expressed in terms of the variance spectrum only.Based on these two assumptions, an approximation is found for the absolute value of the bispectrum, which is multiplied by a parameterized form of the biphase to obtain an estimate for the imaginary part of the bispectrum [29]: in which α LTA is a calibration coefficient, and Q consists of variance cross-products: β is a parameterized form of the biphase, based on the spectral Ursell number: with This parametrization ensures a smooth transition from β = 0 (no energy transfer) for low Ur to a maximum value of β = −π/2 for high Ur-numbers.Ur crit is a calibration coefficient that controls how fast β evolves.A wide variety of Ur crit -values is found in the literature, e.g., 0.20 by Eldeberky and Battjes [37] and 0.63 by Doering and Bowen [38].Similarly, many variations for α LTA are found in the literature, (e.g., 0.05 [39], 0.1 [40], 0.25 [25], 0.5 [41] 1.0 [29] and, most recently, 0.87 [31]).In the SWAN version applied in this study (version 41.31), the default values are Ur crit = 0.2 and α LTA = 0.87.

SPB Method
An alternative approach is the Stochastic Parameterized Boussinesq (SPB) method.In contrast to the LTA method, the SPB method takes all co-linear sum and difference interactions into account, and relies on Holloway's closure approximation [42], which assumes that the trispectrum is expressed in terms of the variance spectrum and the bispectrum.The resulting expression is: with α SPB is a proportionality factor, which is 1 for unidirectional waves and lower then 1 for directional waves, to compensate for not all interactions being co-linear.∆k m,p−m = k p − k m − k p−m is the wave number mismatch, in which the individual wave numbers are computed with the linear dispersion relationship as being freely propagating components, and K is a calibration factor, with dimension L −1 .
The connection between Equations ( 13) and ( 19) relies on the fact that, for the SPB method, the imaginary part of the bispectrum is expressed as: In the SPB formulation, K is commonly taken as a linear function of the spectral peak wave number k peak : where a and b are calibration constants.Becq-Girard et al. [30] presented a form in which K = 0.95k peak,o f f shore − 0.75, based on a single laboratory calibration study [43].Because it is unclear where the offshore peak wave number k peak,o f f shore should be defined for field cases, and to prevent negative values of K that could arise as b < 0, Salmon et al. [31] proposed using a different expression: K = 0.95k peak .This is presently the default setting in SWAN.

Energy Conservation Correction for the SPB Method
The expression for SPB (Equation ( 19) is given by Becq-Girard et al. [30], and is, as such, implemented in the default SWAN version [31].As discussed by both Eldeberky [29] and Becq et al. [33], a triad source term function based on the Holloway [42] closure approximation is not energy (and energy flux) conservative.Therefore, it can lead to an artificial decay or gain of energy.Although the difference between the frequency integrated positive and negative fluxes is expected to be small, over a longer distance, it can lead to significant changes in energy flux.Therefore, a correction factor, as was suggested by Eldeberky [29], is applied in this study.This factor reduces the positive In the following, the results of the SPB method include this correction factor, even when referring to default SPB settings.The importance of including the energy conservation correction is shown in Appendix A.

Wave Breaking S break
Several source term formulations are available for the dissipation due to wave breaking.Here, the well-known Battjes and Janssen [44] model is applied, in which a constant breaker index γ defines the ratio of the maximum wave height over depth: in which α BJ is a calibration constant taken equal to 1 in the following.Q b is the fraction of breaking waves.The amount of wave breaking for the bound wave height is assumed to be directly related to the total wave breaking source term: in which α break defines how much of the breaking source term affects the bound harmonics and, thus, indirectly, also how much it affects the free harmonics.Here, it is assumed that: such that there is no preference, and the bound harmonics dissipate at the same rate as the free harmonics.

Numerical Implementation
In this study, the spectral wave model SWAN is used to solve the energy evolution equation (Equation ( 11)) and the source terms.Subsequently, an offline post-processing routine is used to compute the bound energy evolution equation (Equation ( 12)).Therefore, there is only one-way feedback: the SWAN model does influence the bound energy evolution equation, but not the other way around.The bound evolution equation for a given f is numerically resolved using a first order upwind scheme: with ∆x as the constant grid spacing used in the SWAN computations.This post-processing routine is referred to as the bound wave evolution (BWE) model.

Test Simulations and Model Set-Up
The performance of the BWE model is quantified by comparing the predicted bound wave height and wave shape to those obtained from the phase-resolving model SWASH [27].
SWASH is a time-domain model, which solves the nonlinear shallow water equations, including the non-hydrostatic pressure.It has been shown to capture short wave propagation and dispersion [27], wave-breaking dissipation [45], nonlinear wave dynamics [46], waveinduced currents [47-49] and infragravity wave dynamics [50][51][52].Of specific interest to this study is the work by [46], who convincingly showed the capability of SWASH to predict the nonlinear wave shape parameters over a planar beach, skewness and asymmetry, as well as the spectral evolution and the triad source term.More recently, ref. [53] showed, for two different laboratory experiments with a planar beach, that SWASH serves as a good benchmark model for nonlinear wave studies.Therefore, SWASH is deemed sufficiently accurate to be used as a benchmark in this study.
The comparisons focus on 1D wave propagation over a constantly sloping beach.For these simulations, the JONSWAP wave spectra are imposed at the offshore boundary with an incoming significant wave height of 1 m and a varying peak period.As the bound wave height evolution is known to be affected by the incoming wave period and the bed slope, the model's performance is evaluated for four different peak wave periods (6, 8, 10, and 12 s) and four different bed slopes (1/20, 1/50, 1/100, and 1/200).The reason to focus on variations in wave period and bed slope, rather than incoming wave height, is that these are expected to predominantly determine the uncertainty in nonlinear wave evolution predictions by spectral triad formulations.This is because, in the literature and all triad formulations, there is agreement that the triad energy transfer is a function of the energies of the interacting components (e.g., [29]).The more complicated and, thus, uncertain part is the performance of the interaction coefficients, tuning parameters, and which interactions to be included.
These inter-model comparisons are used to investigate what the optimal calibration parameters for the S nl3 terms are for the considered parameter space in both the LTA and SPB formulations.While S nl3 -calibration studies have been performed before, they did not result in general guidelines, making it challenging to choose appropriate parameter values from the large range found in the literature (see Sections 3.2.1 and 3.2.2).Furthermore, the SPB-method has only been calibrated for laboratory scale experiments [31,43], limiting the validity of the calibration for the dimensional parameter K.
The performance of the SWAN simulations is quantified with the root mean squared errors (RMSE) of the significant wave height H, bound wave height H b and spectral wave period T m02 .For the evaluation of the S nl3 source term performance, the RMSEs are computed for the region between the offshore boundary and the mean breakpoint location, defined as the point where the wave height reaches its maximum, to exclude the region where S b,break significantly affects the bound wave height prediction.This region is referred to as the shoaling zone in the remainder of this paper, whereas the region beyond the mean breakpoint is referred to as the surf zone.

SWASH Model Set-Up and Output Processing
The SWASH model is set-up with a horizontal computational grid resolution ∆x of 1 m, ensuring a minimum of 48 grid cells per offshore peak wave length.In the vertical, the grid is discretized with three vertical sigma layers.The offshore boundary is located in a depth of 20 m for the cases with T peak = 8, 10, and 12 s, and in a depth of 10 m for cases with T peak = 6 s.The offshore depth for the 6 s period is decreased to prevent high kd-numbers at the boundary, resulting in poor dispersive properties for the number of vertical layers considered [27].Moreover, these depths are large enough to ensure sufficient linearity at the wave maker boundary (Ur < 0.025 for all conditions at the offshore boundary).This is important, because bound super-harmonics are presently not included at the offshore boundary in SWASH.Nevertheless, Fiedler et al. [54] recently showed that ignoring the bound super-harmonics at the wave maker boundary, provided that it is deep enough, does not significantly influence the nearshore sea-swell third order statistics, such as the bound wave energy, skewness and asymmetry.The simulations are performed with the Hydrostatic Front Approximation breaking routine with default parameters, in order to improve the wave breaking dissipation in the case of a coarse vertical discretization [45].
Simulations are performed for 75 min, including 15 min of spin-up time in order to have 60 min of surface elevation output with a sampling frequency of 4 Hz, which is used to generate statistically reliable estimates for the spectrum and bispectrum.The output timeseries are divided in 71 semi-overlapping blocks of 100 s each to generate ensemble-averaged spectra and bispectra with a frequency resolution of 0.01 Hz and 142 degrees of freedom.To separate the sea-swell waves from the infragravity waves, only spectral values between f i min = f peak /2 and f i max = f N are included, with f peak being the peak frequency [55].These spectra and bispectra are subsequently used to obtain the bound super-harmonic wave height using Equations ( 7) and ( 8) with f i b,min = 1.5 f peak and f i b,max = 2.5 f peak .

SWAN and BWE Model Set-Up
A horizontal grid resolution ∆x is used, being 1 m for a 1/20 slope, 2.5 m for a 1/50 slope, 5 m for a 1/100 slope, and 10 m for a 1/200 slope, such that the depth difference between subsequent grid cells is similar for all simulations.The spectrum is discretized using 71 frequencies between 0.01 and 0.5 Hz that are logarithmically distributed.Because SWAN solves for the frequency directional spectra, 45 directional bins of width ∆θ = 4 • , varying between −90 • and 90 • , are used.Quasi-unidirectionality is achieved by applying a cos M directional distribution with M = 300 to the input spectrum.In order to compare simulations with exactly the same offshore spectrum, the SWASH spectrum at x = x boundary is imposed at the boundary for the SWAN simulations.Breaker parameter γ is optimized per simulation, such that the decay of H from the SWAN simulations closely resembles the decay of H in SWASH in the surf zone.This results in γ ranging between 0.52 and 0.81, depending on the bed slope and the wave period.Surface elevation spectra and S nl3 source term outputs are generated at all computational points and used to compute the bound wave energy spectrum and wave height using Equations ( 8) and ( 12) with f i b,min = 1.5 f peak and f i b,max = 2.5 f peak .
The variations in the S nl3 calibration parameters considered in this study (α LTA and Ur crit for LTA and a and b to compute K for SPB) are given in Table 1.As there is no physical reason to add a dimensional offset b in the parameterization of K (Equation ( 22)), it is here chosen to follow the approach proposed by Salmon et al. [31], i.e., to set b = 0 and only vary a, such that K is a function of the local peak wave number only.The final range of values considered is such that it includes values presented in previous calibration studies (see Table 1).In total, 182 variations of the calibration parameters are considered for the SWAN simulations using the LTA source term, and 46 for those using the SPB source term.

Results
The evolution of the (total) significant sea-swell wave height is accurately predicted both using default LTA and SPB triads.Up to the break point, the RMSE H is, at most, 0.03 m for all SWAN simulations, regardless of whether LTA or SPB triads are used.As an example, the wave height is visualized for the case with T peak = 8 s in Figure 1a.Because of this accurate prediction of H, mismatches in the evolution of H b in the shoaling zone discussed in the following are assumed to be associated with the used triad formulations, rather than with a mismatch of the evolution of the primary sea-swell energy.The evolution of H b , predicted by SWASH and SWAN using both LTA and SPB triads with default calibration coefficients, is visualized in Figure 1b for simulations with T peak = 8s.Note that depth is represented on the horizontal axis in order to show simulations with a different slope in the same figure.The evolution of H b using the SPB triads (red lines in Figure 1b) closely resembles that of reference model SWASH (black lines), resulting in an average RMSE H b of 0.041 m.With the LTA triads (blue lines in Figure 1b), on the other hand, H b is initially under-predicted in deeper water, after which it overshoots and, subsequently, leads to a substantial over-prediction.This is confirmed by a substantially higher RMSE H b of 0.11 m.
A remarkable difference between the SWAN simulations and the SWASH reference simulations is the slope dependence of the predicted bound wave height evolution as function of depth.For SWASH, no clear distinction is observed between the simulations with different slopes (i.e., black lines in Figure 1b overlap).In contrast, for the SWAN simulations with both the LTA and SPB triads, a faster growth of H b is observed for mild slopes (thick lines) than for the steeper slopes (thin lines), indicating that the effect of bottom slope is not properly included in the triad source terms.

Optimized LTA Triads
The optimal LTA calibration parameters α LTA and Ur crit are first determined for each reference simulation (i.e., the combination of the peak period and the bed slope) as the parameter values that minimize the shoaling zone RMSE H b for this specific simulation (Figure 2a).Notably, the optimal Ur crit is low for all conditions, and even zero for 8 out of 16 conditions (Figure 2b).This means that the parameterized biphase leading to the best results is independent of the local Ursell number, and equals −π/2 (Equation ( 17)), corresponding to a fully imaginary bispectrum.Since the nonlinear energy transfer is proportional to the imaginary part of the bispectrum (see also sin(β)-dependence in Equation ( 15)), this maximizes the energy transfers.This, in turn, explains the low values of the proportionality coefficient α LTA (Figure 2c) needed to compensate for the overestimation of the imaginary part of the bispectrum.In essence, this means that the uncertainty of Ursellbased biphase parametrization in the form presented in Equation ( 17) is so high that it is better to use a spatially constant value.
As the optimal values are more sensitive to changes in slope than in period, only the slope-dependent optimized values for α LTA and Ur crit are retained (see Table 2).This leads to an improved match with the SWASH results, compared to using the default values (compare blue lines in Figure 1b,c).This improvement reduces the average RMSE H b from 0.109 m to 0.031 m (see Table 3).Interestingly, the RMSE for T m02 , bulk wave statistics used in past S nl3 -calibration studies appears relatively insensitive to the change in parameters (a 6.7% increase in RMSE T m02 , while RMSE H b decreases by more than 70%).This highlights the added-value of the current approach, where the performance of the triad source terms is evaluated based on their skill to predict bound energy rather than relying on bulk parameters derived from the total (bound + free) spectrum.To optimize the SPB triads, the shoaling zone RMSE H b is calculated for all simulations with varying a-values (see Equation ( 22)).As was found for the LTA optimization, the error curves are similar for simulations with the same slope, but different period (see Figure 3).Therefore, the slope-dependent optimal a-values are obtained from the wave period-averaged error curves (black dashed lines in Figure 3).Whereas, for the simulations with a slope of 1/20, a clear single minimum is observed at a = 0.45, for the other slopes, two local minima are observed in the error curves that move further apart for increasingly milder slopes (see Figure 3 and Table 2).As the global minimum in RMSE H b is obtained for the larger a-values, these will be used in the following.The significance of the secondary minimum will be discussed in Section 6.1.2.By applying the slope-dependent, optimized a-values, the prediction of H b in the shoaling zone is improved (compare the red lines in Figure 1b,c), with a reduction of the average RMSE H b from 0.041 m to 0.024 m compared with the default values (see summary of error metrics in Table 3).

Wave Shape Prediction
The wave shape is directly obtained from SWASH with Equation ( 6) and estimated from the optimized SWAN simulations with Equations ( 9) and (10).As an example the evolution of S is presented for the case with T peak = 8 s and a slope of 1/50, as this is the case with the best prediction of H b for both LTA and SPB (see Figure 4a).The better prediction of H b by SPB can be explained by the fact that it includes all of the sum and difference nonlinear triad interactions, whereas LTA solely relies on the self-self sum interactions.The development of S is accurately predicted by simulations with both triad formulations in the shoaling zone (see Figure 4b).In contrast, within the surf zone (d < 3 m), S is clearly over-predicted by both LTA and SPB formulations, even though using SPB triads results in a close match with H and H b (Figure 4a) from SWASH in the surf zone.The over-predictions are related to the differences in the spectral shapes between the SWAN and SWASH predictions (further discussed in Section 6.2) and to the over-prediction of H b in the case of the LTA triads (Figure 4a).Next, the SWAN-predicted wave shape S, using the optimized slope dependent parameters (see Table 2), are compared with the SWASH benchmark for all slopes and wave periods.Using LTA, the wave shape in the shoaling zone is under-predicted for low S, but over-predicted for high S (Figure 5a).This results in an RMSE S of 0.08 and a R 2 of 0.91.Improved results are obtained when the SPB triads are used (Figure 5b) with a reduction in the RMSE S to 0.05 and a R 2 of 0.94.In the surf zone, both the LTA and SPB triads over-predict S (Figure 5c,d) with a concurrent increase in the RMSE S to 0.46 and 0.21 for LTA and SPB, respectively.This over-prediction is partly caused by the aforementioned spectral differences between SWAN and SWASH, and their impact on Ψ.The larger mismatch for the LTA results are related to the persistent over-prediction of H b within the surfzone (see Figure 1c).6. Discussion 6.1.Uncertainties Associated with the S Nl3 Triad Formulations 6.1.1.Parametrization of the Biphase in LTA In Section 5.1.2,it is shown that, by using the optimal calibration settings for SWAN simulations with LTA triads, the biphase estimates (Equation ( 17)) are almost independent on Ursell, given the very low values of Ur crit .For SWAN simulations with a steeper slope, Ur crit is, in fact, zero (see Table 2), corresponding to an Ursell-independent biphase of − π 2 .To further investigate this Ursell-dependence, the behavior of the biphase is examined using the SWASH simulations.A reliable estimate of the biphase at ( f peak , f peak ) is obtained from the SWASH bispectra by applying an energy-weighted average of all biphases (Equation ( 5)) that satisfy 0.5 f peak < f m , f n < 1.5 f peak .
The biphase for all SWASH simulations is visualized in Figure 6.The large scatter in biphase for very low Ursell values (Ur < 0.05) is explained by the fact that, for linear waves, the biphase is expected to be randomly distributed between −π and π (because the real and imaginary part of the bispectrum are both very close to 0).For higher Ur-numbers, the SWASH-derived biphase estimates become increasingly negative with increasing values of the Ursell number.Consistent with previous research [20,21,56,57], a clear bed slope dependence is observed for the biphase from the SWASH simulations.More specifically, for a given Ursell value, the biphase increases in absolute sense for simulations with steeper slopes (Figure 6).In contrast, the Ursell-based biphase parameterizations (e.g., the ones of Eldeberky and Battjes [37] and Doering and Bowen [38]; see the dashed and solid black lines in Figure 6) are independent of the bed slope.This partly explains the wide range of values for Ur crit proposed in the literature.The SWASH results additionally suggest that the period dependence of the biphase is not accurately represented by the local Ursell number, thereby adding to the scatter (see the vertical spread in datapoints of the same color for a given Ursell-value in Figure 6).17) with Ur crit = 0.68 (solid black line, Doering and Bowen [38]) and Ur crit = 0.20 (dashed black line, Eldeberky and Battjes [37]) as a function of the Ursell number.Each color corresponds to a given bed slope, and includes the results for all peak wave periods considered.

Parametrization of K in SPB
In the SPB formulation, the imaginary part of the bispectrum ℑ(B), which ultimately controls the strength of the nonlinear triad interactions, is expressed as a function of the calibration parameter K and the spectral cross-products Q (Equation ( 21)).We have shown, in Section 5.1.3,that a slope-dependent parametrization of a, and thus of K, was needed to optimize bound wave height predictions.This explicitly introduces a slope dependence into ℑ(B).In the following, we examine to what extent the proposed slope-dependent SPB formulation is consistent with the evolution of the imaginary part of the bispectrum predicted by SWASH.
We focus more specifically on the self-interaction ( f peak , f peak ), and examine the evolution of ℑ(B( f peak , f peak )) and Q( f peak , f peak ) predicted by SWASH for a fixed peak period T peak = 8 s but varying bed slopes (Figure 7).Reliable estimates of ℑ(B) are obtained by ensemble-averaging the bispectra from 50 SWASH simulations that only differ by the random wave phases imposed at the offshore boundary.The resulting evolution of ℑ(B( f peak , f peak )) as a function of depth is presented in Figure 7a.ℑ(B( f peak , f peak )) exhibits much stronger spatial variations for steeper slopes than for milder slopes, with large differences between the slopes over most of the shoaling and surf zones (0 < d < 10 m in Figure 7a).Q( f peak , f peak ), on the other hand, only depends on the slope at the shallowest locations (d < 5 m in Figure 7b).In the SPB formulation, ℑ(B( f peak , f peak )) is assumed to be proportional to Q( f peak , f peak ) with a proportionality factor that depends on K and ∆k peak,peak only (Equation ( 21)).Because we compare simulations with the same peak period in Figure 7, the ∆k peak,peak for a given depth is the same for all simulations and, hence, only the dependence on K remains to explain differences in behavior between the different slopes.This confirms that a proper prediction of ℑ(B( f peak , f peak )), according to the SPB formulation, requires a slope-dependent K factor.This implies that K cannot be proportional to the local peak wave number only, as expressed in Equation (22).Alternative approaches proposed in the literature to parameterize K as a small constant value [29], as a function of the offshore peak wave number [30,33], or as a function of the minimum of the three interacting wave numbers [32], for instance, do not solve this issue, as none of them involve a slope dependence.
Furthermore, two local minima in RMSE H b were identified for the milder slopes (see Figure 3).Taking the simulation with T peak = 8 s and slope = 1/200 as an example, Figure 8a shows the prediction of H b for both of the local minima.It shows that, depending on the choice of a, H b is either overestimated between d = 10 m and d = 4 m and underestimated between d = 4 to d = 3 m, or the other way around.This suggests that, even for a given slope, a spatially varying value of a and, thus, K should be used which is not linearly proportional to the local k peak .

Bound Wave Prediction in the Surf Zone
Where S is accurately predicted in the shoaling zone using SPB triads, an overprediction is found in the surf zone, even for the cases where the evolution of H and H b is accurately predicted in the surf zone (see, for example, Figure 4 with T = 8 s and a slope of 1/50).This indicates that the over-prediction is caused by an over-prediction of Ψ which, itself, depends on the prediction of the spectral shape.This is confirmed by calculating Ψ for the SWASH and the SWAN spectra for this case, showing that, between d = 4 and d = 1 m, the Ψ from SWAN is, on average, 23% higher, with a maximum over estimation of 44%.The same simulation with the LTA triads also considerably overestimates Ψ (average 17% and maximum 24%).
The mismatch in Ψ using the SPB triads is examined next by looking in more detail at the spectra in the surf zone for this case (see Figure 9).It is seen that, before breaking (d = 4 m), there is a good agreement between the SWASH and SWAN spectral shape, as the secondary super harmonics and the high frequency tail are properly resolved (compare the magenta and black lines).When the water depth decreases the SWAN spectra maintain their initial shape while the total energy reduces as a result of dissipation due to breaking.In contrast, the SWASH spectrum gradually transforms into a saturated spectrum towards shallower water.So, even though the bulk energy transfer and decay due to triads and breaking is properly captured by SPB (see Figure 8b), the distribution over the frequencies is not.As Ψ is based on the spectral shape (see Equation ( 10)), this leads to errors in the estimated wave shape (Figures 4 and 5).Interestingly, the spectra predicted by SWAN using the a-value corresponding to the secondary local minimum in RMSE H b (see Table 3), become more saturated with decreasing depth (compare blue and black lines in Figure 9).From Equation (11), it can be seen that in absence of breaking the spatial flux gradient and the triad source term should be equal.Both the SWASH and SWAN simulations indicate that the influence of breaking is negligible at the surfzone edge at 4 m water depth.Therefore, the differences in the spectral evolution for the two a-values can be investigated by comparing the SWAN triad source terms with the SWASH-inferred spatial gradient in flux dEc g /dx (see Figure 10).In cases of a larger a-value, the source term is relatively narrow, with two distinct peaks moving energy from the primary spectral peak to a well-defined secondary peak at twice the primary frequency.As S nl3 is relatively weak, and dissipation by breaking in the current formulation is linearly proportional to the energy density, the spectral shape is conserved throughout the surfzone.For the smaller a-value, the non-linear transfer to the secondary super-harmonic spectrum is wider and significantly stronger.This results in a broadening and subsequent saturation of the high frequency part of the spectrum ( f > f peak ), which is more consistent with the SWASH results (Figure 9).However, this is accompanied by a significant transfer of energy to the lower frequencies ( f < f peak ), erroneously adding energy at half of the peak frequency (compare the blue and black lines at f = 0.063 Hz in Figure 10).This effect is exacerbated for milder slopes.Again, it is likely that the mismatch within the surfzone between the SPB-SWAN and SWASH spectra can be reduced by using a locally optimized a value.
Figure 10.SWAN S nl3 source term and SWASH flux gradient (dEc g /dx) for the simulation with T peak = 8 s and a slope of 1/50 at a depth of 4 m.The cyan and magenta lines represent SWAN simulations using SPB triads with a-values of 0.12 and 0.90, respectively.

Outlook
Since K is the only tunable parameter in the SPB method in this study, mismatches in H b and, thus, S are mainly associated with the choice for K. Currently, K is optimized for the shoaling zone to limit the influence of the breaking parametrization on the results, but also because this is the area where the wave shape plays a dominant role in the sediment transport [7,11,58].This approach leads to sub-optimal results for the surf zone (see Figure 8).Ideally, a generic parameterization of K performs well in both.It should be kept in mind, however, that the assumptions on the nonlinear interaction coefficients (Equation ( 14)) and bispectral estimates also contribute to this mismatch.
The fact that spatially varying and slope-dependent calibration factors seem to be required might indicate that the triad formulation itself should be improved.In the benchmark model (SWASH), the bound wave height as a function of depth (so dE b c g /dd) seems to be near-independent of the bed slope.The BWE SWAN models show a substantial bed slope dependence (see Figure 1), which can be explained because the triad term S nl3 = dE b c g /dx is slope-independent.In this study, this mismatch is reduced by using the calibration parameters, where the optimal values might be outside of the physically explainable range (e.g., a fully imaginary bispectrum for LTA and very low or high values for K in SPB for 1/200 slope).A potential future solution could be to directly include the slope dd/dx in the triad formulation, and re-calibrate tuning factor α SPB .This would reduce the triad strength for the milder slopes and increase the strength for the steeper slopes.If we revisit Figure 1), this would improve the default solutions.Subsequently, K can be re-calibrated to fine-tune the predictions, rather than using it to compensate for the missing slope dependence.
A way forward for addressing all of these aspects is by creating ensemble averages using multiple SWASH simulations with different random phases, from which the accurate stochastic estimates can be retrieved for the variance spectrum, bispectrum and even the trispectrum.From these (higher order) spectra, the assumptions on the closure approximation can be validated per triad and location.These findings can subsequently be used to obtain improved expressions for the triad formulations, the nonlinear interaction coefficients, and the imaginary part of the bispectrum and, thereby, extend the domain for which reliable predictions for the bound harmonics and the wave shape can be obtained.
Finally, it is noted that this study is restricted to simulations with a wave height of 1 m.For higher incoming waves, the wave transformation will be different.The waves will start breaking in deeper water, and the higher nonlinearity will result in a stronger evolution of the higher harmonics, and thus higher bound wave height, skewness and asymmetry.As all triad formulation are relatively straightforward related to the incoming energy, it is expected that this will be well-predicted, provided that proper slope and period dependent triad formulations are used.Nevertheless, this should still be verified in a future study.

Conclusions
In this paper, a new method is presented that allows for phase-averaged spectral wave models to predict a nonlinear wave shape.The spatial evolution of the bound variance density spectra is resolved using source term functions extracted from the spectral wave model (SWAN) that account for the influence of nonlinear energy transfers and wave breaking.By integrating the bound variance density spectra, the bound wave height is obtained.Using the proportionality factor, derived in Equations ( 9) and (10), this bound wave height can serve as a proxy for the nonlinear wave shape.As the BWE model's skill to predict the bound wave height and wave shape is predominantly determined by the accuracy of the triad source term, the performance of the LTA and SPB triad formulations are evaluated.This is performed by comparing the bound wave height predictions from the integrated bound evolution equation to those obtained from a detailed phase-resolving wave model (SWASH) using bispectral analysis.Although computationally more expensive, the SPB method proved to be significantly more accurate in predicting the bound wave height for all test simulations that cover a range of bed slopes and peak periods.The performance of the SPB source term can be further improved by optimizing calibration parameter K.This study already shows that different optimal K-values are found for different bed slopes, but more research is required to find a generally applicable formulation for K. Nevertheless, using the SPB method, the predicted wave shape agrees very well with the wave shape directly obtained from the reference model (SWASH) in the shoaling zone (R 2 = 0.96 and RMSE S = 0.07).In the surf zone, however, an over-prediction of the wave shape is observed.This arises when the proportionality factor is computed in shallow water, and because there is no source term accounting for the release of bound harmonics.Overall, this paper shows the promising capability of a wave-averaged spectral wave model to predict the nonlinear wave shape.Due to the limited additional computational effort used by the added evolution equation, this could be used in large-scale morphological models to improve the wave-shape-induced sediment transport formulations.
Author Contributions: All co-authors contributed to the developed methodology.F.d.W. performed and post-processed the numerical simulations and wrote the manuscript.M.T. and A.R. contributed to supervising, editing, and improving this paper.All authors have read and agreed to the published version of the manuscript.

Figure 1 .
Figure 1.Evolution of H and H b using SWASH (black lines) and SWAN with LTA (blue lines) and SWAN with SPB (red lines), where different bed slopes are shown with different line widths for T peak = 8s.The top panel (a) shows H results for default triad settings, the middle panel (b) shows Hb for default triad settings, whereas the bottom panel (c) shows results for optimized triad settings.5.1.Bound Wave Height Prediction 5.1.1.Default S Nl3 Settings

Figure 2 .
Figure 2. Overview of the errors and calibration parameters for SWAN LTA simulations with optimized settings for the different slopes and conditions.Panel (a) shows RMSE H b for the simulation with the lowest RMSE.Panels (b,c) show the corresponding optimal values for Ur crit and α LTA .

Figure 3 .
Figure 3. RMSE H b for simulations over a constant bed slope of 1/20 (panel a), 1/50 (panel b), 1/100 (panel c), and 1/200 (panel d) as a function of the SPB calibration parameter a. Different colors indicate simulations with different incoming peak wave periods; the black dashed line indicates the average over the four simulations with different periods.

Figure 4 .
Figure 4. Spatial evolution of the bound wave height (panel a) and normalized wave shape (panel b) using reference model SWASH (black lines), SWAN with LTA triads with optimized settings (blue lines), and SWAN with SPB triads with optimized settings (red lines).

Figure 5 .
Figure 5. SWAN-predicted normalized wave shape using LTA (a) and SPB (b) triads in the shoaling zone and LTA (c) and SPB (d) triads in the surf zone as a function of normalized wave shape from reference model SWASH.

Figure 6 .
Figure 6.Estimated biphase β( f peak , f peak ) from SWASH simulations (colored dots) and parameterized biphase using Equation (17) with Ur crit = 0.68 (solid black line, Doering and Bowen[38]) and Ur crit = 0.20 (dashed black line, Eldeberky and Battjes[37]) as a function of the Ursell number.Each color corresponds to a given bed slope, and includes the results for all peak wave periods considered.

Figure 7 .
Figure 7. Evolution of ℑ(B( f peak , f peak )) in panel (a) and Q( f peak , f peak ) in panel (b) as a function of depth for the SWASH simulation with T peak = 8 s and varying bed slopes.

Figure 8 .
Figure 8. Spatial evolution of the bound wave height for T peak = 8 s and a slope of 1/200 (panel a) and a slope of 1/50 (panel b) for reference model SWASH (black line) and SWAN with SPB triads for the a-values corresponding to the two local error minima listed in Table 2 (colored lines).

Figure 9 .
Figure9.Variance density spectra for the simulation with T peak = 8 s and a slope of 1/50 at depths of 4, 3, 2, 1 m, in panels (a-d), respectively.Black lines represent the SWASH simulations, whereas the blue and red lines represent SWAN using SPB triads with a-values of 0.12 and 0.90, respectively.

Table 1 .
Variations in conditions and calibration values applied for the SWAN calibration simulations for the S nl3 source term formulations using the LTA and SPB methods.Here, default refers to the default parameter values in SWAN version 41.31.

Table 3 .
RMSE H , RMSE H b and RMSE T m02 in the shoaling zone for LTA and SPB with default and optimized calibration parameters per slope (see Table2for their values).