3.1 Incorporation of first-order spatial derivatives and dissipation effects into the nonlinear source term

This subsection focuses on the inclusion of spatial inhomogeneities into deterministic and stochastic models. For conciseness a 1-D formulation is used in the derivation. In all former deterministic and stochastic models (e.g. Kaihatu & Kirby Reference Kaihatu and Kirby1995; Agnon & Sheremet Reference Agnon and Sheremet1997; Vrecica & Toledo Reference Vrecica and Toledo2016) the nonlinear interaction term ( $W$ ) was defined in terms of wavenumbers and angular frequencies, and the first and higher order spatial derivatives ( $k^{^{\prime }}$ ) were discarded. As the wave field enters shallower waters, slopes grow larger and the discarded terms can become significant. These effects can be included in the governing equations by following the same derivation procedure as shown in Vrecica & Toledo (Reference Vrecica and Toledo2016), with the exceptions presented below. This again results in (2.1) and (2.2), but with a new formulation of the $W$ -terms.

The derivation procedure of Bredmose et al. (Reference Bredmose, Agnon, Madsen and Schaffer2005) is followed up to the combined bottom boundary condition. A linear dispersion relation is defined as

(3.1) $$\begin{eqnarray}\displaystyle \left(\tan (h\unicode[STIX]{x1D735})\unicode[STIX]{x1D735}+{\displaystyle \frac{\unicode[STIX]{x1D714}_{p}^{2}}{g}}\right)\hat{\unicode[STIX]{x1D719}}_{p}=0, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D735}$ is differential operator on $\hat{\unicode[STIX]{x1D719}}_{p}=a_{p}\text{e}^{-\int k_{p}\,\text{d}x}$ , which is a Fourier amplitude of velocity potential. Equation (3.1) is solved to yield a known linear dispersion relation, and inserted into the full nonlinear equation. Here, a different approach is followed. The shoaling term and the imaginary part of the nonlinear term are retained, and a linear dissipation term ( $\unicode[STIX]{x1D6FC}_{p}$ ) is added to the equation to yield

(3.2) $$\begin{eqnarray}\displaystyle \left(\tan (h\unicode[STIX]{x1D735})\unicode[STIX]{x1D735}+{\displaystyle \frac{\unicode[STIX]{x1D714}_{p}^{2}}{g}}\right)\hat{\unicode[STIX]{x1D719}}_{p} & = & \displaystyle -\sec (h\unicode[STIX]{x1D735})(\unicode[STIX]{x1D735}h)\boldsymbol{\cdot }\left\{(\cos (h\unicode[STIX]{x1D735})\unicode[STIX]{x1D735}\hat{\unicode[STIX]{x1D719}}_{p})-\sin (h\unicode[STIX]{x1D735}){\displaystyle \frac{\unicode[STIX]{x1D714}_{p}^{2}}{g}}\hat{\unicode[STIX]{x1D719}}_{p}\right\}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{p}\hat{\unicode[STIX]{x1D719}}_{p}-\unicode[STIX]{x1D700}\overset{N}{\underset{s=p-N}{\sum }}F_{s,p-s}^{(2)}\hat{\unicode[STIX]{x1D719}}_{s}\hat{\unicode[STIX]{x1D719}}_{p-s}+O(\unicode[STIX]{x1D700}^{2}).\end{eqnarray}$$

The real part of the right-hand side (only the nonlinear term has a significant real component) modifies the solution of the dispersion relation. This can account for bound wave properties, and gives second-order correction of the dispersion relation for most energetic harmonics. However, this has limited effect on the overall evolution, and as it requires iterative computation it is not further considered. Therefore, the solution of (3.1) is used for the real part of the wavenumber. The imaginary terms on the right-hand side act as a forcing, causing the wavenumber to become complex. It can be defined as

(3.3) $$\begin{eqnarray}\displaystyle \tilde{k}_{p}=k_{p}+\text{i}\bar{k}_{p}, & & \displaystyle\end{eqnarray}$$

where $\bar{k}_{p}$ is imaginary component, and $\tilde{k}_{p}$ is the complex wavenumber. The $\bar{k}_{p}$ is obtained from solving

(3.4) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}(a_{s}(x)\text{e}^{-\text{i}\int k_{s}\,\text{d}x})=-\text{i}k_{s}a_{s}(x)\text{e}^{-\text{i}\int k_{s}\,\text{d}x}+a_{s}^{^{\prime }}(x)\text{e}^{-\text{i}\int k_{s}\,\text{d}x}, & & \displaystyle\end{eqnarray}$$

as

(3.5) $$\begin{eqnarray}\displaystyle \bar{k}_{p}=-{\displaystyle \frac{|a_{s}^{\prime }(x)|}{|a_{s}(x)|}}=-G_{s}. & & \displaystyle\end{eqnarray}$$

In some cases fast oscillations of the $G_{s}$ -term can cause singularities in the calculation, therefore it is recommended to smoothen out the $G_{s}$ -term by evaluating the change of amplitude over longer distances. A finite difference method is used for estimating the  $\bar{k}_{p}$ .

The change of wave amplitudes is then described through an imaginary part of the phase function ( $a_{p}(x)\text{e}^{\text{i}\int k\,\text{d}x}=a_{p}\text{e}^{\text{i}\int \tilde{k}\,\text{d}x}$ ). For example, solving the dispersion relation which accounts for the change of bathymetry, results in the imaginary component which exactly corresponds to the solution of the mild slope equation with a bed slope term. Dissipation and change due to nonlinear terms behave in a similar manner.

The equation (3.2) is then multiplied with the $H$ -operator defined as

(3.6) $$\begin{eqnarray}\displaystyle H={\displaystyle \frac{\unicode[STIX]{x1D735}+\text{i}k_{p}}{\tan (h\unicode[STIX]{x1D735})\unicode[STIX]{x1D735}+{\displaystyle \frac{\unicode[STIX]{x1D714}_{p}^{2}}{g}}}} & & \displaystyle\end{eqnarray}$$

to yield

