Modelling statistical wave interferences over shear currents

Abstract

Wave forecasting in ocean and coastal waters commonly relies on spectral models based on the spectral action balance equation. These models assume that different wave components are statistically independent and as a consequence cannot resolve wave interference due to statistical correlation between crossing waves, as may be found in, for instance, a focal zone. This study proposes a statistical model for the evolution of wave fields over non-uniform currents and bathymetry that retains the information on the correlation between different wave components. To this end, the quasi-coherent model (Smit & Janssen, J. Phys. Oceanogr., vol. 43, 2013, pp. 1741–1758) is extended to allow for wave–current interactions. The outcome is a generalized action balance model that predicts the evolution of the wave statistics over variable media, while preserving the effect of wave interferences. Two classical examples of wave–current interaction are considered to demonstrate the statistical contribution of wave interferences: (1) swell field propagation over a jet-like current and (2) the interaction of swell waves with a vortex ring. In both examples cross-correlation terms lead to development of prominent interference structures, which significantly change the wave statistics. Comparison with results of the SWAN model demonstrates that retention of cross-correlation terms is essential for accurate prediction of wave statistics in shear-current-induced focal zones.

1 Introduction

Wind-generated waves play an important role in the dynamics of oceanic and coastal waters. In the upper ocean, surface waves can force large-scale circulations (e.g. Craik & Leibovich 1976), whereas near the shore they can drive alongshore currents (e.g. Bowen 1969; Longuet-Higgins 1970; Reniers and Battjes 1997; Ruessink et al. 2001), return flow (e.g. Dyhr-Nielsen & Sørensen 1970; Stive and De Vriend 1995) and associated sediment transport processes (e.g. Deigaard et al. 1992; Van Rijn 1993). Furthermore, waves control shipping operations and associated downtime as well as coastal safety through beach and dune erosion and potential inundation (e.g. Vellinga 1982; Roelvink et al. 2009).

The common approach to predicting statistical parameters of wind waves is via operational (phase-averaged) wave models, e.g. WAM model (Group 1988), WAVEWATCH model (Tolman 1991) and SWAN model (Booij et al. 1999). These models solve numerically the so-called spectral action balance equation that can be written in the following form:

(1.1)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}N+\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }(\boldsymbol{C}_{x}N)+\unicode[STIX]{x1D735}_{\boldsymbol{k}}\boldsymbol{\cdot }(\boldsymbol{C}_{k}N)=S,\end{eqnarray}$$

where $N$ represents the spectrum of the action density, being equal to the spectrum of the energy density, $E$, divided by the intrinsic frequency, $\unicode[STIX]{x1D70E}$. The propagation part, on the left-hand side, describes the kinematic behaviour of the field as it propagates through slowly varying current, $\boldsymbol{U}$, and bathymetry, with propagation velocities $\boldsymbol{C}_{k}$ and $\boldsymbol{C}_{x}$ over wavenumber space, $\boldsymbol{k}=(k_{1},k_{2})$, and physical space, $\boldsymbol{x}=(x_{1},x_{2})$, respectively. On the right-hand side, the equation is forced by source terms, $S$, to account for processes of wave generation (by wind), dissipation (e.g. due to whitecapping) and wave–wave interactions.

The statistical assumptions underlying the derivation of (1.1) are that the wave field can be regarded as Gaussian and quasi-homogeneous. The former suggests that the field is completely defined by its correlation function (assuming a zero-mean field), while the latter proposes that the correlation between any two distinct wave components equals zero. Based on these assumptions, variation of the field statistics is governed completely by variations of the wave variances (which are represented by $N$), as indeed described by (1.1).

In most circumstances at sea, the parameters of the wave field (e.g. wave amplitudes) are evolving slowly over spatial scales of $O$(10–100 km) due to the action of wind, slow medium changes and weak nonlinearity. Under these conditions, the assumption of quasi-homogeneity is easily met, and (1.1) remains valid. However, there might be situations where the field encounters medium variability over much smaller scales $O$(0.1–1 km). Such situations can occur quite frequently in coastal regions, where currents and bathymetry can vary rapidly (e.g. Chen et al. 1999; Ardhuin et al. 2003). Furthermore, following recent studies (e.g. Poje et al. 2014; McWilliams 2016), they may also occur in the open ocean over small-scale currents (e.g. submesoscale currents). Physically, in these situations, waves are rapidly scattered into multiple directions, and consequently can form focal zones which give rise to wave interferences. Well-known examples of such wave–media interactions are given by the evolution of waves over a submerged shoal (e.g. Vincent & Briggs 1989) or over a vortex ring (e.g. Yoon & Liu 1989). Statistically, the interference effects that arise in such cases are described by cross-correlations between different wave components of the scattered field and may result in significant and rapid variations of the wave statistics (Janssen et al. 2008; Smit & Janssen 2013; Smit et al. 2015a). The quasi-homogeneous assumption excludes the contribution of the cross-correlation terms, and therefore equation (1.1) cannot describe the effect of wave interferences arising in interactions between waves and rapidly varying media.

The ability to account for the effect of wave interferences in these situations is important, since they can alter dramatically the spatial distributions of wave parameters (e.g. the significant wave height), which serve as input for numerous applications in coastal zones. In addition, through the interaction of waves with small-scale ocean currents, generated interference structures may also introduce leading-order statistical contributions for applications in the open ocean. For example, they may contribute to changes driven by waves of submesoscale currents (McWilliams 2018), or the interpretation of noise obtained (due to the presence of waves) in measurements of the sea surface, revealing the evolution of small-scale circulations (e.g. Ardhuin et al. 2017), and they may also enhance and alter the spatial distribution of extreme elevations in energetic focal regions (e.g. Metzger et al. 2014; Fedele et al. 2016).

In order to take into account the statistical effect of wave interference, Smit & Janssen (2013) and Smit et al. (2015a) have recently developed an evolution equation that allows for the generation and evolution of correlations between different wave components when interacting over small-scale bathymetry changes. This newly developed stochastic model is called the quasi-coherent model (QCM). The main aim of the present study is to extend the capabilities of the QCM so it can handle the interaction between waves and ambient currents. The derivation of the extended QCM is detailed in § 2. The model is verified in § 3 through the problem of interaction between a swell field and a jet-like current (e.g. Janssen and Herbers 2009). Then, the model is used to study the statistical mechanism for the generation of wave interferences in § 4, through the classical problem of interaction between swell waves and a vortex ring (e.g. Yoon & Liu 1989). Finally, conclusions are drawn in § 5.

2 Stochastic model for linear waves over varying current and bathymetry

Generally speaking, stochastic wave models are derived based on deterministic equations that physically describe the evolution of wave fields. This approach of deriving a stochastic formulation is also adopted here. Therefore, the derivation starts with a physical description of the wave field which is effectively represented by the so-called action variable. Section 2.1 introduces the definition of the action variable and its governing equation. As discussed in § 2.2, the second-order statistics of the wave field, including the statistics of wave interferences, are fully described through the correlation function or the spectral distribution function of the action variable. These starting points are used in § 2.3 to formulate a stochastic model that takes into account the generation and transportation of wave interference contributions. Finally, the numerical implementation of the model and an overview of the considered simulations are described in § 2.4 and § 2.5, respectively.

2.1 The action variable and its evolution equation

The formulation starts by considering the evolution of a random linear wave field through a variable medium that can be represented by its surface potential and surface elevation, $\unicode[STIX]{x1D719}(\boldsymbol{x},t)$ and $\unicode[STIX]{x1D702}(\boldsymbol{x},t)$. It is assumed that the medium changes slowly so that the ratio, $\unicode[STIX]{x1D716}=L/L_{m}$, between the characteristic wavelength, $L$, and the characteristic length scale of medium variation, $L_{m}$, is small ($\unicode[STIX]{x1D716}\ll 1$). Accordingly, the field can locally be approximated as a summation of plane waves with slowly varying phase and amplitude, which to leading order in $\unicode[STIX]{x1D716}$ obey to the following general dispersion relation (e.g. Dingemans 1997):

