Implementation of QC approximation

One observes that the phase-space equation (2.225)$-$(2.230) is formulated in $(\vec{x},\vec{k})-$space, whereas SWAN operates in frequency-direction space. In particular due to the QC scattering term that involves a Fourier transform from physical space to a wave number space, one should stick to the $\vec{k}-$space in solving the governing equation for the statistics of the wave field (as represented by the Wigner distribution). Examples of such an approach are presented in Smit et al. (2015a) and Akrish et al. (2020).


Here, we take a different route by applying a transformation to $(\sigma,\theta)-$space. We first define the Wigner distribution in the frequency-direction space and subsequently transform it to the $\vec{k}-$space. Next, we compute the scattering term in the phase space using the transformed Wigner spectrum, and lastly, the result is transformed back to the frequency-direction space.


For the purpose of the implementation in SWAN we thus consider the following equation

  $$\frac{\partial W}{\partial t} + \nabla_{\vec{x}} \cdot [({\vec{c}}_g + \vec{u}) W] = J^{-1}\,S_{\rm qc}$$ (3.86)

with $W(\vec{x},t;\sigma,\theta)$ the Wigner distribution in the appropriate spectral space and $J = \vert\vec{c}_g\vert/\vert\vec{k}\vert$ the Jacobian. Note that for evaluating the source term $S_{\rm qc}(\vec{x},\vec{k})$ we thus priorly need to compute $W(\vec{x},t;\vec{k}) = J\,W(\vec{x},t;\sigma,\theta)$.


It is important to note that the Wigner distribution extends the notion of the action density spectrum $N(\vec{x},t;\sigma,\theta)$ by means of the introduction of the cross correlations. It basically implies conservation in the presence of ambient current. So for computing the spectral moments and, in turn, the integral parameters (e.g. significant wave height and mean period) one must consider the product $\sigma\,W$ with which the integration over spectral space can be performed. Additionally, recall that $W$ can take on negative values. This implies, for example, that the conservative elimination approach as outlined in Section 3.2.4 should not be applied because otherwise the moments will be wrongly computed. Note that by definition, $\int \int\, W(\sigma,\theta)d\sigma d\theta >0$.


There are, however, a few issues that require attention. Below, we discuss each of these issues in some more depth.


There may arise a problem of mapping the cross-variance contributions to the frequency domain using the dispersion relation (2.3) (assuming for the moment no currents). Generally, such contributions do not obey the dispersion relation. This is relevant because of the transformation from $\vec{k}-$space to $(\sigma,\theta)-$space, and vice versa. For instance, two interacting waves $\vec{k}_1$ and $\vec{k}_2$ having the same frequency $\sigma(\vert\vec{k}_1\vert) = \sigma(\vert\vec{k}_2\vert) = \sigma$ create a cross contribution at frequency $\sigma(\vert(\vec{k}_1+\vec{k}_2)/2\vert) \neq \sigma$, with the result that a mapping error occur when one is transforming $W(\vec{k})$ to $W(\sigma,\theta)$. Although the cross correlations may thus be aliased at offset frequencies, we expect that such aliasing errors would prove to be small enough to be accepted (see also Appendix C of Smit et al., 2015b).


Another issue concerns the inclusion of the statistical effects of wave refraction and Doppler shifting induced by medium variations in the source term $S_{\rm qc}$. As a result, the modelling of these wave processes do not appear explicitly as conventional transport terms in Eq. (3.86) (viz. the third and fourth terms of Eq. (2.16)). There are two reasons for this approach. First, it circumvents the stiffness of the equation that otherwise would occur owing to the relatively rapid variations in the medium, thus leading to unnecessary and excessive resolution in the spectral space (Akrish et al., 2020). Second, by means of the local plane approximation of Smit et al. (2015a) the refraction term $\nabla_{\vec{x}}\,\omega\cdot\nabla_{\vec{k}}\,W$ could have been extracted out of the integral in Eq. (2.230), which is, however, formulated in $\vec{k}-$space rather than in $(\sigma,\theta)-$space. This would mean an additional error-prone coordinate mapping through the Jacobian transformation (and interpolation).