(3.7) $$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D735}+\text{i}k_{p})\hat{\unicode[STIX]{x1D719}}_{p}=D_{p,l}\hat{\unicode[STIX]{x1D719}}_{p}-{\displaystyle \frac{\unicode[STIX]{x2202}C_{g,p}/\unicode[STIX]{x2202}x}{2C_{g,p}}}\hat{\unicode[STIX]{x1D719}}_{p}-\unicode[STIX]{x1D700}\overset{N}{\underset{s=p-N}{\sum }}HY_{s,p-s}\hat{\unicode[STIX]{x1D719}}_{s}\hat{\unicode[STIX]{x1D719}}_{p-s}. & & \displaystyle\end{eqnarray}$$

It is assumed that the wave field does not vary in time and lateral direction, the shoaling term reduces to $(\unicode[STIX]{x2202}C_{g,p}/\unicode[STIX]{x2202}x)/(2C_{g,p})$ , and $D_{p}=H\unicode[STIX]{x1D6FC}_{p}$ .

The inclusion of an imaginary wavenumber component has little effect on the imaginary component of the shoaling term, based on preliminary simulations. However, the $W$ and $Y$ terms are much more sensitive to the inclusion of imaginary components. The effect of amplitude slopes can be much more significant than that of four wave interactions in certain cases, such as strong nonlinearity.

In previous works, the first-order spatial derivatives of wavenumbers and bottom were neglected in the nonlinear part. However, the inclusion of real components of derivatives again results in additional imaginary components, and as the term is sensitive to them they are retained in this paper. The $Y$ -term is defined as

(3.8) $$\begin{eqnarray}\displaystyle Y_{s,p-s} & = & \displaystyle {\displaystyle \frac{\text{i}}{2g^{3}}}\unicode[STIX]{x1D714}_{s}^{2}\unicode[STIX]{x1D714}_{p-s}^{2}\unicode[STIX]{x1D714}_{p}-{\displaystyle \frac{\text{i}}{2g^{3}}}\unicode[STIX]{x1D714}_{s}\unicode[STIX]{x1D714}_{p-s}\unicode[STIX]{x1D714}_{p}^{3}\nonumber\\ \displaystyle & & \displaystyle -\,{\displaystyle \frac{\text{i}}{g}}\unicode[STIX]{x1D714}_{p}\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{p-s}-{\displaystyle \frac{\text{i}}{2g}}\unicode[STIX]{x1D714}_{p-s}\unicode[STIX]{x1D6FB}_{s}^{2}-{\displaystyle \frac{\text{i}}{2g}}\unicode[STIX]{x1D714}_{s}\unicode[STIX]{x1D6FB}_{p-s}^{2}.\end{eqnarray}$$

See Bredmose et al. (Reference Bredmose, Agnon, Madsen and Schaffer2005) and Vrecica & Toledo (Reference Vrecica and Toledo2016) for details. If $O(h^{\prime })$ are retained in the $Y$ -term, and if $\unicode[STIX]{x1D735}$ -terms are substituted with the solution of (3.2), the $Y$ -term is defined as

(3.9) $$\begin{eqnarray}\displaystyle {\tilde{Y}}_{s,p-s,u,l-u}^{(2)} & = & \displaystyle {\displaystyle \frac{\text{i}}{2g^{3}}}\unicode[STIX]{x1D714}_{s}^{2}\unicode[STIX]{x1D714}_{p-s}^{2}\unicode[STIX]{x1D714}_{p}-{\displaystyle \frac{\text{i}}{2g^{3}}}\unicode[STIX]{x1D714}_{s}\unicode[STIX]{x1D714}_{p-s}\unicode[STIX]{x1D714}_{p}^{3}+{\displaystyle \frac{\text{i}}{g}}\unicode[STIX]{x1D714}_{p}(\tilde{k}_{s,u}^{x}\tilde{k}_{p-s,l-u}^{x}+k_{s,u}^{y}k_{p-s,l-u}^{y})\nonumber\\ \displaystyle & & \displaystyle +\,{\displaystyle \frac{\text{i}}{2g}}\unicode[STIX]{x1D714}_{p-s}\tilde{k}_{s,u}^{2}+{\displaystyle \frac{\text{i}}{2g}}\unicode[STIX]{x1D714}_{s}\tilde{k}_{p-s,l-u}^{2}-{\displaystyle \frac{1}{2g}}\unicode[STIX]{x1D714}_{p-s}k_{s,u}^{^{\prime }}-{\displaystyle \frac{1}{2g}}\unicode[STIX]{x1D714}_{s}k_{p-s,l-u}^{^{\prime }}.\end{eqnarray}$$

The $H$ -term operates both on the $Y$ -term and amplitudes. The $Y$ -term is first expressed in an exponential form as

(3.10) $$\begin{eqnarray}\displaystyle Y_{s,p-s}(x)=Y_{s,p-s}(0)\text{e}^{-\text{i}\int _{0}^{x}V_{s,p-s}\,\text{d}x^{\prime }}. & & \displaystyle\end{eqnarray}$$

The $V$ -term is composed of a real and an imaginary part. The real part of the $V$ -term is the derivative of phase of the $Y$ -term and represents the variation of the ratio of real and imaginary components of the $Y$ -term. The imaginary part is the local growth of the $Y$ -term. It is defined as

(3.11) $$\begin{eqnarray}\displaystyle V_{s,p-s}={\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}\ln \left({\displaystyle \frac{Y_{s,p-s}(x)}{Y_{s,p-s}(0)}}\right)\text{i}={\displaystyle \frac{Y_{s,p-s}^{^{\prime }}(x)}{Y_{s,p-s}(x)}}\text{i}. & & \displaystyle\end{eqnarray}$$

The above equation can be easily derived by solving (3.10). In this way the $V$ -term fully describes changes of the $Y$ -term. The new formulation of the $H$ -term is then given as

(3.12) $$\begin{eqnarray}\displaystyle \tilde{H}_{s,p-s}={\displaystyle \frac{\text{i}k_{p,l}^{x}+V_{s,p-s}-\text{i}\tilde{k}_{s,u}-\text{i}\tilde{k}_{p-s,l-u}}{-\text{i}(-V_{s,p-s}+\tilde{k}_{s,u}+\tilde{k}_{p-s,l-u})\tan (-\text{i}(-V_{s,p-s}+\tilde{k}_{s,u}+\tilde{k}_{p-s,l-u})h)+\unicode[STIX]{x1D714}_{p}^{2}/g}}.\qquad & & \displaystyle\end{eqnarray}$$