(2.1)$$\begin{eqnarray}\unicode[STIX]{x1D714}=\boldsymbol{U}\boldsymbol{\cdot }\boldsymbol{k}+\unicode[STIX]{x1D70E}.\end{eqnarray}$$

Variations in the medium are introduced by the ambient current, $\boldsymbol{U}(\boldsymbol{x})$, and by the water depth, $h(\boldsymbol{x})$. Using the medium information and the definition of the intrinsic frequency, $\unicode[STIX]{x1D70E}(\boldsymbol{x},\boldsymbol{k})=\sqrt{|\boldsymbol{k}|g\tanh (|\boldsymbol{k}|h)}$, the value of the absolute frequency, $\unicode[STIX]{x1D714}$, is obtained through (2.1), where $|\boldsymbol{k}|$ is the magnitude of the local wavenumber, defined as $|\boldsymbol{k}|=\sqrt{k_{1}^{2}+k_{2}^{2}}$, and $g$ is the gravitational acceleration. Finally, from the statistical point of view, the field is assumed to be zero-mean, Gaussian and quasi-stationary.

Under this statistical and physical framework, it will be convenient to use the so-called action variable (e.g. Besieris & Tappert 1976; Krasitskii 1994), $\unicode[STIX]{x1D713}$, which is defined as

(2.2)$$\begin{eqnarray}\unicode[STIX]{x1D713}={\displaystyle \frac{1}{\sqrt{2g}}}[g\,a^{-1}\unicode[STIX]{x1D702}+\text{i}a\,\unicode[STIX]{x1D719}],\end{eqnarray}$$

where $a(\boldsymbol{x},-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}})$ is a pseudo-differential operator that is associated with the symbol $a(\boldsymbol{x},\boldsymbol{k})=\sqrt{\unicode[STIX]{x1D70E}(\boldsymbol{x},\boldsymbol{ k})}$ (see detailed definition of this operator in appendix A).

The convenience of working with the action variable, $\unicode[STIX]{x1D713}$, becomes significant in the formulation of the second-order statistics of the field, since, following its definition, second-order statistical functions of $\unicode[STIX]{x1D713}$ (e.g. the correlation function) are inherently related to the definition of the wave action (Bretherton & Garrett 1968). As a consequence, the action variable, $\unicode[STIX]{x1D713}$, is intimately related to the mean action density and the mean energy density through the following expressions:

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\langle |\unicode[STIX]{x1D713}|^{2}\rangle =m_{0}/\unicode[STIX]{x1D70E}+O(\unicode[STIX]{x1D716}), & \displaystyle\end{eqnarray}$$

(2.4)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70C}\langle |a\unicode[STIX]{x1D713}|^{2}\rangle =m_{0}+O(\unicode[STIX]{x1D716}), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}$ is the water mass density and the angular brackets, $\langle \cdots \rangle$, should be read as an ensemble average. The variable $m_{0}$ provides a leading-order estimation (in $\unicode[STIX]{x1D716}$) of the mean energy density (also known as the zero-order moment of the spectral energy density) and it is defined as follows:

(2.5)$$\begin{eqnarray}m_{0}=\unicode[STIX]{x1D70C}\left\langle {\displaystyle \frac{1}{2}}g\unicode[STIX]{x1D702}_{0}^{2}+{\displaystyle \frac{1}{2g}}(\unicode[STIX]{x1D70E}\unicode[STIX]{x1D719})_{0}^{2}\right\rangle ,\end{eqnarray}$$

where now (in (2.3) and (2.5)) $\unicode[STIX]{x1D70E}(\boldsymbol{x},-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}})$ represents a pseudo-differential operator that is associated with the intrinsic frequency, $\unicode[STIX]{x1D70E}(\boldsymbol{x},\boldsymbol{k})$, and the subscript $0$ indicates $O(1)$ terms (refer to appendix A for the definition of $\unicode[STIX]{x1D70E}(\boldsymbol{x},-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}})$ and its leading-order operation, e.g. $(\unicode[STIX]{x1D70E}\unicode[STIX]{x1D719})_{0}$). Further details explaining why the expression in (2.5) defines the leading-order estimation of the mean energy density are given in appendix B.

An additional motivation for the definition of $\unicode[STIX]{x1D713}$ (see (2.2)) is that, under the physical assumptions made here, the governing equation of the wave field can be written as (Besieris & Tappert 1976; Besieris 1985)

(2.6)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}=-\text{i}\unicode[STIX]{x1D714}(\boldsymbol{x},-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}})\unicode[STIX]{x1D713},\end{eqnarray}$$

where $\unicode[STIX]{x1D714}(\boldsymbol{x},-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}})$ is a pseudo-differential operator that is associated with the dispersion relation, $\unicode[STIX]{x1D714}(\boldsymbol{x},\boldsymbol{k})$ (see appendix A). This equation form is convenient since it can be transformed directly into an evolution equation of the correlation function, which under the assumption of Gaussian statistics provides a complete statistical description of the wave field. A verification of this equation for homogeneous and weakly inhomogeneous media is described in appendix B. For homogeneous media, equation (2.6) exactly describes the evolution of the considered linear field. For weakly inhomogeneous media, equation (2.6) reduces for each wave component to the local dispersion relation, equation (2.1) (or the eikonal equation, which governs the evolution of the wavenumber) at leading order, and the well-known transport equation for the mean action density, $\langle |\unicode[STIX]{x1D713}|^{2}\rangle$ at $O(\unicode[STIX]{x1D716})$. This indicates that at $O(\unicode[STIX]{x1D716})$, equation (2.6) provides the correct representation of the evolution of the field.

To summarize, the formulation presented here considers a random, linear and slowly varying wave field, which is concisely represented by the action variable, $\unicode[STIX]{x1D713}$. The definition of this action variable introduces convenient properties which will eventually lead to a derivation of a generalized action balance equation that accounts for the effect of wave interferences. As a first step on this path, the next subsection aims to demonstrate that the statistical information about wave interferences is naturally included in the representative second-order statistical functions (i.e. the correlation function and the Wigner distribution).

2.2 Second-order statistics

Following the statistical assumptions for the surface variables, $\unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D719}$, and following the linearity of the definition (2.2), the action variable $\unicode[STIX]{x1D713}(\boldsymbol{x},t)$ is said to be a zero-mean, complex Gaussian and quasi-stationary field (e.g. Soong 1973). The statistics of such a random field are defined completely by the following correlation function:

(2.7)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}(\boldsymbol{x},\boldsymbol{x}^{\prime },t)=\langle \unicode[STIX]{x1D713}(\boldsymbol{x}+\boldsymbol{x}^{\prime }/2,t)\unicode[STIX]{x1D713}^{\ast }(\boldsymbol{x}-\boldsymbol{x}^{\prime }/2,t)\rangle .\end{eqnarray}$$

The statistical information carried by the correlation function is better seen using its spectral form, written as

(2.8)$$\begin{eqnarray}\unicode[STIX]{x1D6E4}(\boldsymbol{x},\boldsymbol{x}^{\prime },t)=\int \text{d}\boldsymbol{k}\exp (\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}^{\prime })\int \hat{\unicode[STIX]{x1D6E4}}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)\exp (\text{i}\boldsymbol{k}^{\prime }\boldsymbol{\cdot }\boldsymbol{x})\,\text{d}\boldsymbol{k}^{\prime },\end{eqnarray}$$

