Triads

This section outlines the computation of nonlinear interactions at shallow water by means of two methods.


The first method is based on a quadratic wave theory that express the nonlinear energy transfer using the interaction coefficient. With this method, three different formulations will be discussed: the Full Triad Interaction Model (FTIM), the Stochastic Parametric model based on Boussinesq equations (SPB) and the Lumped Triad Approximation (LTA).


The second method is called the Discrete Collinear Triad Approximation (DCTA) method and heuristically captures the $k^{-4/3}$ spectral tail at shallow water depths, while generating all transient sub and super harmonics.


Quadratic formulations


The starting point of the quadratic wave theory is to describe the free surface elevation of random spread but long-crested waves propagating over a mildly sloping bed $d(x)$ using the following one-dimensional expression (see, e.g. Freilich and Guza, 1984; Eldeberky, 1996; Herbers and Burton, 1997; Akrish et al, 2024)

  $$\eta(x,t) = \sum_{p=-\infty}^{\infty} A_p(x)\,\exp{[{\rm i} (\omega_p t - \psi_p(x))]}$$ (2.91)

where $A_p$ is the complex Fourier amplitude of the $p-$th harmonic, $\omega_p$ is the $p-$th angular frequency and $\psi_p(x)$ is the linear phase of $p-$th wave component with $d\psi_p/dx = k_p$ the wave number. The latter two numbers satisfy the linear dispersion relation: $\omega^2_p = g k_p \tanh{(k_p\,d)}$. Note that $A_{-p} = A^{\star}_p$ with $\star$ indicating the complex conjugate and also $\omega_{-p} = -\omega_p$ and $\psi_{-p} = -\psi_p$.


The amplitude $A_p$ is a slow function of $x$ due to shoaling and nonlinear interactions. Its evolution equation is given by (Madsen and Sørensen, 1993; Eldeberky ,1996)

  $$\frac{dA_p}{dx} = -S_p A_p - {\rm i} \sum_{m=-\infty}^{\infty} R_{(m,p-m)}\,A_m\,A_{p-m}\,\exp{[-{\rm i}(\psi_m+\psi_{p-m}-\psi_p)]}$$ (2.92)

with $S_p$ the shoaling coefficient and $R_{(m,p-m)}$ the interaction coefficient. Note that $R_{(k,l)}$ is real and symmetric: $R_{(l,k)} = R_{(k,l)}$ with the subscripts $l$ and $k$ denoting dummy indices. Various expressions for the interaction coefficient will be discussed later. The first and second term on the right hand side represent linear shoaling and nonlinear triad interactions, respectively. In the following we focus on the interaction term and therefore omit the shoaling part.


Let us introduce the following complex Fourier amplitude

  $$C_p (x) = A_p (x) \exp{[-{\rm i}\,\psi_p(x)]}$$ (2.93)

then one can derive the following equation

  $$\frac{dC_p}{dx} = -{\rm i} \, k_p C_p - {\rm i} \sum_{m=-\infty}^{\infty} R_{(m,p-m)}\,C_m\,C_{p-m}$$ (2.94)

where the first term on the right hand side describes the linear phase evolution and the second term the quadratic interactions. Usually, this second term can be splitted into two terms that describe the energy transfer to component $p$ from sum interactions of two components $m$ and $p-m$ and difference interactions between $p+m$ and $m$, respectively, as follows

  $$\frac{dC_p}{dx} = -{\rm i} \, k_p C_p - {\rm i} \sum_{m=1}^{p-1} ... ...\,C_{p-m} - 2\, {\rm i} \sum_{m=1}^{\infty} R_{(p+m,-m)}\,C^{\star}_m\,C_{p+m}$$ (2.95)

With the above deterministic evolution equation for the complex amplitude we are able to derive a set of stochastic equations that describe the spatial evolution of the high-order moments (Eldeberky, 1996). In this section we restrict ourselves to the second-order and third-order moments and their evolution equations. The second-order moment is defined as

  $$E_p = \langle C_p\,C^{\star}_p \rangle$$ (2.96)