From here, the equations are redefined in terms of surface elevation amplitudes following the approach established in Bredmose et al. (Reference Bredmose, Agnon, Madsen and Schaffer2005) and Vrecica & Toledo (Reference Vrecica and Toledo2016), and are not given in detail. The only deviation in the procedure is that the first-order spatial derivatives are retained, to finally yield the nonlinear interaction term $W$ as

(3.13) $$\begin{eqnarray}\displaystyle \tilde{W}_{s,p-s}=-{\displaystyle \frac{\unicode[STIX]{x1D714}_{p}}{g}}(\tilde{k}_{s,u}^{x}+\tilde{k}_{p-s,l-u}^{x}-k_{p,l}^{x})\tilde{T}_{s,p-s}^{(2)}+g{\displaystyle \frac{\unicode[STIX]{x1D714}_{p}}{\unicode[STIX]{x1D714}_{s}\unicode[STIX]{x1D714}_{p-s}}}\tilde{H}_{s,p-s}{\tilde{Y}}_{s,p-s}^{(2)}, & & \displaystyle\end{eqnarray}$$

with

(3.14) $$\begin{eqnarray}\displaystyle \tilde{T}_{s,p-s}^{(2)}={\displaystyle \frac{g^{2}\tilde{k}_{s}\tilde{k}_{p-s}}{2\unicode[STIX]{x1D714}_{p}\unicode[STIX]{x1D714}_{s}\unicode[STIX]{x1D714}_{p-s}}}-{\displaystyle \frac{1}{2}}\unicode[STIX]{x1D714}_{p}+{\displaystyle \frac{\unicode[STIX]{x1D714}_{s}\unicode[STIX]{x1D714}_{p-s}}{2\unicode[STIX]{x1D714}_{p}}}+{\displaystyle \frac{\text{i}\left({\displaystyle \frac{g^{2}k_{s}^{^{\prime }}\tilde{k}_{p-s}}{2\unicode[STIX]{x1D714}_{p}\unicode[STIX]{x1D714}_{s}\unicode[STIX]{x1D714}_{p-s}}}+{\displaystyle \frac{g^{2}\tilde{k}_{s}k_{p-s}^{^{\prime }}}{2\unicode[STIX]{x1D714}_{p}\unicode[STIX]{x1D714}_{s}\unicode[STIX]{x1D714}_{p-s}}}\right)}{\tilde{k}_{s,u}^{x}+\tilde{k}_{p-s,l-u}^{x}-k_{p,l}^{x}}}. & & \displaystyle\end{eqnarray}$$

The imaginary components of the wavenumber are absorbed into the formulation of wave amplitudes at this point, yielding a new deterministic evolution equation which is identical to (2.1), with the exception of nonlinear interaction term ( $\tilde{W}$ ). Throughout this paper a $W$ -term notation is used as before. The two different formulations (original $W$ , and new $\tilde{W}$ ) are interchangeable in the equations.

The inclusion of spatial inhomogeneities in the mild slope equation represents another main achievement of this work. The consequence of the above is that the nonlinear interaction term becomes complex resulting in a shift of phase. The real solution for each wave amplitude is obtained by summing the evolution equation for $a_{p}$ and $a_{-p}$ . Nonlinear dissipation was previously implemented in the bispectrum evolution equation in Janssen (Reference Janssen2006) (specifically § 4.4), based on the model presented in Janssen, Herbers & Battjes (Reference Janssen, Herbers and Battjes2008). Here, the nonlinear dissipation is included into the deterministic model, which then also gets into the stochastic model. As a result nonlinear transfer to higher frequencies is significantly reduced for both deterministic and stochastic models, yielding a better agreement with the measurements (see § 4), although the stochastic model still over-predicted the energy transfer. Further discussion of the implementation of this improvement to the deterministic and stochastic models is given in appendix C.

3.2 Improved nonlinear coupling coefficients for strong spectral changes

One of the main assumptions used in Agnon & Sheremet (Reference Agnon and Sheremet1997) is that the spectrum components and the nonlinear interaction term vary slowly in space and can be taken outside of the integral in (2.6), without significant loss of accuracy. In Agnon & Sheremet (Reference Agnon and Sheremet2000) it was found that the interaction term changes rapidly as waves enter shallow water, and it cannot be simply assumed to be constant. In the work of Stiassnie & Drimer (Reference Stiassnie and Drimer2006), the model of Agnon & Sheremet (Reference Agnon and Sheremet1997) was localized by taking into account the spatial changes of the interaction term. This approach was further improved in Toledo & Agnon (Reference Toledo and Agnon2012) and Vrecica & Toledo (Reference Vrecica and Toledo2016). All of the above models assume that the spectrum components can be taken out of the integral. This assumption holds in intermediate waters. Nevertheless, in situations where wave breaking or significant nonlinear energy transfer occurs, it is not appropriate to use this assumption.

The new formulation proposed here also takes into account the change of the spectrum components over the wave path (2.6). First, equation (2.6) is rewritten as

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle Q_{s,p-s,r,l-r}=\text{e}^{-\int _{0}^{x}\unicode[STIX]{x1D708}(x^{\prime })\,\text{d}x^{\prime }}\int _{0}^{x}E_{p-s,l-u}E_{s,r}W_{p,s,p-s}\text{e}^{\int _{0}^{x^{\prime }}\unicode[STIX]{x1D708}(x^{\prime \prime })\,\text{d}x^{\prime \prime }}\,\text{d}x^{\prime }, & \displaystyle\end{eqnarray}$$

(3.16) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D708}(x)=\text{i}(k_{s,r}^{x}+k_{p-s,l-r}^{x}-k_{p,l}^{x})+\left({\displaystyle \frac{C_{g,p}^{\prime }}{2C_{g,p}}}+{\displaystyle \frac{C_{g,s}^{\prime }}{2C_{g,s}}}+{\displaystyle \frac{C_{g,p-s}^{\prime }}{2C_{g,p-s}}}\right)-(D_{p}+D_{s}+D_{p-s}).\qquad & \displaystyle\end{eqnarray}$$