where $\boldsymbol{k}$ and $\boldsymbol{k}^{\prime }$ are defined as the average and difference of two interacting wavenumbers, namely $\boldsymbol{k}=(\boldsymbol{k}_{1}+\boldsymbol{k}_{2})/2$ and $\boldsymbol{k}^{\prime }=\boldsymbol{k}_{1}-\boldsymbol{k}_{2}$. In addition $\hat{\unicode[STIX]{x1D6E4}}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)$ is defined as $\hat{\unicode[STIX]{x1D6E4}}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)=\langle \hat{\unicode[STIX]{x1D713}}(\boldsymbol{k}+\boldsymbol{k}^{\prime }/2,t)\hat{\unicode[STIX]{x1D713}}^{\ast }(\boldsymbol{k}-\boldsymbol{k}^{\prime }/2,t)\rangle$. The expression above reveals the spectral content of the correlation function. It shows that, in general, $\unicode[STIX]{x1D6E4}$ oscillates with a wavenumber difference $\boldsymbol{k}^{\prime }$ over the space $\boldsymbol{x}$. Such an oscillation occurs when wave components with two different wavenumbers are statistically correlated, and thus creating a spatially dependent pattern of wave interference.

The assumption that the wave field is quasi-homogeneous trims the spectral information provided by $\hat{\unicode[STIX]{x1D6E4}}$ with respect to $\boldsymbol{k}^{\prime }$ and accounts only for a narrow window around $\boldsymbol{k}^{\prime }=0$, which consists of the components that characterize the slow changes of the medium. Therefore, under this assumption, the spectrum obtained by the Fourier transform of $\unicode[STIX]{x1D6E4}$ from $\boldsymbol{x}^{\prime }$ to $\boldsymbol{k}$ only allows for slow variations of the variance terms of the field. This spectrum is the conventional action density spectrum, $N(\boldsymbol{x},\boldsymbol{k},t)$. In this study, however, statistical inhomogeneity of the wave field is taken into account by considering the full spectrum provided by $\hat{\unicode[STIX]{x1D6E4}}$ with respect to $\boldsymbol{k}^{\prime }$. In this general case, the corresponding spectral representation of wave action follows the definition of the Wigner distribution, $W(\boldsymbol{x},\boldsymbol{k},t)$:

(2.9)$$\begin{eqnarray}W(\boldsymbol{x},\boldsymbol{k},t)=\int \unicode[STIX]{x1D6E4}(\boldsymbol{x},\boldsymbol{x}^{\prime },t)\exp (-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}^{\prime })\,\text{d}\boldsymbol{x}^{\prime }.\end{eqnarray}$$

Therefore, the Wigner distribution of $\unicode[STIX]{x1D713}$ captures the same information as the correlation function and basically generalizes the concept of the action density spectrum by including the cross-correlation terms that correspond to wave interferences (also see, e.g. Hlawatsch & Flandrin (1997)). As such, the Wigner distribution provides a complete spectral description of the second-order statistics of the field. Finally note that, as implied by (2.9), the zero-order moment of $W$ equals the variance of $\unicode[STIX]{x1D713}$, and therefore following (2.3) gives a leading-order evaluation of the mean action density.

Practically speaking, one would eventually be interested in certain field parameters (e.g. characteristic wave height and period) for engineering applications. These parameters are commonly estimated based on the spectral moments of the energy density (Rice 1945). Most importantly is the zero-order moment, $m_{0}$, which is used to estimate, for example, the so-called ‘significant wave height’ $H_{s}$ (defined as the mean height of the highest one-third of the waves in the field) through the following formula:

(2.10)$$\begin{eqnarray}H_{s}(\boldsymbol{x},t)=4\sqrt{m_{0}^{\prime }},\end{eqnarray}$$

where $m_{0}^{\prime }=m_{0}/(\unicode[STIX]{x1D70C}g)$. Therefore, in order to estimate $H_{s}$ using (2.10) one is required to calculate the transformation from the spectral representation of the action density to $m_{0}$. Using the conventional spectrum of the action density, $N(\boldsymbol{x},\boldsymbol{k},t)$, $m_{0}$ is easily obtained as

(2.11)$$\begin{eqnarray}m_{0}=\unicode[STIX]{x1D70C}\int \unicode[STIX]{x1D70E}(\boldsymbol{x},\boldsymbol{k})N(\boldsymbol{x},\boldsymbol{k},t)\,\text{d}\boldsymbol{k}.\end{eqnarray}$$

However, if cross-correlation terms are taken into account, equation (2.11) is no longer adequate since the cross terms at $(\boldsymbol{x},\boldsymbol{k})$ should not be multiplied by $\unicode[STIX]{x1D70E}(\boldsymbol{x},\boldsymbol{k})$. In order to multiply each term stored at $(\boldsymbol{x},\boldsymbol{k})$ by the correct factor, one must distinguish between the variance term and the cross-correlation terms. Therefore, for cases where cross-correlation terms (e.g. interference terms) are important, a direct substitution of $W(\boldsymbol{x},\boldsymbol{k},t)$ instead of $N(\boldsymbol{x},\boldsymbol{k},t)$ in (2.11) would be inaccurate. A modified formula to calculate $m_{0}$ based on $W(\boldsymbol{x},\boldsymbol{k},t)$ is given as follows:

(2.12)$$\begin{eqnarray}m_{0}=\unicode[STIX]{x1D70C}\int \int \sqrt{\unicode[STIX]{x1D70E}(\boldsymbol{x},\boldsymbol{ k}+\boldsymbol{ k}^{\prime }/2)}\sqrt{\unicode[STIX]{x1D70E}(\boldsymbol{x},\boldsymbol{ k}-\boldsymbol{ k}^{\prime }/2)}\hat{\unicode[STIX]{x1D6E4}}(\boldsymbol{k}^{\prime },\boldsymbol{k},t)\exp (\text{i}\boldsymbol{k}^{\prime }\boldsymbol{\cdot }\boldsymbol{x})\,\text{d}\boldsymbol{k}^{\prime }\,\text{d}\boldsymbol{k},\end{eqnarray}$$

where

(2.13)$$\begin{eqnarray}\hat{\unicode[STIX]{x1D6E4}}(\boldsymbol{k},\boldsymbol{k}^{\prime },t)=\int W(\boldsymbol{x},\boldsymbol{k},t)\exp (-\text{i}\boldsymbol{k}^{\prime }\boldsymbol{\cdot }\boldsymbol{x})\,\text{d}\boldsymbol{x}.\end{eqnarray}$$

Appendix C details the derivation of (2.12) and also provides a simple example that explains why the cross-correlation terms should be scaled differently.

To conclude, the Wigner distribution, $W$, of the action variable, $\unicode[STIX]{x1D713}$, generalizes the concept of the action density spectrum (i.e. $N$), by including the cross-correlation terms that correspond to wave interferences. Once $W$ is known, local field parameters (e.g. $H_{s}$) can be derived and used for practical applications. The last step of the formulation should therefore be devoted to the derivation of a stochastic model for computing the evolution of $W$.

2.3 Evolution equation for the Wigner distribution

The procedure to derive the evolution equation for $W$ is analogous to the procedure presented in Smit & Janssen (2013) and Smit et al. (2015a), and is briefly presented below. Starting with the governing equation of the action variable, equation (2.6), the evolution equation for the correlation function is derived (see e.g. Papoulis and Pillai (2002)) by noting first that

