Quadruplets

In this section two methods are described for the computation of nonlinear interactions at deep water. The first method is called the DIA method and is relatively crude in the approximation of the Boltzmann integral. The second one is called the XNL approach and is implemented in SWAN by G. Ph. van Vledder.


DIA


The quadruplet wave-wave interactions are computed with the Discrete Interaction Approximation (DIA) as proposed by Hasselmann et al. (1985). Their source code (slightly adapted by Tolman, personal communication, 1993) has been implemented in the SWAN model. In the DIA two quadruplet wave number configurations are considered, both with frequencies:

  $$\sigma_1 = \sigma_2 = \sigma\, , \sigma_3 = \sigma (1+\lambda) = \sigma^+ \, , \sigma_4 = \sigma (1 - \lambda) = \sigma^{-}$$ (2.76)

where $\lambda$ is a coefficient with a default value of 0.25. To satisfy the resonance conditions for the first quadruplet, the wave number vectors with frequency $\sigma_3$ and $\sigma_4$ lie at an angle of $\theta_3 = -11.48^o$ and $\theta_4 = 33.56^o$ to the angle of the wave number vectors with frequencies $\sigma_1$ and $\sigma_2$. The second quadruplet is the mirror image of the first quadruplet with relative angles of $\theta_3 = \theta^+ = 11.48^o$ and $\theta_4 = \theta^- = -33.56^o$. An example of this wave number configuration is shown in Figure 2.2. See Van Vledder et al. (2000) for further information about wave number configurations for arbitrary values of $\lambda$.

Figure 2.2: Wave number configuration for $\lambda$=0.25 and its position in a discrete frequency- direction spectrum (from Van Vledder et al., 2000).

Within this discrete interaction approximation, the source term $S_{\rm nl4}(\sigma,\theta)$ for the nonlinear transfer rate is given by

  $$S_{\rm nl4} (\sigma,\theta) = S^*_{\rm nl4} (\sigma,\theta) + S^{**}_{\rm nl4} (\sigma,\theta)$$ (2.77)

where $S_{\rm nl4}^{*}$ refers to the first quadruplet and $S_{\rm nl4}^{**}$ to the second quadruplet (the expressions for $S_{\rm nl4}^{**}$ are identical to those for $S_{\rm nl4}^{*}$ for the mirror directions). The DIA exchanges wave variance at all three wave number vectors involved in a quadruplet wave number configuration. The rate of change of wave variance due to the quadruplet interaction at the three frequency-direction bins can be written as

$$\left( \begin{array}{l} \delta S^*_{\rm nl4} (\sigma,\theta) \\ ... ...gma^+,\theta^+) \\ \delta S^*_{\rm nl4} (\sigma^-,\theta^-) \end{array}\right) = \left( \begin{array}{c} 2 \\ -1 \\ -1 \end{array}\right) C_{\rm nl4} (2\pi)^2 g^{-4} \left(\frac{\sigma}{2\pi} \right)^{11} \times$$

$$\left[E^2 (\sigma,\theta) \left\{\frac{E(\sigma^+,\theta^+)}{(1+\lambda)^4} + \frac{E(\sigma^{-},\theta^-)}{(1-\lambda)^4} \right\} \right.$$

$$\left. \left. -2\frac{E(\sigma,\theta) E(\sigma^+,\theta^+) E(\sigma^-,\theta^-)}{(1-\lambda^2)^4} \right. \right]$$ (2.78)

where $C_{nl4} = 3 \times 10^7$ by default. Eq. (2.78) conserves wave variance, momentum and action when the frequencies are geometrically distributed (as is the case in the SWAN model). The wave variance density at the frequency-direction bins $E(\sigma^+,\theta^+)$ and $E(\sigma^-,\theta^-)$ is obtained by bi-linear interpolation between the four surrounding frequency-direction bins. Similarly, the rate of change of variance density is distributed between the four surrounding bins using the same weights as used for the bi-linear interpolation.


In the DIA algorithm, Eq. (2.78) (and its mirror image) is applied to all spectral bins in a discrete frequency-direction spectrum. Figure 2.2 shows an example of one wave number configuration and its mirror image in a discrete spectrum. An extended spectral grid is applied to compute the interactions in the frequency range affected by the parametric spectral tail.


SWAN has an option to replace bi-linear interpolation of wave variance density using the nearest bin approach using a weight equal to 1.


Following the WAM group (WAMDI, 1988), the quadruplet interaction in shallow water with depth $d$ is obtained by multiplying the deep water nonlinear transfer rate by a scaling factor $R(k_pd)$:

  $$S_{\rm nl4}^{\rm finite~depth} = R(k_pd)S_{\rm nl4}^{\rm deep~water}$$ (2.79)

where $R$ is given by

  $$R(k_pd) = 1 + \frac{C_{sh1}}{k_pd} (1-C_{sh2}k_pd) e^{C_{sh3}k_pd}$$ (2.80)