Assuming that the wave field is stationary in time, the spectrum evolution equation (2.2) becomes

(3.17) $$\begin{eqnarray}\displaystyle & & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}(C_{g,p}E_{p,l})-2D_{p,l}E_{p,l}\nonumber\\ \displaystyle & & \displaystyle \quad =\mathop{\sum }_{s=p-N}^{N}\mathop{\sum }_{r=max\{l-M,-M\}}^{min\{l+M,M\}}\text{Re}[Q_{s,p-s,r,l-r}+Q_{p,s,l,t-r}+Q_{p,p-s,l,l-r} ]\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,W_{s,p-s,r,l-r}C_{g,p}=S_{nl,3,p}.\end{eqnarray}$$

Using the integrating factor method, equation (3.17) can be solved to yield

(3.18) $$\begin{eqnarray}\displaystyle E_{p,l}(x)=\text{e}^{\int _{0}^{x}D_{p,l}\,\text{d}x^{\prime }}{\displaystyle \frac{C_{g,p}(0)}{C_{g,p}(x)}}\left(E_{p,l}(0)+{\displaystyle \frac{1}{C_{g,p}(0)}}\int _{0}^{x}\text{e}^{-\int _{0}^{x^{\prime }}D_{p,l}\,\text{d}x^{\prime }}S_{nl,3,p}\,\text{d}x^{\prime }\right), & & \displaystyle\end{eqnarray}$$

where $x$ is the point at which the nonlinear coupling term $Q$ is evaluated. The objective of this transformation is to define the localized coupling term $Q$ only in terms of the local components.

The energy spectra components inside of the integral in (2.7) are replaced either by an approximated first-order linear solution (accounting only for linear shoaling), or by an exact second-order nonlinear solution. If the approximate solution is used, the localization process is the same as in Vrecica & Toledo (Reference Vrecica and Toledo2016). For the exact solution it is not possible to precalculate the values of coupling coefficients. The nonlinear corrections to the coupling term depend on the spectrum itself, and if accounted for a form of non-local approach must be used. The coupling term represents the resonant part of the bispectrum. It is possible to compute a correction at each step using an iterative approach (until the solution converges), however, a different approach is followed. It is chosen, in this paper, to transport values of the coupling coefficients.

To overcome the Nyquist limitation and save computational resources the following procedure is employed. The change of the coupling term ( $Q$ -term) is precalculated between grid points using the localization procedure. In the actual run the coupling term is transported between grid points, and its evolution is obtained by multiplying the $Q$ -term with appropriate spectral components. This approach still enables much larger grid spacing than for two-equation models, avoiding the Nyquist limitation. This advantage is negated when waves propagate into the breaking zone, as the spatial step needs to be small enough to accurately grasp wave breaking. Still, if the coupling terms are partially precalculated (see appendix C), the numerical cost of running the model in the breaking zone is essentially equal to the cost of solving the bispectrum evolution equation.

In cases of weak nonlinearity, it is possible to discard the nonlinear and dissipation terms. It is recommended to use them only in the cases of strong nonlinearity and strong breaking events, as they have little effect in mild sea conditions. The formulation of the equation (3.18) which accounts only for linear shoaling effects is defined as

(3.19) $$\begin{eqnarray}\displaystyle E_{p,l}(x)={\displaystyle \frac{C_{g,p}(0)}{C_{g,p}(x)}}E_{p,l}(0). & & \displaystyle\end{eqnarray}$$

Equation (3.19) gives a relation describing the energy function ( $E_{p,l}$ ) in terms of a initial wave energy $E_{p,l}(0)$ . Hence, by its substitution to the integral of (2.7), the integrand will contain only local wave energy components, which can be extracted from the integral to yield

(3.20) $$\begin{eqnarray}\displaystyle Q_{s,p-s,r,l-r}(x) & = & \displaystyle E_{p-s,l-r}(x)E_{s,r}(x)\sqrt{{\displaystyle \frac{C_{g,s}C_{g,p-s}}{C_{g,p}}}}\text{e}^{-\text{i}\int _{0}^{x}\unicode[STIX]{x1D709}(x^{\prime })\,\text{d}x^{\prime }}\nonumber\\ \displaystyle & & \displaystyle \times \int _{0}^{x}W_{p,s,p-s}\sqrt{{\displaystyle \frac{C_{g,p}^{x^{\prime }}}{C_{g,s}^{x^{\prime }}C_{g,p-s}^{x^{\prime }}}}}\text{e}^{\text{i}\int _{0}^{x^{\prime }}\unicode[STIX]{x1D709}(x^{\prime \prime })\,\text{d}x^{\prime \prime }}\,\text{d}x^{\prime },\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D709}(x)=k_{s,r}^{x}+k_{p-s,l-r}^{x}-k_{p,l}^{x}. & & \displaystyle\end{eqnarray}$$

The new formulation of the nonlinear coupling term, derived for the CDB model, is the second main result of this work. Even a linear shoaling approximation of the spectrum components is found to yield a significant improvement to the model derived in Vrecica & Toledo (Reference Vrecica and Toledo2016), while retaining the same numerical cost. The formulation of the model which incorporates the breaking and nonlinear terms, with all the necessary details for its implementation, is given in appendix C.

3.3 Summary of the localized triad source term formulations

In this subsection the localized versions of the nonlinear triad source term, fitting for wave forecasting models, are outlined. One is appropriate for intermediate waters. The other one is suitable for shallower waters, where wave triads approach resonance. The main model is (2.2), which is defined as

(3.22) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}E_{p,l}}{\unicode[STIX]{x2202}t_{1}}}+{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}}(C_{g,p}E_{p,l}) & = & \displaystyle 4\mathop{\sum }_{s=p-N}^{N}\mathop{\sum }_{r=max\{l-M,-M\}}^{min\{l+M,M\}}\nonumber\\ \displaystyle & & \displaystyle \times \,\text{Re}[P_{s,p-s,r,l-r}E_{s,r}E_{p-s,l-r}+P_{p,s,l,r}E_{p,l}E_{s,r}\nonumber\\ \displaystyle & & \displaystyle +\,P_{p,p-s,l,l-r}E_{p,l}E_{p-s,l-r} ]W_{s,p-s,r,l-r}C_{g,p}=S_{nl,3},\end{eqnarray}$$