(2.14)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D6E4}(\boldsymbol{x}_{1},x_{2},t)=\langle \unicode[STIX]{x1D713}^{\ast }(\boldsymbol{x}_{2},t)\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}(\boldsymbol{x}_{1},t)+\unicode[STIX]{x1D713}(\boldsymbol{x}_{1},t)\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}^{\ast }(\boldsymbol{x}_{2},t)\rangle ,\end{eqnarray}$$

then, by substituting the governing equation of $\unicode[STIX]{x1D713}$ into the above equation, and using the variable transformation $\boldsymbol{x}_{1}=\boldsymbol{x}+\boldsymbol{x}^{\prime }/2$ and $\boldsymbol{x}_{2}=\boldsymbol{x}-\boldsymbol{x}^{\prime }/2$, one obtains

(2.15)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D6E4}(\boldsymbol{x},\boldsymbol{x}^{\prime },t)=-\text{i}[\unicode[STIX]{x1D714}(\boldsymbol{x}+\boldsymbol{x}^{\prime }/2,-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{ x}^{\prime }}-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}}/2)-\unicode[STIX]{x1D714}(\boldsymbol{x}-\boldsymbol{x}^{\prime }/2,-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{ x}^{\prime }}+\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}}/2)]\unicode[STIX]{x1D6E4}(\boldsymbol{x},\boldsymbol{x}^{\prime },t).\end{eqnarray}$$

The corresponding evolution equation for the Wigner distribution is derived through the Fourier transformation, equation (2.9), and associating the factor $\boldsymbol{x}^{\prime }$ with $i\unicode[STIX]{x1D735}_{\boldsymbol{k}}$ and the operator $-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}}$ with $\boldsymbol{k}$, as

(2.16)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}W(\boldsymbol{x},\boldsymbol{k},t)=-\text{i}\unicode[STIX]{x1D714}(\boldsymbol{x}+\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{k}}/2,\boldsymbol{k}-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}}/2)W(\boldsymbol{x},\boldsymbol{k},t)+\text{c.c.},\end{eqnarray}$$

where $\text{c.c.}$ stands for complex conjugate. For the purpose of interpreting the operation $\unicode[STIX]{x1D714}$ upon $W$, equation (2.16) is written in the following, more explicit, form (see details in appendix D):

(2.17)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}W(\boldsymbol{x},\boldsymbol{k},t)=-\text{i}\unicode[STIX]{x1D714}(\boldsymbol{x},\boldsymbol{k})\exp [\text{i}\overleftarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{x}}\boldsymbol{\cdot }\overrightarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{k}}/2-\text{i}\overleftarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{k}}\boldsymbol{\cdot }\overrightarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{x}}/2]W(\boldsymbol{x},\boldsymbol{k},t)+\text{c.c.},\end{eqnarray}$$

where the arrows indicate the function on which the differential operator should operate, i.e. $\unicode[STIX]{x1D714}$ or $W$.

Formally, equation (2.17) defines the evolution of $W$. Smit & Janssen (2013) showed that essentially two parameters, $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D707}$, govern the order of approximation introduced by a truncated version of the exponential operator in (2.17). The parameter $\unicode[STIX]{x1D6FD}$ arises due to the operation of the first term in the exponential operator (i.e. $\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{k}}/2$), and it represents the ratio between the correlation length scale $L_{c}$ and the medium variation scale $L_{m}$, namely $\unicode[STIX]{x1D6FD}=L_{c}/L_{m}$. The parameter $\unicode[STIX]{x1D707}$ arises due to the operation of the second term in the exponential operator (i.e. $\unicode[STIX]{x1D735}_{\boldsymbol{k}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}/2$) and it is equal to the ratio between the wavelength $L$ that corresponds to $\boldsymbol{k}$ and the characteristic length scale of the interference structures stored in $\boldsymbol{k}$, $L_{W}$, i.e. $\unicode[STIX]{x1D707}=L/L_{W}$. Accordingly, Taylor expansion may applied to define the operator in (2.16) by requiring that both $\unicode[STIX]{x1D6FD}\ll 1$ and $\unicode[STIX]{x1D707}\ll 1$. Under these conditions, the general evolution equation, equation (2.16), can be approximated to $O(\unicode[STIX]{x1D6FD},\unicode[STIX]{x1D707})$ by

(2.18)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}W+\unicode[STIX]{x1D735}_{\boldsymbol{k}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}W-\unicode[STIX]{x1D735}_{\boldsymbol{x}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{k}}W=0,\end{eqnarray}$$

which is exactly the transport equation employed in most commonly used third-generation spectral wave models (e.g. SWAN). Therefore, the conventional transport equation, equation (2.18), is only valid for certain sea conditions for which $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D707}$ are small.

Assuming that the incident wave field is statistically homogeneous, Smit & Janssen (2013) demonstrated that generated cross-correlations (and, therefore, wave interferences) may have an important contribution for cases where the variation scale of the medium is of the same order or smaller than the scale of the correlation length, namely for cases in which $\unicode[STIX]{x1D6FD}\geqslant O(1)$. Obviously, for such cases, the interpretation of the operator in (2.16) using a Taylor expansion is no longer valid. Alternatively, the operator can be partially defined using a Fourier integral (Smit & Janssen 2013), leading to an integro-differential form, which remains valid also for cases in which $\unicode[STIX]{x1D6FD}\geqslant O(1)$, but retains the assumption of weak spatial variability of the statistics of the field ($\unicode[STIX]{x1D707}\ll 1$). This form of the operator is defined as (see appendix D for details)

(2.19)$$\begin{eqnarray}\unicode[STIX]{x1D714}(\boldsymbol{x}+\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{k}}/2,k-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}}/2)W(\boldsymbol{x},\boldsymbol{k},t)=\int \hat{\unicode[STIX]{x1D714}}(\boldsymbol{q},\boldsymbol{k},\boldsymbol{x})(1-\text{i}\overleftarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{k}}\boldsymbol{\cdot }\overrightarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{x}}/2)W(\boldsymbol{x},\boldsymbol{k}-\boldsymbol{q}/2,t)\,\text{d}\boldsymbol{q},\end{eqnarray}$$

where $\hat{\unicode[STIX]{x1D714}}(\boldsymbol{q},\boldsymbol{k},\boldsymbol{x})$ is the Fourier transform of the dispersion relation around the point $\boldsymbol{x}$. Additionally, the part of the operator that results in the common transport terms of (2.18) can be extracted out of the integral in (2.19) (see Smit et al. 2015a). However, for cases in which $\unicode[STIX]{x1D6FD}\geqslant O(1)$, it will be convenient to extract only the spatial transport term ($\unicode[STIX]{x1D735}_{\boldsymbol{k}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}W$) and to leave the refraction term ($\unicode[STIX]{x1D735}_{\boldsymbol{x}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{k}}W$) inside the integral. This is because such cases involve relatively rapid variations in the medium and also narrow spectrum, and therefore require not only high resolution in the spatial space, but also high resolution in the spectral space. Leaving the refraction term inside the integral eliminates the need to evaluate the derivative of $W$ with respect to $\boldsymbol{k}$, and thus prevents excessive resolution in the spectral space. As a consequence, the integral of (2.19) can be computed much more efficiently. To this end, the local value of the dispersion relation at the point $\boldsymbol{x}$ is subtracted from the original dispersion relation and the remainder is defined as $\unicode[STIX]{x0394}\unicode[STIX]{x1D714}(\boldsymbol{x}+\overline{\boldsymbol{x}},\boldsymbol{k})=\unicode[STIX]{x1D714}(\boldsymbol{x}+\overline{\boldsymbol{x}},\boldsymbol{k})-\unicode[STIX]{x1D714}(\boldsymbol{x},\boldsymbol{k})$ (where, due to computational considerations, $\overline{\boldsymbol{x}}$ is defined as $\overline{\boldsymbol{x}}=\boldsymbol{x}^{\prime }/2$; see details in appendix E). With this decomposition, the evolution equation can be rewritten as