The quasi-coherent modelling of wave processes accounting for depth-induced wave breaking (Smit et al., 2015b) and nonlinear (triad) interactions (Smit and Janssen, 2016) are not considered in the current QC framework. These developments are in progress.


With respect to the numerical implementation of Eq. (3.86), the left hand side is approximated in the usual way as outlined in Sections 3.2 and 3.3 or, in case of the use of unstructured meshes, in Section 8.3. However, because of the interaction with small-scale medium variations, the Wigner distribution can narrow and broaden rapidly during the propagation in the physical space. Hence, it is recommended to apply at least second order accuracy of spatial discretization to resolve the shoaling term. In case of structured grids the SORDUP scheme (or the Stelling and Leendertse scheme for non-stationary runs) should then be selected (see Section 3.2.1). In the same vein, a more accurate unstructured mesh discretization would be preferred over the currently implemented first order upwind finite difference scheme. The work on this extension is in progress.


The rest of this section addresses the details of implementation issues with respect to the right hand side of Eq. (3.86). In this regard, we first define three discrete grids in the $\vec{k}-$, $\vec{x}^\prime-$ and $\vec{q}-$spaces, respectively, and after this we provide an outline of the QC implementation.


Evaluation of the QC scattering source term calls for a grid in $\vec{k}-$space in which the effects of wave interference are resolved. The origin and the size of its wave number domain $D_{\vec{k}} = \{(k_x,k_y)\,\vert\,k_{x,0} \leq k_x \leq k_{x,{\rm max}}, \, k_{y,0} \leq k_y \leq k_{y,{\rm max}}\}$ is determined by the transformation of the user-defined spectral $(\sigma,\theta)-$domain (commonly defined as a directional sector). Apart from this, the determination of the wave number grid spacing is a crucial factor as explained in the following.


The coherent interference patterns in wave fields depend mainly on the nature of the medium variations and the spectral width of the Wigner distribution. In particular, narrow-banded wave groups entering shallow water with an uneven bottom can cause substantial variations in the wave statistics over many wave lengths and, in turn, the wave field remains thus correlated over such distances. Hence, to determine the importance of the coherent effects and also the resolution of the wave number grid, knowledge of the spectral width of the incident waves is required. In SWAN, it is computed as the standard deviation of the incoming spectrum of the action density $N(\vec{k})$. The grid size as denoted by $\Delta k$ is then set to this standard deviation. In relation to this, the coherent length scale of the correlated wave field is defined as $L_{\rm c} = 2\pi/\Delta k$ (Smit et al., 2015a). Note that the two-dimensional $\vec{k}$ grid thus obtained is equidistant in both directions.


The computation of the QC source term also requires a discrete Fourier transform of $\Delta \omega(\vec{x}^\prime)$ to obtain $\Delta {\hat \omega}(\vec{q})$. To this end, the so-called coherent region in $\vec{x}^\prime-$space is defined as a square with the origin in the center and, owing to its compact support, its side length is half the correlation length, thus $D_{\vec{x}^\prime} = \{(x^\prime,y^\prime)\,\vert\,-\frac{1}{4}L_c \leq x^\prime \leq \frac{1}{4}L_c, \,-\frac{1}{4}L_c \leq y^\prime \leq \frac{1}{4}L_c\}$. Furthermore, this domain is equipped with a grid consisting of $(M+1)\times(M+1)$ equally spaced points separated by a distance $\Delta x^\prime = \Delta y^\prime$ in each direction. The specification of $M$ will be discussed below.