where by following the equation (3.20) the coupling term is expressed as

(3.23) $$\begin{eqnarray}\displaystyle Q_{s,p-s,r,l-r}(x)=P_{s,p-s,r,l-r}(x)E_{s,r}(x)E_{p-s,l-r}(x)=P_{s,p-s,r,l-r}E_{s,r}E_{p-s,l-r}. & & \displaystyle\end{eqnarray}$$

Formulation of the model is slightly changed in this subsection in order to formulate the localized coupling coefficients solely as a function of the local depth and slope. In order to do so the bathymetry can be assumed to be of a uniform slope. For a slowly varying bathymetry this is found to be a plausible approximation. A more accurate formulation fitting for non-uniform bathymetries (for example containing underwater bars), can be found in appendix C. A gate term defined in appendix B governs which of the two regimes is used, by giving a depth ( $h_{g}$ ) at which the shallow water localization activates. Note that as the gate term never activates for flat bottom, it is allowed to divide by $h^{\prime }$ .

This formulation corresponds to the model derived in Vrecica & Toledo (Reference Vrecica and Toledo2016), where the first-order approximation (3.19) from subsection 3b is used to describe the change of the energy spectra inside of the integral. The effects of the bottom slope, derived in subsection 3a, can also be included. The localized coupling term for intermediate water is derived by applying the same intermediate water localization used in Vrecica & Toledo (Reference Vrecica and Toledo2016) to (3.20) yielding

(3.24) $$\begin{eqnarray}\displaystyle P_{s,p-s,u,l-u}(x) & = & \displaystyle \left[{\displaystyle \frac{-\text{i}W_{-s,-(p-s),u,l-u}}{\unicode[STIX]{x1D709}+\text{i}\tilde{J}_{s,u,p-s,l-u}}}+{\displaystyle \frac{W_{-s,-(p-s),u,l-u}^{^{\prime }}}{(\unicode[STIX]{x1D709}+\text{i}\tilde{J}_{s,u,p-s,l-u})^{2}}}\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\,{\displaystyle \frac{W_{-s,-(p-s),u,l-u}}{(\unicode[STIX]{x1D709}+\text{i}\tilde{J}_{s,u,p-s,l-u})^{3}}}{\displaystyle \frac{\text{d}}{\text{d}x}}(\unicode[STIX]{x1D709}+\text{i}\tilde{J}_{s,u,p-s,l-u})\right]^{h(x)}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x0394}P\text{i}_{s,p-s,u,l-u}(h(x),h^{\prime }(x)).\end{eqnarray}$$

The initial step of the integral corresponds to deep water, where the $P$ -term is assumed to be zero. The above expression is evaluated solely based on the local depth and local bottom slope.

For shallow water the value of the coefficient is obtained by applying the shallow water localization process of Vrecica & Toledo (Reference Vrecica and Toledo2016) to (3.20) yielding

(3.25) $$\begin{eqnarray}\displaystyle P_{s,p-s,u,l-u}(x) & = & \displaystyle \unicode[STIX]{x0394}P\text{i}_{s,p-s,u,l-u}(h_{g}(x),h^{\prime }(x))\text{e}^{-\text{i}\int _{h_{g}(x)}^{h(x)}\unicode[STIX]{x1D709}(x^{\prime })/h^{\prime }(x)\,\text{d}x^{\prime }}\nonumber\\ \displaystyle & & \displaystyle \times \,\sqrt{{\displaystyle \frac{C_{g,p}^{h_{g}(x)}C_{g,s}C_{g,p-s}}{C_{g,p}C_{g,s}^{h_{g}(x)}C_{g,p-s}^{h_{g}(x)}}}}+\unicode[STIX]{x0394}Ps_{s,p-s,u,l-u}(h_{g}(x),h(x),h^{\prime }(x)).\end{eqnarray}$$

In order to use local quantities, the bottom slope in the localization procedure is assumed to be linear. It is possible to adjust the model to account for other formulations of the slope, such as the power formulation ( $h=\sum _{i=0}^{N}A_{i}x^{i}$ , where $A_{i}$ are constants).

Definitions of the $\unicode[STIX]{x0394}Ps$ -term, $W$ -term, $h_{g}$ -term and $\tilde{J}$ -term are given in appendix B. To further simplify the model the exponential term containing wavenumber mismatch $(\text{e}^{-\text{i}\int _{h_{g}(x)}^{h(x)}\unicode[STIX]{x1D709}(x^{\prime })/h^{\prime }(x)\,\text{d}x^{\prime }})$ can be dropped, as the exact phasing of the $P$ -term at the gate is unknown, and its contribution is usually negligible.

4.1 Shoaling of a unidirectional spectrum over a slope

The results of the numerical simulations are compared to the field measurement of Freilich & Guza (Reference Freilich and Guza1984). A large number of point measurements was performed for laterally uniform beach. The value of the spectrum for a specific depth was obtained by averaging the measurements performed at this depth. The bathymetry from the field measurement was digitized, and is defined as

(4.1) $$\begin{eqnarray}\displaystyle h(x)=\left\{\begin{array}{@{}cc@{}}14-0.0151x & 0<x<265\\ 10-0.0238(x-265) & 65<x<370\\ 7.5-0.0266(x-370) & 370<x<490\\ 4.3-0.0162(x-490) & 490<x<570\\ 3 & 570<x\end{array}\right\}, & & \displaystyle\end{eqnarray}$$

with values expressed in metres. The incident spectrum has a variance of about 0.05  $\text{m}^{2}$ and is discretized with a bin width of 0.0025 Hz for a total of 80 bins. The initial values of the spectrum are obtained by digitizing the first station’s measurements.