and is interpreted as the discrete spectral energy. The operator $\langle\cdot\rangle$ denotes the expected value. The third-order moment is known as the discrete bispectrum (Hasselmann et al, 1963) and its definition is given by

  $$B_{m,p-m} = \langle C_m\,C_{p-m}\,C^{\star}_p \rangle$$ (2.97)

The bispectrum accounts for the phase coupling among the three components involved in a triad.


Since

  $$\frac{dE_p}{dx} = \biggl \langle C^{\star}_p\,\frac{dC_p}{dx} \biggr \rangle + \biggl \langle C_p\,\frac{dC^{\star}_p}{dx} \biggr \rangle$$ (2.98)

the reader may verify by substituting Eq. (2.95) that the evolution equation of the discrete energy spectrum is given by

  $$\frac{dE_p}{dx} = 2\,\sum_{m=1}^{p-1} R_{(m,p-m)} \, {\rm Im}(B_{m,p-m}) - 4\,\sum_{m=1}^{\infty} R_{(p+m,-m)} \, {\rm Im}(B_{m,p})$$ (2.99)

In a similar manner, an evolution for the bispectrum is obtained by taking the derivative of $B_{m,p-m}$ to $x$, as follows

  $$\frac{dB_{m,p-m}}{dx} = \biggl \langle C_{p-m}\,C^{\star}_p\, \fr... ...\rangle + \biggl \langle C_m\,C_{p-m}\,\frac{dC^{\star}_p}{dx} \biggr \rangle$$ (2.100)

Substitution of Eq. (2.94) yields (Herbers and Burton, 1997)

$$\begin{eqnarray*} \frac{dB_{m,p-m}}{dx} = &-& {\rm i}\,\Delta k\,B_{m,p-m} \non... ..._{(n,p-m-n)}\,T_{m,n,p-m-n} - R_{(n,m-n)}\,T_{n,m-n,p-m} \right] \end{eqnarray*}$$

where $\Delta k = k_{p-m}+k_m-k_p$ is the wave number mismatch and $T_{m,n,p-m-n}$ is the fourth-order moment or the discrete trispectrum and is defined as

  $$T_{m,n,p-m-n} = \langle C_m\,C_n\,C_{p-m-n}\,C^{\star}_p \rangle$$ (2.101)

Instead of deriving an evolution equation for the trispectrum that contains fifth-order moments, a closure hypothesis is invoked that close the coupled set of evolution equations of high-order moments.


Basically, a trispectrum can generally be expressed in terms of quadratic products of energy spectrum plus a so-called fourth-order cumulant (Eldeberky, 1996), as follows

  $$T_{m,n,p-m-n} = E_m E_n \delta_{np} + E_n E_p \delta_{mp} + E_m E_p \delta_{n(-m)} + c_4$$ (2.102)

with $\delta_{lk}$ the Kronecker delta. With respect to the cumulant term $c_4$, the following two hypotheses are introduced: i) the quasi-Gaussian hypothesis (or the quasi-normal closure) and the closure hypothesis of Holloway (1980).


Under the assumption of weak nonlinearity, the sea state tends to be near-Gaussian and, consequently, the fourth-order cumulant can be neglected (Eldeberky, 1996), that is, $c_4 = 0$. With Eq. (2.102), the evolution equation of $B_{m,p-m}$ reduces to

  $$\frac{dB_{m,p-m}}{dx} = -{\rm i}\,\Delta k\,B_{m,p-m} + 2\,{\rm i... ...p-m)}\,E_m\,E_{p-m} - R_{(p,-m)}\,E_m\,E_p - R_{(m-p,p)}\,E_{p-m}\,E_p \right]$$ (2.103)

The closure hypothesis of Holloway (1980) assumes that $c_4$ is proportional to $B_{m,p-m}$, implying the following evolution equation for the bispectrum