(2.20)$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}W+\unicode[STIX]{x1D735}_{\boldsymbol{k}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}W=S_{QC},\end{eqnarray}$$

where $S_{QC}$ is a scattering source term that takes into account the statistical effects of wave refraction and interference induced by medium variations. The expression that defines this source term is given by

(2.21)$$\begin{eqnarray}\displaystyle S_{QC} & = & \displaystyle -\text{i}\int \unicode[STIX]{x0394}\hat{\unicode[STIX]{x1D714}}(\boldsymbol{q},\boldsymbol{k},\boldsymbol{x})(1-\text{i}\overleftarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{k}}\boldsymbol{\cdot }\overrightarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{x}}/2)W(\boldsymbol{x},\boldsymbol{k}-\boldsymbol{q}/2,t)\,\text{d}\boldsymbol{q}\nonumber\\ \displaystyle & & \displaystyle +\,\text{i}\int \unicode[STIX]{x0394}\hat{\unicode[STIX]{x1D714}}(\boldsymbol{q},\boldsymbol{k},\boldsymbol{x})(1+\text{i}\overleftarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{k}}\boldsymbol{\cdot }\overrightarrow{\unicode[STIX]{x1D735}}_{\boldsymbol{x}}/2)W(\boldsymbol{x},\boldsymbol{k}+\boldsymbol{q}/2,t)\,\text{d}\boldsymbol{q}.\end{eqnarray}$$

Note that the subscript $QC$, which indicates this source term, stands for ‘quasi-coherent’ approximation (Smit & Janssen 2013). The notion ‘quasi-coherent’ refers to the assumption that $\unicode[STIX]{x1D707}\ll 1$. Assuming that $\unicode[STIX]{x1D707}$ is small, the model can accurately resolve only the interference patterns with spatial variation, $L_{W}$, larger than the length of the considered wave, $L$.

The transport equation of $W$, equation (2.20), provides a generalization of the conventional transport model (2.18), by allowing statistical interferences to be generated due to the interaction of the wave field with variable bathymetry and currents. In that sense, equation (2.20) can be seen as a generalized action balance equation. In the following, the numerical implementation of (2.20) is discussed.

2.4 Numerical implementation

The numerical implementation of (2.20) is confined to steady-state solutions, for which spatial and spectral discretizations are required. A detailed explanation of the discretization process and how $S_{QC}$ is implemented numerically is given in appendix E. The discretization process results in a coupled system of algebraic equations that is characterized by a matrix of size $N_{x1}N_{x2}N_{k1}N_{k2}\times N_{x1}N_{x2}N_{k1}N_{k2}$, where $N_{j}$ is the number of grid points in the direction $j$. As a consequence of the implicit approach adopted here, where the spatial derivatives and the terms that construct $S_{QC}$ are evaluated at the same spatial point, the coupled system of algebraic equations must be solved iteratively. This is performed using the Gauss–Seidel method, where the rows of the matrix are arranged in accordance with the sweeping approach as detailed in Zijlema & van der Westhuysen (2005). Once a steady-state solution of $W$ is reached, the evaluation of $m_{0}$ which is required for the estimation of certain statistical field parameters is computed through (C 12) (see appendix C for details). The next subsection describes the numerical simulations which are considered in this study.

Figure 1. Wave rays due to $\boldsymbol{k}_{0}$ over a jet-like current field indicated by the solid lines. The rays at $x_{1}=0$ are obliquely incident with an angle of $15^{\circ }$. In addition, the ambient current is marked by arrows. Finally, the dashed vertical lines are sections along which the results of the significant wave height will be displayed.

2.5 Set-up and overview of the considered numerical simulations

Two classical examples of wave–current interactions are considered. The first concerns the evolution of an incoming wave field over a jet-like current. This example is used to verify the model in § 3. In the second example, the field interacts with a vortex-ring current. This example is used in § 4 to study the statistical condition for the effect of wave interferences to appear. A visual description of the spatial variation of the considered current fields is presented by the arrows in figure 1 for the jet-like current and in figure 4 for the vortex ring. Mathematically, these current fields are defined as follows. The jet is defined as

(2.22)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{U}(x_{1},x_{2})=[U_{x_{1}},0],\\ U_{x_{1}}=C_{1}f[\tanh [(x_{2}+R)/(C_{2}R)]-\tanh [(x_{2}-R)/(C_{2}R)]],\\ f=1+\tanh [(x_{1}-R)/(C_{2}R)],\end{array}\right\}\end{eqnarray}$$

where $R=200~\text{m}$, $C_{1}=-0.1~\text{m}~\text{s}^{-1}$ and $C_{2}=0.5$. In this case, the maximum opposing current value is $|U_{x_{1}}\!|_{max}=0.38~\text{m}~\text{s}^{-1}$. Using cylindrical coordinates, the definition of the vortex ring is given by (Mapp et al. 1985)

(2.23)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{U}(r,\unicode[STIX]{x1D703})=[0,U_{\unicode[STIX]{x1D703}}],\\ U_{\unicode[STIX]{x1D703}}=\left\{\begin{array}{@{}l@{}}C_{1}(r/R_{1})^{2},~~r\leqslant R_{1},\\ C_{2}\exp [-(R_{2}-r)^{2}/R_{3}^{2}],\,r\geqslant R_{1},\\ \end{array}\right.\end{array}\!\!\right\}\end{eqnarray}$$

for which the values of $R_{1}$, $R_{2}$, $R_{3}$, $C_{1}$ and $C_{2}$ were chosen identical to those detailed in Belibassakis et al. (2011). In this case, the maximum current value is $|U_{\unicode[STIX]{x1D703}}|_{max}=1.00~\text{m}~\text{s}^{-1}$.

Both of these examples are formulated over a spatial domain of $4000~\text{m}\times 4000~\text{m}$ and a constant depth of $h=10~\text{m}$. Waves enter the domain along the left-hand boundary, on $x_{1}=0$. This is simulated by prescribing an incoming energy density, $E_{0}=E(x_{1}=0,x_{2},k_{1},k_{2})$. Note that, as the incoming wave field is assumed to be statistically homogeneous, the corresponding boundary condition of the Wigner distribution is readily obtained as follows: $W_{0}=E_{0}/\unicode[STIX]{x1D70E}$. Finally, note that the lateral boundaries are treated as periodic.

Table 1. An overview of the considered simulations in terms of their physical, statistical and numerical parameters.