We finally introduce the medium wave number domain $D_{\vec{q}}$ over which the integral of Eq. (2.230) is to be performed. This domain is conjugated to $D_{\vec{x}^\prime}$ and, since the Nyquist wave length is $2\Delta x^\prime$, it is truncated to $\{(q_x,q_y)\,\vert\,-q_{\rm max} \leq q_x \leq q_{\rm max}, \,-q_{\rm max} \leq q_y \leq q_{\rm max}\}$ with $q_{\rm max} = \pi/\Delta x^\prime$. This maximum scattering wave number can be set directly by the user or, alternatively, prescribed as a multiple of the peak (or mean) wave number at the wavemaker $\vert\vec{k}_p\vert$, $q_{\rm max} = \alpha\,\vert\vec{k}_p\vert$ with $\alpha$ usually 0.5, 1 or 2.


The domain $D_{\vec{q}}$ is subsequently discretized into uniform grid of $N\times N$ bins with grid size $\Delta q$. We select $\Delta q = 2\,\Delta k$ so as to avoid interpolation errors in computing the convolution integral of Eq. (2.230). Consequently, $N = q_{\rm max} / \Delta k$. From the above we have that $M = \frac{1}{2}\,L_{\rm c} / \Delta x^\prime = L_{\rm c}\,q_{\rm max} / 2\pi = L_{\rm c}\,N\,\Delta k / 2\pi = N$. So, there is a one-to-one mapping between the grid of $D_{\vec{x}^\prime}$ and the grid of $D_{\vec{q}}$, thus requiring no interpolation.


With the above definitions of discrete grids and their resolution, the discrete form of the convolution integral and subsequently the QC scattering term (2.230) can be readily obtained as follows

$$S_{\rm qc} (x,y,k_x,k_y) = -{\rm i}\sum_{q_x = -q_{\rm max}}^{q_{\rm max}} \, \sum_{q_y = -q... ...}}\, \biggl [\Delta {\hat \sigma} + \vec{k}\cdot \Delta {\hat{\vec{u}}} \biggr.$$

$$\biggl. -\frac{{\rm i}}{2}\Bigl( \Delta {\hat{c_g}}\, \frac{\vec{... ...bla_{\vec{x}} \biggr]\, W(x,y,k_x-\frac{q_x}{2},k_y-\frac{q_y}{2}) + \rm {c.c.}$$ (3.87)

with $\Delta {\hat \sigma}$, $\Delta {\hat{c_g}}$ and $\Delta {\hat{\vec{u}}}$ the discrete Fourier transform of, respectively,

$$\Delta \sigma = {\cal W}(x^\prime)\,{\cal W}(y^\prime) \left [ \,\sigma(x+x^\prime,y+y^\prime,k_x,k_y) - \sigma(x,y,k_x,k_y) \,\right ]$$

$$\Delta c_g = {\cal W}(x^\prime)\,{\cal W}(y^\prime) \left [ \,c_g(x+x^\prime,y+y^\prime,k_x,k_y) - c_g(x,y,k_x,k_y) \,\right ]$$

$$\Delta \vec{u} = {\cal W}(x^\prime)\,{\cal W}(y^\prime) \left [ \,\vec{u}(x+x^\prime,y+y^\prime) - \vec{u}(x,y) \,\right ]$$

Here, ${\cal W}(z)$ is a window function which tapers the function to be transformed near the boundaries of the $z-$interval (and so to reduce or possibly remove any jump discontinuity of the outcome). A standard Tukey window is utilized and is given by

  $${\cal W}(z) = \left\{ \begin{array}{ll} \frac{1}{2}\left\{ 1 +... ...ight ]\right ) \right\}, & (1-\gamma)\,\ell < z \leq \ell \end{array} \right.$$ (3.88)

with $\ell = L_{\rm c}/2$ the side of the square coherent zone centered at $\vec{x}$ and $0 \leq \gamma \leq 0.5$ a dimensionless taper parameter. In SWAN, it is set hardcoded to $\gamma = 0.2$.


As a final note on the implementation of the convolution sum (3.87), the differential operator $\nabla_{\vec{x}}\,W$ is discretized by means of second order central differences on structured grids. In case of an unstructured mesh, the Green-Gauss theorem with the assumption of a constant gradient over the centroid dual volume is then applied; viz. Eq. (8.35).