$$\begin{eqnarray*} \frac{dB_{m,p-m}}{dx} = &-& {\rm i}\,\Delta k\,B_{m,p-m}\nonu... ...{(m-p,p)}\,E_{p-m}\,E_p \right] \nonumber \\ &-& K\, B_{m,p-m} \end{eqnarray*}$$

with $K > 0$ an empirical parameter (unit: m$^{-1}$) which allows for relaxing the bispectrum towards a Gaussian state. Note that this is rather a crude approximation since it does not account for the influence of the interaction coefficients.


The approximate equation of the bispectrum can be written as

  $$\frac{dB_{m,p-m}}{dx} = -{\rm i}\,\Psi\,B_{m,p-m} + 2\,{\rm i} \... ...p-m)}\,E_m\,E_{p-m} - R_{(p,-m)}\,E_m\,E_p - R_{(m-p,p)}\,E_{p-m}\,E_p \right]$$ (2.104)

with

  $$\Psi = \Delta k \qquad \quad \textnormal{(quasi-Gaussian hypothesis)}$$ (2.105)

or

  $$\Psi = \Delta k -{\rm i}\,K\qquad \textnormal{(Holloway's hypothesis)}$$ (2.106)

Following the solution procedure of Eldeberky (1996), the steady-state solution of the bispectrum is given by

  $$B_{m,p-m} = \frac{2}{\Psi}\,\left[ R_{(m,p-m)}\,E_m\,E_{p-m} - R_{(p,-m)}\,E_m\,E_p - R_{(m-p,p)}\,E_{p-m}\,E_p \right]$$ (2.107)

It should be noted that the bispectrum depends only on the local spectrum which is reasonable given the weak nonlinearity and small bed slope (Herbers and Burton, 1997).


We first proceed with the quasi-Gaussian hypothesis. Since the bispectrum is generally a complex number we may write it as follows

  $$B_{l,k} = \vert B_{l,k}\vert \, \exp{[-{\rm i} \beta_{l,k}]}$$ (2.108)

with $\vert B_{l,k}\vert$ the bicoherence and $\beta_{l,k}$ the biphase. Consequently, the evolution equation (2.99) for the discrete spectrum is rewritten as

  $$\frac{dE_p}{dx} = 2\,\sum_{m=1}^{p-1} R_{(m,p-m)} \, \vert B_{m,p... ...\,\sum_{m=1}^{\infty} R_{(p+m,-m)} \, \vert B_{m,p}\vert\,\sin{(-\beta_{m,p})}$$ (2.109)

The bicoherence is obtained directly from Eq. (2.107) and is given by

  $$\vert B_{m,p-m}\vert = \frac{2}{\vert k_{p-m}+k_m-k_p\vert}\,\big... ...\,E_m\,E_{p-m} - R_{(p,-m)}\,E_m\,E_p - R_{(m-p,p)}\,E_{p-m}\,E_p \biggr \vert$$ (2.110)

(note that the wave number mismatch $\Delta k$ is a real number). The above approach using the quasi-Gaussian hypothesis serves as the basis for modelling the spectral source term for triad interactions. We come back to this point later. This also includes the parametrization of the biphase.


Next, we consider the closure hypothesis of Holloway (1980). Now with $\Psi = \Delta k -{\rm i}\,K$ the steady-state solution (2.107) is a complex number and hence

  $${\rm Im}(B_{m,p-m}) = \frac{2K}{(\Delta k)^2+K^2}\,\left[ R_{(m,p-m)}\,E_m\,E_{p-m} - R_{(p,-m)}\,E_m\,E_p - R_{(m-p,p)}\,E_{p-m}\,E_p \right]$$ (2.111)

Using this expression and the spectrum evolution equation (2.99), we obtain another method for the computation of the triad source term in spectral models. This will be discussed in the next section.


Spectral formulations for the triad source term $S_{\rm nl3}$


The quadratic formulations presented in the previous section are expressed in terms of the complex Fourier amplitude $C_p$. To derive the spectral source term for triads a relation between $C_p$ and the variance density spectrum $E(f_p)$ must be established. The single-sided variance density spectrum at frequency $f_p$ is given by

  $$E(f_p) = \frac{2\,\langle C_p\,C^{\star}_p \rangle}{\Delta f_p} = \frac{2\,E_p}{\Delta f_p}\,,\quad p=1,2,\cdots$$ (2.112)

with $\Delta f_p$ the $p-$th frequency step. Note that the frequency resolution may not be constant.


With the above first method based on the quasi-Gaussian hypothesis, one can derive the evolution equation of spectral energy density due to triad interactions from Eqs. (2.109) and (2.110), as follows

$$\begin{eqnarray*} \frac{dE(f_p)}{dx} &=& 2\,\sum_{m=1}^{p-1} R_{(m,p-m)} \, \fr... ...ert k_{p}+k_m-k_{p+m}\vert}\, Q(f_m,f_p) \,\sin{(-\beta_{m,p})} \end{eqnarray*}$$

with $Q(f_l,f_k)$ consisting of quadratic products of variance density at different frequencies:

  $$Q(f_m,f_{p-m}) = R_{(m,p-m)}\,E(f_m)\,E(f_{p-m}) - R_{(p,-m)}\,E(f_m)\,E(f_p) - R_{(m-p,p)}\,E(f_{p-m})\,E(f_p)$$ (2.113)

and

  $$Q(f_m,f_p) = R_{(m,p)}\,E(f_m)\,E(f_p) - R_{(p+m,-m)}\,E(f_m)\,E(f_{p+m}) - R_{(-p,p+m)}\,E(f_p)\,E(f_{p+m})$$ (2.114)

Here, two remarks needs to be made. First, the above derivation is an approximation unless the frequency steps are equidistant. However, the resulting model error may be compensated by means of a calibration as introduced below. Second, in case of $Q < 0$ subharmonic transfers take place with a mismatch of 180$^o$ in the biphase. This can be taken into account by just ignoring the absolute sign of the $Q$ function.


The final expression for the $S_{\rm nl3}$ source term that conserves the energy flux $c_g(f_p)E(f_p)$ is obtained after multiplying the above result with a calibration factor $\alpha$ and the group velocity $c_g(f_p)$, as follows

$$S_{\rm nl3}(f_p) = 2\, \alpha_{\mbox{\tiny FTIM}} \, c_{g,p} \left[\sum_{m=1}^{p-1} ... ...}{\vert k_{p-m}+k_m-k_p\vert}\, Q(f_m,f_{p-m}) \,\sin{(-\beta_{m,p-m})} \right.$$

$$\left. - 2\,\sum_{m=1}^{\infty} R_{(p+m,-m)} \, \frac{\Delta f_m}{\vert k_{p}+k_m-k_{p+m}\vert}\, Q(f_m,f_p) \,\sin{(-\beta_{m,p})} \right]$$ (2.115)

This model accounts for all collinear sum and difference interactions and is referred to as the Full Triad Interaction Model (FTIM). This method is implemented in SWAN version 41.51. Note that it is expected that the calibration factor $\alpha_{\mbox{\tiny FTIM}}$ will be ${\cal O}(1)$.


In a similar way, a model based on the closure hypothesis of Holloway (1980) can be obtained called the SPB model as proposed by Becq-Girard et al (1999). The associated net source term due to triads is given by

$$S_{\rm nl3}(f_p) = 2\, \alpha_{\mbox{\tiny SPB}} \, K\, c_{g,p} \left[\sum_{m=1}^{p-... ...{(m,p-m)} \, \frac{\Delta f_m}{\Delta k^2_{m,p-m}+K^2}\, Q(f_m,f_{p-m}) \right.$$

$$\left. - 2\,\sum_{m=1}^{\infty} R_{(p+m,-m)} \, \frac{\Delta f_m}{\Delta k^2_{m,p}+K^2}\, Q(f_m,f_p) \right]$$ (2.116)

with $\Delta k_{m,p-m} = k_{p-m}+k_m-k_p$, etc. Again, $\alpha_{\mbox{\tiny SPB}} = {\cal O}(1)$.


Based on a calibration study using a series of laboratory datasets, Becq-Girard et al (1999) proposed the following expression for tuning parameter $K$:

  $$K = 0.95\,k_{\rm op} - 0.75$$ (2.117)

where $k_{\rm op}$ is the offshore peak wave number. However, Salmon (2016) argued the difficulty of defining $k_{\rm op}$ in case of realistic applications and suggested another expression that also prevents $K<0$, namely, $K = 0.95\, k_{\rm peak}$ with $k_{\rm peak}$ the local peak wave number. This is implemented in SWAN.


LTA model


The Lumped Triad Approximation (LTA) of Eldeberky (1996), which is a slightly adapted version of the Discrete Triad Approximation (DTA) of Eldeberky and Battjes (1995), is effectively a simplification of FTIM (2.115) with the aim to reduce the computational effort. In this regard, the following four assumptions are made (Eldeberky, 1996):

1. The three different interaction coefficients in Eq. (2.113) are equal to each other: $R_{(m,p-m)}= R_{(p,-m)}= R_{(m-p,p)}$.
2. The set of sum and difference interactions is replaced by the self-self interaction only. Also, the frequency increment $\Delta f$ associated with each individual interaction is replaced by an effective interaction bandwidth $\delta f$.
3. The wave number mismatch $\Delta k$ is scaled with $k(f_p)$ and $\delta f$ is scaled with $f_p$. Hence, their ratio $\delta f / \vert\Delta k\vert$ scales with the phase speed $c_p$.
4. Only transfers from lower to higher harmonics are considered.

Based on the first assumption, Eq. (2.113) becomes

  $$Q(f_m,f_{p-m}) = R_{(m,p-m)} \Bigl [E(f_m)\,E(f_{p-m}) - E(f_m)\,E(f_p) - E(f_{p-m})\,E(f_p)\Bigr ]$$ (2.118)

Similarly,

  $$Q(f_m,f_p) = R_{(m,p)} \Bigl [ E(f_m)\,E(f_p) - E(f_m)\,E(f_{p+m}) - E(f_p)\,E(f_{p+m}) \Bigr ]$$ (2.119)

Furthermore, the first summation term of Eq. (2.115) is replaced by one single term that corresponds to the self interaction at frequency $f_{p/2}$ ($m = p/2$) while replacing $\Delta f_m$ by the effective interaction bandwidth, as follows

  $$S^+_{\rm nl3} (f_p) = R_{(p/2,p/2)} \, \frac{\delta f}{\vert 2\,k_{p/2}-k_p\vert}\,Q(f_{p/2},f_{p/2}) \,\sin{(-\beta_{p/2,p/2})}$$ (2.120)

Together with the first, third and fourth assumptions, this leads to

  $$S^+_{\rm nl3} (f_p) = c_p \, R^2_{(p/2,p/2)} \, \max \Bigl [ \, 0,\,E^2(f_{p/2}) - 2\,E(f_{p/2})\,E(f_p) \, \Bigr ] \,\sin{(-\beta_{p/2,p/2})}$$ (2.121)

Owing to the fourth assumption, the resulting term represents only a positive contribution to $E(f_p)$ due to self-self interaction at frequency $f_{p/2}$.


Similarly, the second summation term of Eq. (2.115) is simplified to ($m = p$)

  $$S^-_{\rm nl3} (f_p) = c_{2p} \, R^2_{(p,p)} \, \max \Bigl [ \, 0,\,E^2(f_p) - 2\,E(f_p)\,E(f_{2p}) \, \Bigr ] \,\sin{(-\beta_{p,p})}$$ (2.122)

that describes the interaction of $E(f_p)$ with itself leading to energy transfer to $E(f_{2p})$. Note that $S^-_{\rm nl3} (f_p) = S^+_{\rm nl3} (f_{2p})$ which is the result of the first assumption. This further enhances the computational efficiency.


The final expression for LTA reads

  $$S_{\rm nl3}(f_p) = \alpha_{\mbox{\tiny LTA}} \, c_{g,p} \left[ S^+_{\rm nl3} (f_p) - 2\, S^-_{\rm nl3} (f_p) \right]$$ (2.123)

where $\alpha_{\mbox{\tiny LTA}}$ is a tunable scaling factor that controls the strength of triad interactions.


The LTA model generates second and possibly fourth and eighth higher harmonics, which are, however, persistent over large distances away from the surf zone. Therefore, the source term $S_{\rm nl3}(f_p)$ is computed only for frequencies $f_p < 2.5 \tilde{f}$ with $\tilde{f}$ the mean frequency (see Eq. 2.67). A way to prevent this unwanted situation is to add another triad interaction, as outlined in the next section.


Extended LTA model


Since version 41.51, an extended form of the LTA is implemented. Here, another triad interaction is added that involves the following three components: $p/3$, $2p/3$ and $p$. So, with $m=2p/3$ and according to the first summation term of Eq. (2.115), the associated contribution of the sum interaction at $f_p$ reads

  $$S^{++}_{\rm nl3} (f_p) = c_p \, R^2_{(2p/3,p/3)} \, \max \Bigl [ ... ...eft ( E(f_{p/3}) + E(f_{2p/3}) \right ) \, \Bigr ] \,\sin{(-\beta_{2p/3,p/3})}$$ (2.124)

which represents the superharmonic transfer from $E(f_{p/3})$ and $E(f_{2p/3})$ to $E(f_p)$. The contribution associated with the second term of Eq. (2.115) (difference interaction) is given by ($m=2p$)

  $$S^{--}_{\rm nl3} (f_p) = c_{3p} \, R^2_{(2p,p)} \, \max \Bigl [ \... ...3p}) \, \left ( E(f_p) + E(f_{2p}) \right ) \, \Bigr ] \,\sin{(-\beta_{2p,p})}$$ (2.125)