An overview of the simulations considered in this study is given in table 1. These simulations differ by the current type (indicated by the name of the simulation in the first column of the table) and by the parameters characterizing the incoming spectrum, $E_{0}$. In all the simulations the incoming spectrum, $E_{0}$, is defined as a two-dimensional Gaussian centred around $\boldsymbol{k}_{0}$. The incoming spectrum is therefore defined completely by the significant wave height ${H_{s}}_{0}$, the carrier wave period and direction $T_{0}$ and $\unicode[STIX]{x1D703}_{0}$ (which provide the centre point $\boldsymbol{k}_{0}$ through the linear dispersion relation) and the standard deviation $S_{d}^{(k)}$, which are given in the second, third, fourth and fifth column of table 1, respectively. In order to give a more intuitive physical interpretation of the width of the spectrum, the table also provides the corresponding standard deviations of the transformed spectrum written in terms of frequency and direction, $S_{d}^{(f)}$ and $S_{d}^{(\unicode[STIX]{x1D703})}$, given in the sixth and seventh columns. Numerically, $E_{0}$ is represented over the grid $N_{k}$ with a resolution that is determined by $S_{d}^{(k)}$ and the resolution parameter $\unicode[STIX]{x1D6FC}$ (see appendix E) given in the eighth column. The value of $\unicode[STIX]{x0394}x$, on the other hand, cannot be deduced from the table; appendix E guides how to choose a reasonable value for $\unicode[STIX]{x0394}x$. This value is fixed to $\unicode[STIX]{x0394}x=25$  m for all the simulations. In addition, the table also provides the correlation length, $L_{c}$, and the statistical parameter, $\unicode[STIX]{x1D6FD}$, in the ninth and tenth columns. As outlined in appendix E, $L_{c}$ is evaluated using $S_{d}^{(k)}$. This can also be expected by the scaling property of the Fourier transform ($O(L_{c})=O(2/S_{d}^{(k)})$). Throughout the analysis of the results in the following subsections, the value of $L_{c}$ (as opposed to the value taken into account in the numerical model; see appendix E) is defined as $L_{c}=4/S_{d}^{(k)}$. For the considered Gaussian initial distribution, this value equals the so-called $1/\text{e}^{2}$ width that provides the diameter connecting the two points with $1/\text{e}^{2}$ times the maximum value of the correlation function. Finally, the order of magnitude of $\unicode[STIX]{x1D6FD}$ is obtained following its definition, $\unicode[STIX]{x1D6FD}=L_{c}/L_{m}$.

3 Model verification

The main aim of this section is to verify the performance of the QCM. The model is verified through a comparison with REF/DIF 1 (Kirby and Dalrymple 1986), which solves a parabolic approximation of the well-known mild-slope equation (e.g. Dingemans 1997). Since REF/DIF 1 allows for monochromatic, unidirectional forcing at the incident boundary, statistics for multi-directional and irregular incident waves are constructed by superposition of variances, under the assumption that waves at the incident boundary are statistically independent (see details in Chawla et al. (1998)). Additionally, to demonstrate the statistical contribution of the interference terms, the results of the QCM are also compared to the results of the SWAN model (Booij et al. 1999). To this end, the first two simulations detailed in table 1, namely $\text{Jet}_{1}$ and $\text{Jet}_{2}$, are considered.

The simulations $\text{Jet}_{1}$ and $\text{Jet}_{2}$ describe the evolution of waves over the jet-like current field. Ray tracing results (figure 1) show that for this jet-like current the waves refract and form a focal zone close to $x_{1}=2000~\text{m}$, beyond which interference structures may emerge.

The physical pattern described by the rays in figure 1 is also reflected statistically in the results of figures 2 and 3. While the results of the QCM and of REF/DIF 1 agree well and share a similar evolution pattern before and after the crossing zone in both of the simulations, the SWAN results increasingly deviate beyond the crossing zone, where interference effects emerge (see figures 2a–c and 3a,b) (note that the small differences that arise at the lateral boundaries, as for instance that appear in the results of figure 3(b) are due to different boundary conditions assumed in each of the models). The results also show that interference effects are not confined to caustic regions (where geometric optics break down), but rather spread over much greater distances in the down-wave direction, beyond the crossing zone (e.g. figure 2a–c). The differences between the models are less pronounced in the results of simulation $\text{Jet}_{2}$ which is initiated using a broader spectrum (see figures 2d–f and 3c,d). In this case, all three models qualitatively predict a similar spatial structure of $H_{s}$ throughout the domain.

Figure 2. Comparison between QCM, REF/DIF 1 and SWAN in terms of the spatial distribution of the significant wave height. (a–c) The results of the simulation $\text{Jet}_{1}$;(d–f) the results due to $\text{Jet}_{2}$.

Figure 3. Comparison between QCM, REF/DIF 1 and SWAN in terms of the significant wave height along the sections that are indicated in figure 1. (a,b) The results of the simulation $\text{Jet}_{1}$; (c,d) the results due to $\text{Jet}_{2}$.

Model differences are principally due to the statistical contribution of wave interferences. The transport equation employed by third-generation spectral models (e.g. SWAN), equation (2.18), disregards the contribution of cross-correlations (correlations of different wave components), which contain the information about wave interference. The QCM, on the other hand, does account for this information, and therefore, as the statistical contribution of wave interference becomes significant, the discrepancies between the results of the QCM (or REF/DIF 1) and SWAN are more pronounced. Therefore, it is necessary to understand under which conditions the effect of wave interferences is important.

Generally speaking, the importance of the interference effects reduces as the spectrum of the incoming field becomes wider (e.g. Vincent & Briggs 1989). Effectively, the multiple out-of-phase interference patterns generated by each wave component of the incoming field cancel each other out. Consequently, the superposition of the interference patterns becomes smoother as the incoming spectrum becomes wider. This is the reason why differences between the QCM (and REF/DIF 1) and SWAN are larger for $\text{Jet}_{1}$ than $\text{Jet}_{2}$. Whether or not interference effects can be expected may formally be related to the ratio $\unicode[STIX]{x1D6FD}$ between the correlation length scale of the incident wave field, $L_{c}$, and a typical length scale of the medium, $L_{m}$. Interference effects may become significant when $\unicode[STIX]{x1D6FD}\geqslant O(1)$ and are more pronounced for larger values of $\unicode[STIX]{x1D6FD}$ (hence the difference between $\text{Jet}_{1}$ and $\text{Jet}_{2}$). This statistical condition is discussed in detail next.

4 Discussion

The statistical contribution of wave interferences as a function of the parameter $\unicode[STIX]{x1D6FD}$ can be analysed conceptually as follows. Consider a certain point in space beyond the crossing zone, where interference effects are expected to play a role and assuming that the incoming field is monochromatic, for which $\unicode[STIX]{x1D6FD}\rightarrow \infty$. In this case, the correlation function at the considered point will extend over a very large spatial domain ($L_{c}\rightarrow \infty$), and will generally be composed of in-phase variance terms of the scattered field and out-of-phase cross-correlation terms between each pair of scattered waves. The cross-correlation terms include contributions that were generated due to correlation between the incoming field and the interference structures it forms. As the spectrum of the incoming field becomes wider (namely $S_{d}^{(k)}$ becomes larger), the correlation function will extend over smaller domains and, accordingly, $\unicode[STIX]{x1D6FD}$ will take smaller values. The corresponding change in the interference effect can be analysed from the physical point of view, by examining the correlation function, $\unicode[STIX]{x1D6E4}$, or from the spectral point of view, by considering the Wigner distribution, $W$. From the physical point of view, when the incoming spectrum becomes wider and $\unicode[STIX]{x1D6FD}$ reduces, the correlation value between the incoming field and the interference pattern it forms will become smaller, and consequently the contribution of wave interference, at the considered point, reduces as well. In the limit when $\unicode[STIX]{x1D6FD}\rightarrow 0$ (and therefore $L_{c}\rightarrow 0$), this correlation value converges to zero and the contribution of wave interference is eliminated.

Figure 4. Wave rays due to $\boldsymbol{k}_{0}$ over a vortex ring. The rays are indicated by the solid lines, and the ambient current is marked by arrows. Note that in this case, the rays at $x_{1}=0$ are normally incident. Additionally, the dashed vertical lines are sections along which the results of the significant wave height will be displayed. Finally, the dotted lines distinguish between different regions of the wave field.

The spectral point of view examines the representation of the cross-correlation terms in the Wigner distribution. Since the phases of the cross-correlation terms are not necessarily zero, their amplitudes may either be positive or negative, and therefore tend to cancel each other and lose intensity. As a result, when the Wigner distribution at the considered point is integrated over the spectral space for the purpose of computing the total variance, and the corresponding value of, say $H_{s}$, the contribution of the cross-correlation terms will be less pronounced with the increasing of $S_{d}^{(k)}$, and therefore less pronounced with the decreasing of $\unicode[STIX]{x1D6FD}$.