This case was simulated also in Vrecica & Toledo (Reference Vrecica and Toledo2016). The simulation is repeated in this work, but with finer discretization in frequency and space domains. By accounting for spectral changes in the bispectra, the accuracy of the model is significantly improved. The simulation is performed using the linear model, ensemble averaged deterministic model of Bredmose et al. (Reference Bredmose, Agnon, Madsen and Schaffer2005), fully local formulation using (3.24) and (3.25), and solution for the specific bathymetry using (C 3) and (C 6), without accounting for the dissipation and nonlinear effects. Both of the new formulations give similar results, although the partially local formulation performs better overall. It is possible for both of the presented formulations to save values of interaction coefficients at grid points, meaning that the computational costs of running the model are identical. Pre-calculation of the coefficients is more computationally intensive for the partially local formulation though. Results of the linear model are shown in order to demonstrate the magnitude of the nonlinear effects. The performance of LTA and SPB models is also tested, with and without corrections to the tuning parameter $\unicode[STIX]{x1D6FC}$ discussed in Salmon et al. (Reference Salmon, Smit, Janssen and Holthuijsen2016). A formulation which best fitted the measurements was used in each case.

The evolution of the spectrum for the models derived in this and a previous paper (see Vrecica & Toledo Reference Vrecica and Toledo2016), is shown in figure 1. The agreement between the numerical results and the measurements is very good, accounting for the generation of secondary and tertiary peaks, as well as the generation of the infragravity wave field. The fully local model grasps the essential properties of the wave field evolution, but is not as accurate as the partially local one. The SPB model is unable to account for transfer of energy to higher and infragravity frequencies. The original formulation of the LTA model overestimates the energy transfer by a large margin, while the corrected version under-predicts it. Both of these do not account for the generation of the infragravity wave field and the third peak.

Figure 1.

Evolution of the unidirectional spectra over a sloping beach. The initial peak is at 0.06 Hz, and the wave variance $(\sum 1/2F^{2})$ is around $0.05~\text{m}^{2}$ . The spectrum is discretized in 80 bins. The grey dots represents the measurements of Freilich & Guza (Reference Freilich and Guza1984). In panels (a) and (b), the dotted red line represents the linear model, the black solid line is the ensemble averaged deterministic model, the dot-dashed brown line represents the fully local model based on (3.24) and (3.25) and the blue dashed line represents the partially local model based on (C 3) and (C 6). In panels (c) and (d) the dotted red line represents the linear model, the black solid line is the LTA model, the dot-dashed brown line is the LTA model with the correction of Salmon et al. (Reference Salmon, Smit, Janssen and Holthuijsen2016), and the blue dashed line represents the SPB model. The results are shown for depths of 5.25 m (a and c) and 4.1 m (b and d).

One thing to keep in mind is that the original formulation of the SPB model is not physical for certain conditions. This is due to a coefficient used by the model obtained in wave tank experiments $K=0.95k_{peak}-0.75$ , where $k_{peak}$ is peak wavenumber in deep water. It is possible to have a reversed transfer of energy for $k_{peak}>0.75/0.95$ , as in the following example. This issue is mitigated in most cases by using a formulation of the SPB model presented in Salmon et al. (Reference Salmon, Smit, Janssen and Holthuijsen2016), which uses local value of peak wavenumber.

4.2 Breaking of a unidirectional TMA spectrum over a slope

In this section shoaling and breaking of a unidirectional Texel, MARSEL, ARSLOE (TMA) (see Bouws et al. Reference Bouws, Günther, Rosenthal and Vincent1985) spectraum (case 4, table 1) over monotone slope is numerically investigated. The simulations correspond to the experimental case 28 from Smith (Reference Smith2004). The bathymetry is defined as

(4.2) $$\begin{eqnarray}\displaystyle h(x)=\{\begin{array}{@{}cc@{}}0.61-x/30 & 0<x\end{array}\}, & & \displaystyle\end{eqnarray}$$

with values expressed in metres. The initial significant wave height is 0.09 m, amplification factor is 100, and the peak period is 1.75 s. The input spectra and results of the measurements are obtained by digitization of the graphs from the original work.

Table 1.

Cases under consideration for evolution of the 2-D spectra over a monotone slope. The column labelled ‘Case’ defines the number of the simulation. The spectrum for case 3 is defined as the sum of the two spectra 3a and 3b evolving together with their mutual nonlinear interactions. Incidence angle defines angle of peak wave direction to the onshore direction. Peak frequency represents $T_{p}$ . $H_{s}$ is the significant wave height. The column labelled ‘Spreading’ describes how much waves propagate around the mean direction; spreading of $30^{\circ }$ means they are propagating $\pm 15^{\circ }$ around $D_{p}$ following $\cos ^{2}$ distribution.

One of the forms of the dissipation term, fitting for the case under consideration is defined as

(4.3a,b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}a_{p}}{\unicode[STIX]{x2202}x}}=D_{p}a_{p},\quad D_{p}=c_{1}{\displaystyle \frac{D_{br}}{\sum C_{g,p}a_{p}/2}}a_{p}+c_{2}{\displaystyle \frac{D_{br}}{\sum C_{g,p}a_{p}/2}}{\displaystyle \frac{f_{p}^{2}}{f_{m02}}}a_{p}, & \displaystyle\end{eqnarray}$$

(4.4a,b ) $$\begin{eqnarray}\displaystyle D_{br}={\displaystyle \frac{1}{4}}f_{m01}(0.8h)^{2}Q,\quad {\displaystyle \frac{1-Q}{\ln (Q)}}=-{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{H_{s}}{0.8h}}\right)^{2}, & & \displaystyle\end{eqnarray}$$

where $f_{m}$ is the mean weighted frequency, $f_{m2}$ is the mean square weighted frequency, $c_{1}$ and $c_{2}$ are tunable coefficients and $H_{s}$ is the significant wave height. First a total loss of energy $(D_{br})$ is calculated based on the formulation of Battjes & Janssen (Reference Battjes and Janssen1978). The most common approach is to equally distribute loss of energy across all spectral components, such as in Eldeberky & Battjes (Reference Eldeberky and Battjes1996) where $c_{1}=1$ and $c_{2}=0$ . However, the model shows better agreement with the frequency dependent formulation of Mase & Kirby (Reference Mase and Kirby1992) where $c_{1}=0.5$ and $c_{2}=0.5$ . It should be taken into account that the energy dissipation term $D_{p}$ is based on heuristics, therefore some tuning is necessary when changing the base model.