This expression clearly displays the energy transfer from components $p$ and $2p$ to $3p$. Again, notice that $S^{--}_{\rm nl3} (f_p) = S^{++}_{\rm nl3} (f_{3p})$.


Both contributions can be added to the original LTA leading to the final expression:

  $$S_{\rm nl3}(f_p) = \alpha_{\mbox{\tiny LTA+}} \, c_{g,p} \left[ S... ..._{\rm nl3} (f_p) + S^{++}_{\rm nl3} (f_p) - 2\, S^{--}_{\rm nl3} (f_p) \right]$$ (2.126)

with $\alpha_{\mbox{\tiny LTA+}}$ a calibration factor.


Parametrization of the biphase


Both the FTIM method and the (extended) LTA method express the bispectrum in terms of the biphase $-\pi/2 \leq \beta_{l,k} \leq 0$ where $l$ and $k$ are the dummy indices. Instead of deriving an equation for the biphase (see, e.g. Reniers and Zijlema, 2022) a parametrization is proposed. According to Eldeberky (1996), the biphase depends only on the spectral Ursell number, as follows

  $$\beta = -\frac{\pi}{2} + \frac{\pi}{2} \tanh (\frac{m}{Ur})$$ (2.127)

with Ursell number $Ur$ given by

  $$Ur = \frac{g\,H_{m0}}{8\sqrt{2}} \left( \frac{T_{m01}}{\pi\,d} \right)^2$$ (2.128)