Figure 5. The distribution of the significant wave height due to the interaction between waves and a vortex ring. The panels (a,b) present the results of $\text{Ring}_{1}$ and the panels (c,d) present the results due to $\text{Ring}_{3}$. Additionally, the solid lines represent the wave rays due to $\boldsymbol{k}_{0}$. Finally, the three black points denoted by $P_{1}$, $P_{2}$ and $P_{3}$ indicate the spatial path along which the evolution of the correlation function and the Wigner distribution is considered. Point $P_{1}$ is located at $(1000~\text{m},-525~\text{m})$, $P_{2}$ at $(2000~\text{m},-625~\text{m})$ and $P_{3}$ at $(3000~\text{m},-725~\text{m})$.

4.1 The evolution of the cross-correlation terms

This subsection provides a numerical demonstration of the above discussion of the statistical condition to the appearance of interference effects in the scattered field. The interaction problem between waves and a vortex ring is a convenient example for this purpose. This is due to the fact that in this case, the domain essentially consists of two homogeneous regions separated by a scattering region, which are referred to as the ‘incoming field’, the ‘scatterer’ and the ‘scattered field’, respectively (see figure 4). Consequently, the statistical condition to wave interferences, which says that correlation should emerge between the incoming field and the interference structure it forms, is readily demonstrated through this interaction problem, as it can be replaced by the condition that the correlation function should extend over a larger domain than the effective domain of the vortex ring.

The statistical condition to the appearance of interference effects is examined by considering the evolution of the correlation function and the Wigner distribution for simulations $\text{Ring}_{1}$ and $\text{Ring}_{3}$ (which differ in their initial spectrum width; see table 1) over a specific spatial path. The spatial path was selected such that it would pass over an area where $H_{s}$ is significantly affected by wave interferences (see figure 5). Finally, the contribution of the interference terms is emphasized by comparing the results of the QCM to the corresponding results of SWAN.

In order to identify wave interference effects, the manner in which the cross-correlation terms (which represent the contribution of wave interferences) are represented is explained first (refer also to the definitions in (2.8) and (2.9)). Given two correlated wave components, their contribution in the correlation function results in two variance terms with wavenumbers $\boldsymbol{k}_{1}$ and $\boldsymbol{k}_{2}$, and cross-correlation term (or interference term) with a wavenumber $(\boldsymbol{k}_{1}+\boldsymbol{k}_{2})/2$. The amplitude of the cross-correlation term depends on the amplitudes of the two wave components and their phase difference. If the point around which the correlation function is considered is located at the trough of the interference pattern generated by the two waves, then the amplitude of the cross-correlation term will be negative and vice versa. Also recall that the correlation function presented here follows the definition in (2.7). Consequently, $\unicode[STIX]{x1D6E4}(\boldsymbol{x},\boldsymbol{x}^{\prime })$ is the correlation between $\unicode[STIX]{x1D713}(\boldsymbol{x}+\boldsymbol{x}^{\prime }/2)$ and $\unicode[STIX]{x1D713}^{\ast }(\boldsymbol{x}-\boldsymbol{x}^{\prime }/2)$, which is different from the function that defines the correlation between $\unicode[STIX]{x1D713}(\boldsymbol{x})$ and $\unicode[STIX]{x1D713}^{\ast }(\boldsymbol{x}+\boldsymbol{x}^{\prime })$.

Figure 6. The evolution of the correlation function (a–c) and the corresponding Wigner distribution (d–f) as presented by the spatial points $P_{1}$, $P_{2}$ and $P_{3}$. The values of the results are normalized by $|\unicode[STIX]{x1D6E4}(P_{j},\boldsymbol{x}^{\prime })|_{max}$ and $|W(P_{j},\boldsymbol{k})|_{max}$. These results were obtained for the simulation $\text{Ring}_{1}$ using the QCM.

The analysis starts by examining the evolution results of the correlation function and the Wigner distribution for the simulation $\text{Ring}_{1}$. In this case $\unicode[STIX]{x1D6FD}>1$, and therefore the effect of the cross-correlation terms on the structure of the correlation function and the Wigner distribution is likely to be significant. This is indeed evident by comparing the results of the QCM (figure 6) to the results of SWAN (figure 7). Notable differences clearly appear in the results around $P_{2}$ and $P_{3}$, in the ‘scattered field’, where interference effects are significant. Both of these points are located along the trough of the interference pattern (see figure 5). Indeed, the amplitudes of the cross-correlation terms, which are obtained at these points, are negative, as indicated by the blue areas in the Wigner distribution due to the QCM. Note that these blue areas are located exactly between the red areas which relate to the amplitudes of the variance terms. As expected, the blue areas do not appear in the action density spectrum due to SWAN, as it disregards the cross-correlation terms and only accounts for the variance terms. Moreover, in contrast to SWAN’s results which only accounts for variance terms that are crossing close to the considered points, the QCM also includes contribution of variance terms and related cross-correlation terms that are crossing at some distance away from the considered points. This can be seen by comparing the Wigner distribution due to the QCM and the action density spectrum due to SWAN and by referring to the wave rays in figure 4. Finally, note that the variance areas in the Wigner distribution are somewhat more spread than the corresponding variance areas appearing in the action density spectrum.

The negative values of the cross-correlation amplitudes in the results due to the QCM lead to the fact that the correlation function at these points does not provide the maximum correlation value at its centre (i.e. at $\boldsymbol{x}^{\prime }=0$). Conversely, since the SWAN model ignores the cross-correlation terms, the correlation function will always obtain the maximum value at $\boldsymbol{x}^{\prime }=0$. Therefore, the correlation function as defined in (2.7) does not necessarily show the maximum value at $\boldsymbol{x}^{\prime }=0$ for inhomogeneous fields.

Besides changing the correlation value at the central point, it is difficult to identify the cross-correlation terms directly through the correlation function. However, it is clear that these terms significantly change the structure of the correlation function, as reflected by the differences in the results due to the QCM and SWAN (compare figures 6a–c and 7a–c).

Figure 7. The evolution of the correlation function (a–c) and the corresponding action density spectrum (d–f) as presented by the spatial points $P_{1}$, $P_{2}$ and $P_{3}$. The values of the results are normalized by $|\unicode[STIX]{x1D6E4}(P_{j},\boldsymbol{x}^{\prime })|_{max}$ and $|N(P_{j},\boldsymbol{k})|_{max}$. These results were obtained for the simulation $\text{Ring}_{1}$ using SWAN.

The significant contribution of wave interferences appearing in the results around $P_{2}$ and $P_{3}$ implies that correlation emerges between the ‘incoming field’ and the ‘scattered field’. This is indeed shown by the results of the correlation function around $P_{1}$ in figure 6 (or in figure 7). The results show that the correlation function extends over a much larger domain than the effective domain of the vortex ring and that strong correlation values emerge between the incoming and the scattered field. Accordingly, the generated cross-correlation terms at $P_{1}$ have a clear signature on the structure of the correlation function and the Wigner distribution due to the QCM (compare the results of $P_{1}$ in figures 6 and 7). These cross-correlation terms are transported along with the variance terms, altering dramatically the statistics of the scattered field, as shown by the significant differences between the results of the QCM and SWAN around $P_{2}$ and $P_{3}$.