Initially, a secondary peak forms at double the frequency of the initial peak. As the water gets shallower, the waves start to break, leading to a significant loss of energy. The formation of a third peak at triple the frequency of the initial peak is not as clearly pronounced. A very clear peak forms in the infragravity range, with magnitude comparable to that of the secondary peak. Evolution of the initial spectra is modelled using an ensemble average obtained using the deterministic model of Bredmose et al. (Reference Bredmose, Agnon, Madsen and Schaffer2005), an ensemble average of the deterministic model presented in § 3.1, and for the results obtained using partially local formulation of the stochastic model, with the corrections presented in § 3.1.

As the closure of Benney & Saffman (Reference Benney and Saffman1966) over-predicts energy transfer, a stochastic simulation is also performed using a closure of Holloway (Reference Holloway1980). The closure is based on the empirical observations, and adds additional dissipation to the bispectrum evolution equation. Similar formulation to the one presented in Janssen (Reference Janssen2006) § 4.4 is used, however, only additional linear dissipation was added to the bispectrum evolution equation. This is straightforward to implement, only the dissipation coefficients in the bispectrum evolution equation (2.4) are modified (see appendix C). In total there are three different dissipation effects that are considered. The same as in Janssen, Herbers & Battjes (Reference Janssen, Herbers and Battjes2006), we include a dissipation term ( $D$ ) in the deterministic model, and additionally ( $D2$ ) in the bispectrum evolution equation (to model the effects of decorrelation). The difference is that we also consider the effect of dissipation in the nonlinear interaction term ( $W$ ) through implementation as an imaginary part of the wave number (see § 3.1).

Results are shown in figure 2. In figure 3 results obtained using LTA and corrected SPB formulation are shown. Evolution of wave skewness and asymmetry, for an ensemble average of deterministic model, is shown in figure 2. It should be pointed out that although the agreement of the new models with the measurements is good, further work is needed in understanding wave dissipation and statistics in the nearshore.

Figure 2.

Evolution of the unidirectional TMA spectra over the monotone slope, corresponding to the experimental case 28 from Smith (Reference Smith2004) (grey dots). The primary peak is at 0.58 Hz and the initial significant wave height is 9 cm. The ensemble average of deterministic model from § 3.1 is plotted as a solid blue line for the newly derived model, and dashed red colour for the ensemble average of Bredmose et al. (Reference Bredmose, Agnon, Madsen and Schaffer2005); the dash-dash-dot grey line is the linear model, the dot-dash-dot brown line is the stochastic model and the dotted black line is the output of the stochastic model with the closure of Holloway. The spectrum at water depths of 18 cm (a), 15 cm (b), 12 cm (c) and 9 cm (d) is shown. Evolution of ensemble averaged wave skewness (e) and asymmetry (f) is also shown, for both versions of the deterministic model.

Figure 3.

Comparison between existing stochastic models for the same case shown in figure 2. In each panel a solid red line represents the LTA model, a dash-dot-dash grey line is the linear model, a solid blue line is the new formulation of deterministic model and a dashed blue line is the corrected SPB model. The grey dots represent the measurements. The spectrum at water depths of 18 cm (a), 15 cm (b), 12 cm (c), and 9 cm (d) is shown.

The ensemble average of the deterministic mild slope equation derived in § 3.1 shows much better agreement with the measurements in all frequency ranges of the spectra. The original formulation of Bredmose et al. (Reference Bredmose, Agnon, Madsen and Schaffer2005) significantly over-predicts the energy transfers. The stochastic models of the authors grasp the overall evolution of the wave field up to the 15 cm depth. However, as the water gets shallower the output of the stochastic model deviates significantly from measurements. The statistical properties of the deterministic model’s ensemble average shows that the closure of Benney & Saffman (Reference Benney and Saffman1966) fails in the breaking zone. A closure of Holloway (Reference Holloway1980) is applied to the model, as discussed in Janssen (Reference Janssen2006) § 4.4. The closure improved the prediction of the stochastic model, with regard to the primary and secondary peak generation. An initial examination shows that the error made by the closure of Benney & Saffman (Reference Benney and Saffman1966) is more evenly distributed, this is, however, beyond the scope of this work.

The LTA model and the corrected SPB model account for the generation of the secondary peak, but they fail to predict the generation of the infragravity wave peak. They also predict a very narrow spectral shape, as opposed to the measurements. The present CDB model is capable of grasping both.

Evolution of ensemble averaged wave skewness and asymmetry is shown in figure 2(e,f). Both skewness and asymmetry oscillate in deeper water, however, as the water gets shallower the two models diverge. For the model accounting for the slopes of the amplitudes, the skewness tends to asymptotically go to zero, while the asymmetry continues to increase, which is in line with field observations (see Elgar & Guza Reference Elgar and Guza1985). This is not the case with the other formulation of the deterministic model, in which the opposite happens.

4.3 Evolution of directionally spread JONSWAP spectra over a monotonous slope

The previous case shows good agreement of the model with measurements for 1-D approximation. The aim of this subsection is to verify qualitative behaviour of the model, and to investigate essential properties of 2-D wave shoaling. With this aim, four different cases are simulated, and all of them are presented in table 1. For all cases the spectrum is defined using the original formulation of the Joint North Sea Wave Observation Project (JONSWAP) spectra from Hasselmann, Dunckel & Ewing (Reference Hasselmann, Dunckel and Ewing1973) with $\cos ^{2}$ directional distribution, and an amplification factor of 3.3.

The bathymetry is defined as

(4.5) $$\begin{eqnarray}\displaystyle h(x)=\left\{\begin{array}{@{}cc@{}}50-0.05x & 0<x<900\\ 5 & x>900\end{array}\right\}, & & \displaystyle\end{eqnarray}$$

with values expressed in metres. As the goal of this subsection is to gain a better understanding of directional properties of the wave field, a simple local stochastic model (3.24) and (3.25) is used in the simulations.