and $m$ a tunable coefficient. Note that there is no dependence on the wave frequencies, that is, $\beta_{l,k} = \beta$ for all $l$ and $k$.


Eldeberky and Battjes (1995) proposed $m = 0.2$ based on a laboratory experiment. However, our recent experience shows that this relatively low value triggers some instability that artificially amplifies higher harmonics in the triad computation. Yet it will be less prone to error if the value of $m$ is increased. As suggested by Doering and Bowen (1995), the optimal agreement of Eq. (2.127) with the data of some field measurements is obtained with a value of $m=0.63$, which also reflects a robust numerical performance.


Recently, De Wit (2022) proposed a biphase parametrization that is derived from SWASH experiments demonstrating the dependence of biphase on local bed slope and local peak wave period. This parametrization is implemented in SWAN version 41.45. Note that this parametrization allows the biphase values to be positive, which potentially makes it possible to include the effect of recurrence (i.e. to transfer wave energy back to the primary peak).


Interaction coefficient


The second term of the quadratic model (2.94) describes the nonlinear interaction between two shoaling waves with frequencies $f_m$ and $f_{p-m}$ while the third component at $f_p$ is either a superharmonic or a subharmonic due to the (near) resonance condition. The energy transfer is controlled by the interaction coefficient $R_{(m,p-m)}$ which is derived from a time-domain wave model. Various formulations exist for the quadratic model and an overview is provided by Akrish et al (2024).