Figure 8. The evolution of the correlation function (a–c) and the corresponding Wigner distribution (d–f) as presented by the spatial points $P_{1}$, $P_{2}$ and $P_{3}$. The values of the results are normalized by $|\unicode[STIX]{x1D6E4}(P_{j},\boldsymbol{x}^{\prime })|_{max}$ and $|W(P_{j},\boldsymbol{k})|_{max}$. These results were obtained for the simulation $\text{Ring}_{3}$ using the QCM. Note that the scale over which the correlation function is plotted is much smaller than the corresponding scale used to present the results for $\text{Ring}_{1}$.

The differences in the results between the QCM and SWAN for the simulation $\text{Ring}_{3}$ are much less prominent (see figures 8 and 9). The reason for this is that at $P_{1}$, the correlation function extends over a domain with about the same diameter as that of the vortex ring, and only small correlation value arises between the incoming field and the interference structure it forms in the vicinity of the crossing point at $(x_{1},x_{2})=(1365~\text{m},-355~\text{m})$ (see figure 5). As a consequence, at $P_{1}$, the amplitudes of the generated cross-correlation terms are quite low, as shown by the blue area in the Wigner distribution due to the QCM in figure 8. Over the ‘scattered field’ region, at $P_{2}$ and $P_{3}$, the influence of the cross-correlation terms is hardly detected through the correlation function, and, indeed, at these points the correlation functions due to the QCM and SWAN are almost identical. However, the presence of the cross-correlation terms is visible in the Wigner distribution due to the QCM by the blue area located between the variance areas. These cross-correlation terms eventually result in a limited contribution to the statistics of the scattered field, as for instance appears by the spatial distribution of $H_{s}$ in figure 5.

Figure 9. The evolution of the correlation function (a–c) and the corresponding action density spectrum (d–f) as presented by the spatial points $P_{1}$, $P_{2}$ and $P_{3}$. The values of the results are normalized by $|\unicode[STIX]{x1D6E4}(P_{j},\boldsymbol{x}^{\prime })|_{max}$ and $|N(P_{j},\boldsymbol{k})|_{max}$. These results were obtained for the simulation $\text{Ring}_{3}$ using SWAN. Note that the scale over which the correlation function is plotted is much smaller than the corresponding scale used to present the results for $\text{Ring}_{1}$.

To conclude, the examination of the evolution of the correlation function and the Wigner distribution verifies the statistical condition for the generation of cross-correlation as was introduced conceptually in the beginning of this section. Moreover, the examination also demonstrates numerically that the correlation value between the incoming field and the interference structure it forms determines the dominance of the interference patterns in the scattered field.

4.2 The validity of the QCM

The final issue that is discussed here is the validity of the QCM versus the validity of SWAN over the parameter $\unicode[STIX]{x1D6FD}$. As was explained in the derivation of the QCM in § 2.3 and following the presentation of the results so far, the QCM, in contrast to SWAN, seems to remain statistically valid for $\unicode[STIX]{x1D6FD}\geqslant O(1)$. The reason for this was extensively discussed in the previous subsection, and in short is simply because the QCM accounts for statistical inhomogeneity of the wave field, generated due to interference effects.

Figure 10. On the validity of the QCM (a,b) and SWAN (c,d) over the parameter $\unicode[STIX]{x1D6FD}$, shown through the convergence of the significant wave height to the result of REF/DIF 1 with $\text{Ring}_{1}$. The results are given along Sections A and B that are indicated in figure 4.

Figure 11. On the validity of the QCM (a,b) and SWAN (c,d) over the parameter $\unicode[STIX]{x1D6FD}$, shown through the convergence of the significant wave height to the result of REF/DIF 1 with $\text{Ring}_{1}$. The results are given along Sections C and D that are indicated in figure 4.

The validity of the QCM over $\unicode[STIX]{x1D6FD}$ is presented by demonstrating the convergence of its results, obtained with an increasing value of $\unicode[STIX]{x1D6FD}$, to a single result of REF/DIF 1 obtained with a specific high value of $\unicode[STIX]{x1D6FD}$. To this end, the QCM is used to compute $H_{s}$ along the sections shown in figure 4 using simulations $\text{Ring}_{1}$, $\text{Ring}_{2}$ and $\text{Ring}_{3}$ which are defined with a decreasing value of $\unicode[STIX]{x1D6FD}$ (i.e. $\text{Ring}_{1}$ is defined with the highest $\unicode[STIX]{x1D6FD}$ value, whereas $\text{Ring}_{3}$ is defined with the lowest $\unicode[STIX]{x1D6FD}$ value; see also table 1). In addition, the result due to REF/DIF 1 is obtained through $\text{Ring}_{1}$. Finally, the convergence of the QCM results to the result of REF/DIF 1 is shown in figures 10(a,b) and 11(a,b). The same procedure is performed using SWAN and is presented in figures 10(c,d) and11(c,d).

Over the ‘scatterer’ region, before the focusing zones, SWAN seems to remain valid (see figure 10, section A) even for the highest $\unicode[STIX]{x1D6FD}$ considered, which corresponds to the simulation $\text{Ring}_{1}$. However, over the ‘scattered field’ region, where interference effects emerge, SWAN does not converge to REF/DIF 1 when $\unicode[STIX]{x1D6FD}$ increases. On the other hand, the QCM does converge to REF/DIF 1, and seems to remain valid for the scattered field as well.

It is important to remember that the capabilities of the QCM over $\unicode[STIX]{x1D6FD}$ involve a constraint. This constraint is that $\unicode[STIX]{x1D716}\ll 1$, introduced by the deterministic model, equation (2.6), which underlies the development of the QCM. Finally, recall that the QCM is also limited to small values of $\unicode[STIX]{x1D707}$, which basically limits its capabilities to accurately evolve interference terms with a wavelength of $L_{W}\leqslant O(L)$, where $L$ is the wavelength of the considered point $\boldsymbol{k}$ (see details in Smit & Janssen (2013)).

5 Conclusions

This study presents the development of a statistical model for problems of wave–current interaction, taking into account the effect of wave interferences. The theoretical basis of this model lies in the definition of the Wigner distribution, $W$, and of the action variable, $\unicode[STIX]{x1D713}$. This distribution provides a complete spectral description of the second-order statistics of the wave field. It includes cross-correlation terms, which provide the statistical information about wave interferences. As such, $W$ generalizes the concept of the action density spectrum, $N$, which only accounts for the information of wave variances.

Using the procedure described in Smit & Janssen (2013) and Smit et al. (2015a), an evolution model for $W$ (the QCM) is developed. This model provides a generalization of the conventional action balance model (presently employed by third-generation spectral wave models, e.g. SWAN and WAVEWATCH III), by allowing the generation and transportation of statistical wave interferences.

The effect of wave interferences can contribute significantly for cases where the variation scale of the medium is of the same order or smaller than the scale of the correlation length, namely for cases in which $\unicode[STIX]{x1D6FD}\geqslant O(1)$. This statistical condition is explicitly examined for scenarios where the incoming field is statistically homogeneous, but develops inhomogeneity while propagating over ambient currents. Specifically, in order to obtain a statistical signature of wave interferences, the incident and scattered fields should be correlated, with the dominance of the interference effect determined by the correlation value itself.

In cases where this correlation is strong, the interference patterns alter the statistics of the field significantly. The resulting effect on the significant wave height, $H_{s}$, is demonstrated through two examples of wave–current interaction and by a comparison to the SWAN model. It is demonstrated that in such cases, interference effects dramatically change the distribution of $H_{s}$, not only in the vicinity of wave-focusing areas, but also at a significant distance away from the focusing points.

It is therefore concluded that for regions involving rapid variability in medium (e.g. coastal regions or oceanic regions which tend to contain submesoscale currents), consideration of the statistical information of wave interference might by crucial for many applications, such as wave-induced circulation and transport processes in coastal regions or for prediction of extreme elevations in the open ocean.