It is important to note that discretization needs to be done carefully, especially in 2-D models. In order to account for resonant and near resonant interactions, the discretization needs to be fine enough (see Lvov, Nazarenko & Pokorni Reference Lvov, Nazarenko and Pokorni2006), which is fulfilled in this work. An issue that arises as waves enter shallow water is that no discretization is fine enough to account for resonant interactions. However, it was shown in Annenkov & Shrira (Reference Annenkov and Shrira2006) that accounting for near resonant interactions is of even higher importance.

In case 1, a single 2-D JONSWAP spectrum with spreading of $\pm 15^{\circ }$ propagates over a monotone slope. The initial significant wave height is $H_{s}=3.85~\text{m}$ . The evolution of the single system propagating towards shore is shown in figure 4. As expected a secondary peak forms at double the frequency of the initial peak. A small peak at triple the frequency forms near the end of the propagation. It can be seen that the 2-D nonlinear spectral broadening of the peak discussed in Toledo (Reference Toledo2013), can be accounted for using the derived 2-D stochastic model. Similar results are given by Ma et al. (Reference Ma, Chen, Ma and Dong2017) using a Boussinesq model.

Figure 4.

Amplitude spectrum of shoaling 2-D wave field. The system propagates onshore, has a peak at 0.095 Hz and directional spreading of $\pm 15^{\circ }$ . The initial significant height is 3.85 m. $k_{y}$ is the lateral wavenumber. White lines represent constant angle of incidence to the shoreline. The results are shown for water depths of 50 m (a), 20 m (b), 14 m (c) and 8 m (d).

In case 2, a single 2-D JONSWAP spectrum with an initial angle of attack of $35^{\circ }$ and spreading of $\pm 10^{\circ }$ propagates over a monotone slope. The initial significant wave height is $H_{s}=3.75~\text{m}$ . The results are shown in figure 5. It should be noted that the intensity of the nonlinear interactions is reduced by a small amount for the system propagating at an angle to the onshore direction. This can be explained by the smaller effective bottom slope, and by the inability to close the direct resonance conditions as in the 1-D case.

Figure 5.

Amplitude spectrum of shoaling 2-D wave field. The system propagates at $35^{\circ }$ angle to the shore, has a peak at 0.1 Hz and directional spreading of $\pm 10^{\circ }$ . The initial significant height is 3.75 m. $k_{y}$ is the lateral wavenumber. White lines represent constant angles of incidence to the shoreline. The results are shown for water depths of 50 m (a), 20 m (b), 14 m (c) and 8 m (d).

In case 3, a combination of the two previous JONSWAP spectra propagate over the same monotone slope. The two systems interact with themselves and with each other to generate three new secondary peaks with different intensities and directionalities. The spectrum also significantly broadens due to the subharmonic interactions. In order to show the effects of the interaction between the two systems, the difference between the two systems’ combined evolution and a sum of the two systems evolving independently is shown in figures 6 and 7.

Figure 6.

Amplitude spectrum of shoaling 2-D wave field, and a realization of the deterministic model. This figure shows the evolution of the spectrum for the case of directionally spread seas. The initial spectrum is the sum of initial spectra from figures 4 and 5. White lines represent constant angle of incidence to the shoreline. The results are shown for water depths of 50 m (a), 20 m (b), 14 m (c) and 8 m (d). A realization of a deterministic model for a two system case at depth of 8 m is also shown (e).

Figure 7.

Difference of the amplitude spectrum for the evolution of the whole system and sum of evolutions of two independent systems. The first system propagates at $35^{\circ }$ angle to the shore, has a peak at 0.1 Hz and directional spreading of $\pm 10^{\circ }$ . The second system propagates onshore, has a peak at 0.095 Hz and directional spreading of $\pm 30^{\circ }$ . The initial significant heights are 3.85 m and 3.75 m, respectively. White lines represent constant angles of incidence to the shoreline. The results are shown for water depths of 50 m (a), 20 m (b), 14 m (c) and 8 m (d).

Figure 8.

Wave triad interaction schematic in frequency–wavenumber space showing broadening of the spectrum for directionally spread seas. Generation of wave components in the opposite direction is shown using a solid blue line. Interactions between the two systems are shown using a dot-dashed blue line. Interactions of the system with itself are shown using a dashed red line.

Interesting directional properties of the wave field can be deduced from these figures. The directionality of the secondary peaks depends on the sum of the lateral wavenumbers of the primary peaks. More importantly, as a result of the subharmonic interactions between the primary peaks with the secondary peaks, the spectrum significantly broadens. Although this broadening is second-order in energy, it is first-order in direction. A similar behaviour, whose cause was not explained, can be found in the field measurements of Herbers et al. (Reference Herbers, Orzech, Elgar and Guza2003). A known cause of nonlinear spectrum broadening is the scattering due to bottom ripples (see Fabrice & Herbers Reference Fabrice and Herbers2002). Here, it is shown that also a uniform slope can cause nonlinear broadening. This property can also be qualitatively observed in an output of a single deterministic run corresponding to case 3 (see figure 6 e).

New wave components are generated at the frequency of the primary peaks, which have large attack angles compared to the original ones. Perhaps the most interesting effect is the generation of a peak with the mirrored angle of propagation. This is the result of a subharmonic interaction of a primary peak propagating at an angle and a secondary peak propagating directly to the shoreline. A schematic description of this mechanism is given in figure 8. Wave component $a_{1,0}$ generates its superharmonic $a_{2,0}$ through self-interaction. This allows the $a_{2,0}$ component and $a_{1,1}$ component to generate through subharmonic interaction $a_{1,-1}$ , which has an opposite angle of propagation compared to harmonic $a_{1,1}$ . This effect is shown using a solid blue line. Other interactions which only happen as a result of interactions between two wave systems are shown using dashed blue line. Interactions which do not depend on the existence of the other system are shown using a dashed red line. As discussed in Salmon et al. (Reference Salmon, Smit, Janssen and Holthuijsen2016), one of the biggest limitations of SPB and LTA is that they only account for interactions of co-propagating waves. This shows the importance of the capability of the current model for accounting for nonlinear wave interactions of waves propagating in all onshore directions.