in which $k_p$ is the peak wave number of the frequency spectrum. WAMDI (1988) proposes the following values of the coefficients: $C _{sh1}= 5.5$, $C _{sh2} = 5/6$ and $C _{sh3} = -5/4$. In the shallow water limit, i.e., $k_p \rightarrow 0$ the nonlinear transfer rate tends to infinity. Therefore, a lower limit of $k_p = 0.5$ is applied, resulting in a maximum value of $R(k_pd)=4.43$. To increase the model robustness in case of arbitrarily shaped spectra, the peak wave number $k_p$ is replaced by $k_p = 0.75 \tilde{k}$ (cf. Komen et al., 1994).


XNL (G. Ph. van Vledder)


The second method for calculating the nonlinear interactions in SWAN is the so-called Webb-Resio-Tracy method (WRT), which is based on the original six-dimensional Boltzmann integral formulation of Hasselmann (1962, 1963a,b), and additional considerations by Webb (1978), Tracy and Resio (1982) and Resio and Perrie (1991). A detailed description of the WRT method and its implementation in discrete spectral wave models like SWAN is given in Van Vledder (2006). An overview of computational methods for computing the exact nonlinear transfer rate is given in Benoit (2005).


The Boltzmann integral describes the rate of change of action density of a particular wave number due to resonant interactions between pairs of four wave numbers. To interact these wave numbers must satisfy the following resonance conditions

  $$\left . \begin{array}{ccc} \vec{k}_1 + \vec{k}_2 & = & \vec{k}_3... ... \sigma_2 & = & \sigma_3 + \sigma_4 \end{array} \:\:\: \right \rbrace \:\:\: .$$ (2.81)

The rate of change of action density $N_1$ at wave number $\vec{k}_1$ due to all quadruplet interactions involving $\vec{k}_1$ is given by

$$\frac{{\partial N_1 }}{{\partial t}} = \int\!\!\int\!\!\int G\left( \vec{k}_1 ,\vec{k}_2 ,\vec{k}_3 ,\ve... ...ight) \: \delta \left( {\sigma _1 + \sigma _2 - \sigma _3 - \sigma _4 } \right)$$

$$\hspace{7mm}\times \:\:\: \left[ {N_1 N_3 \left( {N_4 - N_2 } \ri... ...N_3 - N_1 } \right)} \right] \: d\vec{k}_2 \: d\vec{k}_3 \: d\vec{k}_4 \:\:\: ,$$ (2.82)

where the action density $N$ is defined in terms of the wave number vector $\vec{k}$, $N = N(\vec{k})$. The term $G$ is a complicated coupling coefficient for which an explicit expression has been given by Herterich and Hasselmann (1980). In the WRT method a number of transformations are made to remove the delta functions. A key element in the WRT method is to consider the integration space for each ( $\vec{k}_1 ,\vec{k}_3$) combination

  $${{\partial N_1 } \over {\partial t}} = 2\int T\left( {\vec{k}_1 ,\vec{k}_3 } \right) \: d\vec{k}_3 \:\:\: ,$$ (2.83)

in which the function $T$ is given by

$$T \left( \vec{k}_1 ,\vec{k}_3 \right) = \int\!\!\int G\left( \vec{k}_1 ,\vec{k}_2, \vec{k}_3 ,\vec{k}_4 \right) \: \delta \left( \vec{k}_1 + \vec{k}_2 - \vec{k}_3 - \vec{k}_4 \right)$$

$$\times \delta \left( {\sigma _1 + \sigma _2 - \sigma _3 - \sigma _4 } \right) \: \theta \left ( \vec{k}_1 , \vec{k}_3 , \vec{k}_4 \right )$$

$$\times \left [ N_1 N_3 \left ( N_4 - N_2 \right ) + N_2 N_4 \left ( N_3 - N_1 \right ) \right ] \: d\vec{k}_2 \: d\vec{k}_4 \:\:\: ,$$ (2.84)

in which

  $$\theta \left( {\vec{k}_1 ,\vec{k}_3 ,\vec{k}_4 } \right) = \left\... ... > \left\vert {\vec{k}_1 - \vec{k}_4 } \right\vert} \\ \end{array} } \right.$$ (2.85)

The delta functions in Eq. (2.84) determine a region in wave number space along which the integration should be carried out. The function $\theta$ determines a section of the integral which is not defined due to the assumption that $\vec{k}_1$ is closer to $\vec{k}_3$ than $\vec{k}_2$. The crux of the Webb method consists of using a local co-ordinate system along a so-named locus, that is, the path in $\vec{k}$ space that satisfies the resonance conditions for a given combination of $\vec{k}_1$ and $\vec{k}_3$. To that end the $(k_{x},k_{y})$ co-ordinate system is replaced by a $(s,n)$ co-ordinate system, where $s$ ($n$) is the tangential (normal) direction along the locus. After some transformations the transfer integral can then be written as a closed line integral along the closed locus