In SWAN 41.51, four different formulations for $R_{(m,p-m)}$ are implemented. In the earlier SWAN versions, the interaction coefficient based on the Boussinesq-wave theory of Madsen and Sørensen (1993) was implemented for the purpose of the LTA as proposed by Eldeberky (1996). The other three formulations are due to Freilich and Guza (1984), Bredmose et al (2005) and, recently published in Akrish et al (2024), the QuadWave model. The four interaction coefficients are summarized below.


Freilich and Guza (1984) were the first to propose a quadratic model based on the classical Boussinesq-wave model of Peregrine (1967). Their interaction coefficient is given by

  $$R_{(m,p-m)} = -\frac{3}{4}\,\frac{\sigma_p}{d\sqrt{gd}}$$ (2.129)

with $\sigma_p = 2\pi f_p$.


The formulation of the interaction coefficient by Madsen and Sørensen (1993) is derived from their weakly nonlinear Boussinesq-wave model with improved wave dispersion and is given by

  $$R_{(m,p-m)} = -\frac{\left( k_m + k_{p-m} \right)^2\, \left ( \fr... ...( k_p d + 2\,B\,k^3_p\,d^3 - (B+\frac{1}{3})\,\sigma^2_p\,k_p\,d^2/g \right )}$$ (2.130)

where $c_l = \sigma_l/k_l$ and $k_l$ are the phase speed and the wave number, respectively, derived from the linear dispersion: $\sigma^2_l = g\,k_l\,\tanh(k_l\,d)$. Furthermore, $B = 1/15$.


Bredmose et al (2005) proposed another interaction coefficient with which the nonlinear wave transfer matches exactly with the second order Stokes theory. The associated expression reads

  $$R_{(m,p-m)} = \frac{N_{(m,p-m)}}{H_{mp}}$$ (2.131)

with

  $$N_{(m,p-m)} = -\frac{1}{2}\,\frac{g}{\sigma_m\,\sigma_{p-m}}\,\le... ...,\frac{\sigma^2_{mp}}{g} \,\left ( \sigma_m\,\sigma_{p-m} - \sigma^2_p\right )$$ (2.132)

and

  $$H_{mp} = \frac{\sigma^2_p - \sigma^2_{mp}}{k_p - k_{mp}}$$ (2.133)

where $k_{mp} = k_m+k_{p-m}$ and $\sigma^2_{mp} = g\,k_{mp}\,\tanh(k_{mp}\,d)$.


Akrish et al (2024) proposed the QuadWave model that preserves full dispersion of the Bredmose et al (2005) model but minimizes the error related to nonlinearity. The formulation of Bredmose et al (2005) is adapted through a parametrization which is obtained by means of an optimization process as described in Akrish et al (2024). The final experession reads

  $$R_{(m,p-m)} = W_{(m,p-m)}\,\frac{N_{(m,p-m)}}{H_{mp}}$$ (2.134)

where $W_{(m,p-m)}$ is the weight function defined as

  $$W_{(m,p-m)} = \exp{\Biggl [-\left( \frac{\chi}{\alpha_3} \right)^{\alpha_2} \Biggr]}$$ (2.135)

and

  $$\chi = \vert k_{mp}\vert\,d\,\left ( \frac{\vert k_{mp}\vert}{\vert k_p\vert} \right )^{\alpha_1}$$ (2.136)

and $\alpha_1$, $\alpha_2$ and $\alpha_3$ are the optimization parameters. This weight function optimizes the predictive accuracy associated with the nonlinear development of waves propagating through the coastal waters. This is obtained in SWAN using the following values: $\alpha_1 = 1$, $\alpha_2 = 0.4$ and $\alpha_3 = 5.5$. (Note that the original $\alpha_2-$value of 1.4 yields too much energy in the high-frequency part of the spectrum.)


DCTA model


Booij et al. (2009) proposed a heuristic triad formulation, in analogy with the quadruplet interaction, that accounts for both the generation of all super harmonics and the transition to a universal tail of $k^{-4/3}$. The original expression for the Distributed Collinear Triad Approximation (DCTA) is given by

  $$S_{\rm nl3} (\sigma_1) = \lambda \, \frac{\sin(-\beta)\,\tilde{k}... ...,2}k_2^p\,N(\sigma_2) - \sigma_1\,c_{g,1}k_1^p\,N(\sigma_1) \Bigr] \,d\sigma_2$$ (2.137)

for the quasi-resonance conditions $\sigma_3 = \vert\sigma_2 - \sigma_1\vert$ (note that the frequencies match but the wave numbers not). Here, $\lambda$ is a calibration coefficient that controls the magnitude of triad interactions, $\beta$ is the parametrized biphase, $\tilde{k} = \tilde{\sigma}/\sqrt{gd}$ with $\tilde{\sigma}$ the mean frequency as given by Eq. (2.67), $p = 4/3$ is a shape coefficient to force the high-frequency tail, and $\overline{k} = (k_1+k_2+k_3)/3$ is a characteristic wave number of the triad. Note that the factor $\tanh(\overline{k}d)/\overline{k}d$ in Eq. (2.137) accounts for the increasing resonance mismatch with increasing wave number (Booij et al., 2009).


However, recent study (Zijlema, 2022) has demonstrated that the following energy-flux conservative expression

  $$S_{\rm nl3} (\sigma_1) = \lambda \, c_{g,1}\frac{\sin(-\beta)\,\t... ...Bigl[ c_{g,2}k_2^p\,E(\sigma_2) - c_{g,1}k_1^p\,E(\sigma_1) \Bigr] \,d\sigma_2$$ (2.138)

yields improved prediction accuracy over the original one (2.137). This revised formulation is implemented in version 41.45.


In the above DCTA formulations, a collinear approximation is applied by which it is assumed that the primary contribution to the triad interactions arises from collinear interactions. Effectively, the triad source term for each directional bin is taken independently of other directional bins. An extension to noncollinear triad interactions is proposed by Benit and Reniers (2022) which includes a transfer reduction scaling based on the angle difference between two noncollinear interacting components. The final expression yields

$$S_{\rm nl3} (\sigma_1,\theta_1) = \lambda \, c_{g,1}\frac{\sin(-\beta)\,\tilde{k}^{2-p}}{\tilde{\sigma}^2\,d^2} \,\, \times$$

$$\int_{0}^{2\pi} \int_{0}^{\infty} \biggl (\frac{\tanh\,\overline{... ..._2,\theta_2) - c_{g,1}k_1^p\,E(\sigma_1,\theta_1) \Bigr] \,d\sigma_2\,d\theta_2$$

with $G(\Delta \theta_{nm})$ the transfer function of Sand (1982) and $\Delta \theta_{nm} = \theta_n - \theta_m$.


As a final note, the triad wave-wave interactions due to either LTA or DCTA are calculated only for $Ur \geq 0.1$.