$$T \left( \vec{k}_1 ,\vec{k}_3 \right) = \oint G \: J \: \theta ( \vec{k}_1 , \vec{k}_3 , \vec{k}_4 )$$

$$\times \left [ N_1 N_3 \left ( N_4 - N_2 \right ) + N_2 N_4 \left ( N_3 - N_1 \right ) \right ] \: d s \:\:\: ,$$ (2.86)

in which $G$ is the coupling coefficient and $J$ is the Jacobian term of a function representing the resonance conditions. The Jacobian term is a function of the group velocities of interacting wave numbers

  $$J = \vert \vec{c}_{g,2} - \vec{c}_{g,4}\vert^{-1}$$ (2.87)

Numerically, the Boltzmann integral is computed as the finite sum of many line integrals $T$ for all discrete combinations of $\vec{k}_1$ and $\vec{k}_3$. The line integral (2.86) is solved by dividing the locus in typically 40 pieces, such that its discretized version is given by

  $$T\left( {\vec{k}_1 ,\vec{k}_3 } \right) \approx \sum_{i=1}^{n_s } G(s_i) J(s_i) P(s_i) \: \Delta s_i \:\:\: ,$$ (2.88)

in which $P(s_i)$ is the product term for a given point on the locus, $n_s$ is the number of segments, $s_i$ is the discrete co-ordinate along the locus, and $\Delta s_i$ is the stepsize. Finally, the rate of change for a given wave number $\vec{k}_1$ is given by

  $$\frac{\partial N(\vec{k}_1)}{\partial t} \approx \sum_{i_{k3} = 1... ...T(\vec{k}_1,\vec{k}_3) \: \Delta k_{i_{k3}} \: \Delta \theta_{i_{\theta 3} } }$$ (2.89)

where $n_k$ and $n_\theta$ are the discrete number of wave numbers and directions in the computational spectral grid, respectively. Note that although the spectrum is defined in terms of the vector wave number $\vec{k}$, the computational grid in a wave model is more conveniently defined in terms of the absolute wave number and wave direction ($k,\theta$) to assure directional isotropy of the calculations. Taking all wave numbers $\vec{k}_1$ into account produces the complete source term due to nonlinear quadruplet wave-wave interactions. Details of the computation of a locus for a given combination of the wave numbers $\vec{k}_1$ and $\vec{k}_3$ can be found in Van Vledder (2006).


It is noted that these exact interaction calculations are extremely expensive, typically requiring $10^3$ to $10^4$ times more computational effort than the DIA. Presently, these calculations can therefore only be made for highly idealized test cases involving a limited spatial grid.


The nonlinear interactions according to the WRT method have been implemented in SWAN using portable subroutines. In this implementation, the computational grid of the WRT method is based to the discrete spectral grid of SWAN. The WRT method uses a $(\vec{k},\theta)$ grid which is based on the $(\sigma,\theta)$ grid of SWAN. In addition, the WRT routines inherit the power of the parametric spectral tail as in the DIA. Choosing a higher resolution than the computational grid of SWAN for computing the nonlinear interactions is possible in theory, but this does not improve the results and is therefore not implemented.


Because nonlinear quadruplet wave-wave interactions at high frequencies are important, it is recommended to choose the maximum frequency of the wave model about six times the peak frequency of the spectra that are expected to occur in a wave model run. Note that this is important as the spectral grid determines the range of integration in Eq. (2.89). The recommended number of frequencies is about 40, with a frequency increment factor 1.07. The recommended directional resolution for computing the nonlinear interactions is about $10^\circ$. For specific purposes other resolutions may be used, and some testing with other resolutions may be needed.


An important feature of most algorithms for the evaluation of the Boltzmann integral is that the integration space can be pre-computed. In the initialization phase of the wave model the integration space, consisting of the discretized paths of all loci, together with the interaction coefficients and Jacobians, are computed and stored in a binary data file. For each discrete water depth such a data file is generated and stored in the work directory. The names of these data files consist of a keyword, "xnl4v5", followed by the keyword "xxxxx", with xxxxx the water depth in a certain unit (meters by default), or 99999 for deep water. The extension of the binary data file is "bqf" (of Binary Quadruplet File). If a BQF file exists, the program checks if this BQF file has been generated with the proper spectral grid. If this is not the case, a new BQF file is generated and the existing BQF file is overwritten. During a wave model run with various depths, the optimal BQF is used, by looking at the 'nearest' water depth $d_N$ for which a valid BQF file has been generated. In addition, the result is rescaled using the DIA scaling (2.80) according to

  $$S_{\rm nl4}^{d} = S_{\rm nl4}^{d_N} \frac{R(k_pd)}{R(k_pd_N)} \:\:\: .$$ (2.90)