Modelling statistical wave interferences over shear currents

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.

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.
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) where N represents the spectrum of the action density, being equal to the spectrum of the energy density, E, divided by the intrinsic frequency, σ . The propagation part, on the left-hand side, describes the kinematic behaviour of the field as it propagates through slowly varying current, U, and bathymetry, with propagation velocities C k and C x over wavenumber space, k = (k 1 , k 2 ), and physical space, 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.

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.

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, φ(x, t) and η(x, t). It is assumed that the medium changes slowly so that the ratio, = L/L m , between the characteristic wavelength, L, and the characteristic length scale of medium variation, L m , is small ( 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 obey to the following general dispersion relation (e.g. Dingemans 1997): (2.1) Variations in the medium are introduced by the ambient current, U(x), and by the water depth, h(x). Using the medium information and the definition of the intrinsic frequency, σ (x, k) = √ |k|g tanh(|k|h), the value of the absolute frequency, ω, is 891 A2-4 G. Akrish, P. Smit, M. Zijlema and A. Reniers obtained through (2.1), where |k| is the magnitude of the local wavenumber, defined as |k| = k 2 1 + 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 socalled action variable (e.g. Besieris & Tappert 1976;Krasitskii 1994), ψ, which is defined as where a(x, −i∇ x ) is a pseudo-differential operator that is associated with the symbol a(x, k) = √ σ (x, k) (see detailed definition of this operator in appendix A). The convenience of working with the action variable, ψ, becomes significant in the formulation of the second-order statistics of the field, since, following its definition, second-order statistical functions of ψ (e.g. the correlation function) are inherently related to the definition of the wave action (Bretherton & Garrett 1968). As a consequence, the action variable, ψ, is intimately related to the mean action density and the mean energy density through the following expressions: where ρ is the water mass density and the angular brackets, · · · , should be read as an ensemble average. The variable m 0 provides a leading-order estimation (in ) of the mean energy density (also known as the zero-order moment of the spectral energy density) and it is defined as follows: where now (in (2.3) and (2.5)) σ (x, −i∇ x ) represents a pseudo-differential operator that is associated with the intrinsic frequency, σ (x, k), and the subscript 0 indicates O(1) terms (refer to appendix A for the definition of σ (x, −i∇ x ) and its leading-order operation, e.g. (σ φ) 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 ψ (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) where ω(x, −i∇ x ) is a pseudo-differential operator that is associated with the dispersion relation, ω(x, k) (see appendix A To summarize, the formulation presented here considers a random, linear and slowly varying wave field, which is concisely represented by the action variable, ψ. 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).

Second-order statistics
Following the statistical assumptions for the surface variables, η and φ, and following the linearity of the definition (2.2), the action variable ψ(x, t) is said to be a zeromean, 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) The statistical information carried by the correlation function is better seen using its spectral form, written as where k and k are defined as the average and difference of two interacting wavenumbers, namely k = (k 1 + k 2 )/2 and k = k 1 − k 2 . In additionΓ (k, k , t) is defined asΓ (k, k , t) = ψ (k + k /2, t)ψ * (k − k /2, t) . The expression above reveals the spectral content of the correlation function. It shows that, in general, Γ oscillates with a wavenumber difference k over the space 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Γ with respect to k and accounts only for a narrow window around k = 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 Γ from x to k only allows for slow variations of the variance terms of the field. This spectrum is the conventional action density spectrum, N(x, k, t). In this study, however, statistical inhomogeneity of the wave field is taken into account by considering the full spectrum provided byΓ with respect to k . In this general case, the corresponding spectral representation of wave action follows the definition of the Wigner distribution, W(x, k, t): (2.9) Therefore, the Wigner distribution of ψ 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 ψ, 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: where m 0 = m 0 /(ρ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(x, k, t), m 0 is easily obtained as m 0 = ρ σ (x, k)N(x, k, t) dk.
(2.11) However, if cross-correlation terms are taken into account, equation (2.11) is no longer adequate since the cross terms at (x, k) should not be multiplied by σ (x, k). In order to multiply each term stored at (x, 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(x, k, t) instead of N(x, k, t) in (2.11) would be inaccurate. A modified formula to calculate m 0 based on W(x, k, t) is given as follows: (2.13) 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, ψ, 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.

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) then, by substituting the governing equation of ψ into the above equation, and using the variable transformation x 1 = x + x /2 and x 2 = x − x /2, one obtains (2.15) The corresponding evolution equation for the Wigner distribution is derived through the Fourier transformation, equation (2.9), and associating the factor x with i∇ k and the operator −i∇ x with k, as where c.c. stands for complex conjugate. For the purpose of interpreting the operation ω upon W, equation (2.16) is written in the following, more explicit, form (see details in appendix D): where the arrows indicate the function on which the differential operator should operate, i.e. ω or W. Formally, equation (2.17) defines the evolution of W. Smit & Janssen (2013) showed that essentially two parameters, β and µ, govern the order of approximation introduced by a truncated version of the exponential operator in (2.17). The parameter β arises due to the operation of the first term in the exponential operator (i.e. ∇ x · ∇ k /2), and it represents the ratio between the correlation length scale L c and the medium variation scale L m , namely β = L c /L m . The parameter µ arises due to the operation of the second term in the exponential operator (i.e. ∇ k · ∇ x /2) and it is equal to the ratio between the wavelength L that corresponds to k and the characteristic length scale of the interference structures stored in k, L W , i.e. µ = L/L W . Accordingly, Taylor expansion may applied to define the operator in (2.16) by requiring that both β 1 and µ 1. Under these conditions, the general evolution equation, equation (2.16), can be approximated to O(β, µ) by which is exactly the transport equation employed in most commonly used thirdgeneration spectral wave models (e.g. SWAN). Therefore, the conventional transport equation, equation (2.18), is only valid for certain sea conditions for which β and µ 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 β 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 β O(1), but retains the assumption of weak spatial variability of the statistics of the field (µ 1). This form of the operator is defined as (see appendix D for details) whereω(q, k, x) is the Fourier transform of the dispersion relation around the point 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 β O(1), it will be convenient to extract only the spatial transport term (∇ k ω · ∇ x W) and to leave the refraction term (∇ x ω · ∇ 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 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 x is subtracted from the original dispersion relation and the remainder is defined as With this decomposition, the evolution equation can be rewritten as 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 Note that the subscript QC, which indicates this source term, stands for 'quasicoherent' approximation (Smit & Janssen 2013). The notion 'quasi-coherent' refers to the assumption that µ 1. Assuming that µ 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.

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 × 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 FIGURE 1. Wave rays due to 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 • . 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.
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.
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 where R = 200 m, C 1 = −0.1 m s −1 and C 2 = 0.5. In this case, the maximum opposing current value is |U x 1 | max = 0.38 m s −1 . Using cylindrical coordinates, the definition of the vortex ring is given by (Mapp et al. 1985)  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 θ | max = 1.00 m s −1 .
Both of these examples are formulated over a spatial domain of 4000 m × 4000 m and a constant depth of h = 10 m. Waves enter the domain along the left-hand boundary, on x 1 = 0. This is simulated by prescribing an incoming energy density, . 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 /σ . Finally, note that the lateral boundaries are treated as periodic.
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 k 0 . The incoming spectrum is therefore defined completely by the significant wave height H s0 , the carrier wave period and direction T 0 and θ 0 (which provide the centre point k 0 through the linear dispersion relation) and the standard deviation S (k) d , 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 ( f ) d and S (θ) d , 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 (k) d and the resolution parameter α (see appendix E) given in the eighth column. The value of x, on the other hand, cannot be deduced from the table; appendix E guides how to choose a reasonable value for x. This value is fixed to x = 25 m for all the simulations. In addition, the table also provides the correlation length, L c , and the statistical parameter, β, in the ninth and tenth columns. As outlined in appendix E, L c is evaluated using S (k) d . This can also be expected by the scaling property of the Fourier transform . 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 (k) d . For the considered Gaussian initial distribution, this value equals the so-called 1/e 2 width that provides the diameter connecting the two points with 1/e 2 times the maximum value of the correlation function. Finally, the order of magnitude of β is obtained following its definition, β = L c /L m .

Model verification
The main aim of this section is to verify the performance of the QCM. (d-f ) the results due to Jet 2 .
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 Jet 1 and Jet 2 , are considered. The simulations Jet 1 and 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 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 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.
Model differences are principally due to the statistical contribution of wave interferences. The transport equation employed by third-generation spectral models  figure 1. (a,b) The results of the simulation Jet 1 ; (c,d) the results due to Jet 2 .
(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 Jet 1 than Jet 2 . Whether or not interference effects can be expected may formally be related to the ratio β 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 β O(1) and are more pronounced for larger values of β (hence the difference between Jet 1 and Jet 2 ). This statistical condition is discussed in detail next.

Discussion
The statistical contribution of wave interferences as a function of the parameter β 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 β → ∞. In this case, the x 2 (m) x 1 (m) FIGURE 4. Wave rays due to 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.
correlation function at the considered point will extend over a very large spatial domain (L c → ∞), 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 (k) d becomes larger), the correlation function will extend over smaller domains and, accordingly, β 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, Γ , 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 β 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 β → 0 (and therefore L c → 0), this correlation value converges to zero and the contribution of wave interference is eliminated.
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 (k) d , and therefore less pronounced with the decreasing of β. 2.0 P 1 P 2 P 3 P 1 P 2 P 3 P 1 P 2 P 3 P 1 P 2 P 3 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 Ring 1 and the panels (c,d) present the results due to Ring 3 . Additionally, the solid lines represent the wave rays due to 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 m, −525 m), P 2 at (2000 m, −625 m) and P 3 at (3000 m, −725 m).

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 Ring 1 and 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)  . 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 |Γ (P j , x )| max and |W(P j , k)| max . These results were obtained for the simulation Ring 1 using the QCM. 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 crosscorrelation 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 k 1 and k 2 , and cross-correlation term (or interference term) with a wavenumber (k 1 + 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, Γ (x, x ) is the correlation between ψ(x + x /2) and ψ * (x − x /2), which is different from the function that defines the correlation between ψ(x) and ψ * (x + x ).
The analysis starts by examining the evolution results of the correlation function and the Wigner distribution for the simulation Ring 1 . In this case β > 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   . 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 |Γ (P j , x )| max and |N(P j , k)| max . These results were obtained for the simulation Ring 1 using SWAN.
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 x = 0). Conversely, since the SWAN model ignores the cross-correlation terms, the correlation function will always obtain the maximum value at x = 0. Therefore, the correlation function as defined in (2.7) does not necessarily show the maximum value at x = 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).
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  . 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 |Γ (P j , x )| max and |W(P j , k)| max . These results were obtained for the simulation 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 Ring 1 .
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 .
The differences in the results between the QCM and SWAN for the simulation 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 m, −355 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.
To conclude, the examination of the evolution of the correlation function and the Wigner distribution verifies the statistical condition for the generation of   . 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 |Γ (P j , x )| max and |N(P j , k)| max . These results were obtained for the simulation 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 Ring 1 .
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.

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 β. 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 β 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.
The validity of the QCM over β is presented by demonstrating the convergence of its results, obtained with an increasing value of β, to a single result of REF/DIF 1 obtained with a specific high value of β. To this end, the QCM is used to compute H s along the sections shown in figure 4 using simulations Ring 1 , Ring 2 and Ring 3 which are defined with a decreasing value of β (i.e. Ring 1 is defined with the highest β value, whereas Ring 3 is defined with the lowest β value; see also table 1). In addition, the result due to REF/DIF 1 is obtained through 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 β considered, which corresponds to the It is important to remember that the capabilities of the QCM over β involve a constraint. This constraint is that 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 µ, which basically limits its capabilities to accurately evolve interference terms with a wavelength of L W O(L), where L is the wavelength of the considered point k (see details in Smit & Janssen (2013)).

Conclusions
This study presents the development of a statistical model for problems of wavecurrent 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, ψ. 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 β 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.
which, after Fourier transform with respect to k, reduces to k=k 0 dq dp, (A 5) and eventually leading to the following asymptotic form: In the following, the operation of the Weyl operator, using its asymptotic form, equation (A 6), is examined for two examples: the case of a plane wave over a homogeneous medium and the case of a plane wave propagating over a slowly varying medium. The operation of the Weyl operator in the first case is easily worked out, as in this case the dispersion relation is not a function of x, and also the spatial derivatives on the representative variable of the field can be resolved directly. This is demonstrated as follows. The plane wave is represented by the surface potential as φ = A exp[i(k 0 · x − ω 0 t)/ ], and consequently the operation ωφ is obtained by the following: which is the desired result as detailed at the beginning of appendix B.
In the second example, the considered wave component is represented by the surface potential as φ = A(x) exp[(S − iω 0 t)/ ]. In this case, the operation of the Weyl operator is not immediately seen. The leading-order (O(1)) term is obtained when the first exponent in (A 6) is taken to be equal to one, and the spatial derivative of the second exponent, D x , operates only on the exponent of φ. This term is the dispersion relation, equation (2.1). Terms of O( ) are obtained for three different sets of conditions. Two of these sets instruct one to take the first exponent in (A 6) to be equal to one, and at each expansion order of the second exponent the spatial derivative, D x , should operate once on A for the one set or once on S, as instructed by the other set. The third term is obtained using the second term of the expansion of the first exponent of the operator and when the spatial derivative of the second exponent, D x , operates only on the exponent of φ. This description is summarized as follows: This closes the formal definition of the Weyl operator, equation (A 3), and its asymptotic form, equation (A 6). The operation of the Weyl operator on representative variables in the homogeneous and weakly inhomogeneous cases is used next to verify the evolution equation of the action variable, ψ.

Appendix B. On the evolution equation of the action variable
This appendix aims to demonstrate that the starting point equation (the evolution equation of the action variable), equation (2.6), is exact for the case of linear wave propagation over a homogeneous medium, and reduces to the correct evolution equations for the weakly inhomogeneous case. In the latter, the correct evolution equations are the dispersion relation, equation (2.1), and the well-known action balance equation (e.g. Dingemans 1997). For these purposes, the starting point equation, equation (2.6), is written once again, using the slow scale coordinates, as follows: where, as in the previous appendix, these slow scale coordinates are indicated using the same variable notation of the fast scale coordinates, x, t, in order to avoid cumbersome formulations. As introduced, the first aim here is to show that the starting point equation, equation (B 1), describes exactly the correct solution of an incoming plane wave over a homogeneous medium. To this end, the following example will be worked out. Considering a specific domain of interest, the example assumes that into one of the domain's boundaries enters a monochromatic wave field with an absolute frequency ω 0 . In this case, the surface variables of the considered wave obey to the following form: where B 0 = iA 0 (ω 0 − U · (D x S))/g, as follows from the linear relation between φ and η (e.g. Dingemans 1997): The wavenumber, D x S, is constant, and the constant amplitude, A 0 , assumed to be a 'proper' (e.g. Lapidoth 2017) Gaussian random variable, namely A 2 0 = 0. The corresponding form of ψ is obtained following its definition, equation (2.2): where the amplitudes C 0 and D 0 are defined as C 0 = ga −1 B 0 + iaA 0 and D 0 = ga −1 B * 0 + iaA * 0 , and, following (A 7), a = √ σ (D x S). Using these starting points, it is now aimed to show that (B 1) produces the correct magnitude of D x S. By substituting (B 4) into (B 1), one obtains that the first solution (with the amplitude C 0 ) produces, as required, the equation ω 0 = ω(D x S). Substituting this result into the definition of B 0 , which becomes B 0 = iA 0 σ (D x S)/g, provides the necessary result that D 0 = 0 (this is necessary, since the second term of ψ is not a solution of (B 1)), whereas the value of C 0 is then given by C 0 = 2iA 0 √ σ (D x S). The second aim of this appendix is to show that the starting point equation (B 1) reduces to the correct evolution equations for the weakly inhomogeneous case. To do this, the same example as in the homogeneous case is considered. Following the Wentzel-Kramers-Brillouin method (e.g. Holmes 2012), and assuming that the bathymetry and the current are not time dependent, the surface velocity potential and elevation, φ and η, can be defined to leading order in by (B 2), where now A 0 and S are functions of x. The corresponding leading-order definition of ψ is given by (B 4), but now, following (A 8), a = √ σ (x, D x S). By substituting this definition into (B 1) and using the result of (A 8), the O(1) eikonal equation is obtained to be ω 0 = ω(x, D x S), which can be written as Using this result, B 0 is obtained to be B 0 = iA 0 σ (x, D x S)/g, leading to the same necessary result as before, that D 0 = 0, whereas C 0 is given by C 0 = 2iA 0 √ σ (x, D x S). The remaining unknown, A 0 , is found through the following O( ) transport equation: which alternatively can be written as In order to verify that (B 7) is the well-known action balance equation (e.g. Dingemans 1997), one should check that ρ |ψ| 2 = ρ |C 0 | 2 indeed defines the mean action density. To see this, one may calculate |ψ| 2 directly, using the definition of ψ, equation (2.2), which results in where the subscript 0 indicates O(1) terms. The expression in the square brackets is exactly the mean energy density m 0 that was introduced in (2.5). Clearly, the first term of the expression represents the mean potential energy density, though the connection of the second term to the density of the kinetic energy might not be immediately obvious and should be clarified. Following (A 8), the first-order operation, (σ φ) 0 , is defined as (σ φ) 0 = σ (x, D x S)φ. Substituting this result into (B 8), the second term becomes equal to the mean kinetic energy density for a homogeneous medium that is defined by the local values, h(x) and U(x) (e.g. § 6.3 in Van Groesen and Molenaar (2007)). Therefore, dividing the expression in the square brackets by the local intrinsic frequency, σ , leads to the definition of the mean action density (Bretherton & Garrett 1968).

Appendix C. From Wigner distribution to local energy
To motivate the necessity for an alternative formula that links the energy density, m 0 , and the Wigner distribution, W(x, k, t), the following example is discussed (note that here the formulation returns to be written in terms of the original spatial coordinates, x and t). The example assumes an idealized sea state composed of three coherent and forward-propagating wave components in a one-dimensional and homogeneous medium. The waves are defined with the following wavenumbers: k 1 , k 2 and k 3 = (k 1 + k 2 )/2. Accordingly, at a certain moment in time t 0 , the wave field may be represented as follows: Modelling statistical wave interferences over shear currents 891 A2-25 where A n and B n are complex random amplitudes. The amplitude B n is given by B n = iA n σ n /g, as follows from the linear relation between φ and η (see equation (B 3) in appendix B). The mean energy density, m 0 in this example, is derived using the expression in (2.5), and can be written as follows: The approach employed in this study to get to (C 2) is via the Wigner distribution of the action variable, ψ, which for this example is given by where δ(k) is now serving as the usual delta function. Now it can be seen explicitly that a simple substitution of (C 3) into (2.11) will not provide the result described in (C 2). As an example, consider the components in W that multiply the function δ[k − (k 1 + k 2 )/2]. Following the definition of k 3 , there are two such components, 2σ 3 |A 3 | 2 /g and 2 x] + c.c}/g, which are related to the variance of the third wave component and to the correlation between the first and the second wave components, respectively. If W is substituted into (2.11), both of these components will be factored by σ 3 , which by referring to (C 2) will not lead to the correct interference term between the first and the second wave components. It is concluded that in order to calculate the mean energy density correctly, one must distinguish between the variance term and the interference terms for each k, and only then multiply by the correct factor. The derivation of an alternative formula to obtain m 0 based on W is detailed below.
The starting point of the following derivation is the relation between the action varible, ψ, and the mean energy density, m 0 , as given in (2.4), recalling that a(x, D x ) is a Weyl operator associated with the square root of the intrinsic angular frequency, σ 1/2 (x, k). Following the definition of the Weyl operator, the expression in (2.4) can be written as , t) dq 1 dq 2 dp 1 dp 2 . (C 4) By substituting the Fourier transform of the correlation function, Γ = ψ(x + p 1 , t)ψ * (x − p 2 , t) , the following expression is obtained: × exp(ik 1 · p 1 + ik 2 · p 2 )Γ (k 1 , k 2 , t) exp[ix · (k 1 − k 2 )] dk 1 dk 2 dq 1 dq 2 dp 1 dp 2 . (C 5) Integrating the above expression with respect to q 1 and q 2 leads to |a(x, D x )ψ| 2 = â(x + p 1 /2, p 1 )â * (x − p 2 /2, p 2 ) exp(ik 1 · p 1 + ik 2 · p 2 ) ] dk 1 dk 2 dp 1 dp 2 . (C 6) By assuming that σ is independent of x, equation (C 6) reduces to the following expression: |a(x, D x )ψ| 2 = σ (k + k /2) σ (k − k /2)Γ (k , k, t) exp(ik · x) dk dk, (C 7) where the change of variables k 1 = k + k /2 and k 2 = k − k /2 was applied. Otherwise, equation (C 6) can be integrated once more. Now the integration is performed with respect to p 1 and p 2 , leading to the following result: Ultimately, if terms of O( ) are omitted, and by applying the same change of variables over the wavenumber space (as indicated above), the following approximation for the zero-order moment of the energy density spectrum, m 0 , is derived: m 0 ∼ ρ σ (x, k + k /2) σ (x, k − k /2)Γ (k , k, t) exp(ik · x) dk dk, (C 9) which generalizes (C 7) for the case with a slowly varying medium. The numerical implementation of (C 9) is not straightforward. It requires performing a Fourier transform of the Wigner distribution for each wavenumber, k, and around each spatial location, x. In addition, one should distinguish, at each location x, between the Fourier components that relate to the slow variation of the variances and the Fourier components of the cross-correlation terms.
A different direction to obtain an evaluation of m 0 stems from an alternative formulation of (C 9). This formulation is detailed as follows. At first, the Fourier components,Γ (k , k, t) exp(ik · x), are replaced by the Fourier transform of the Wigner distribution around x: m 0 ∼ ρ f (x, k , k)W(x + x, k, t) exp(−ik · x) dx dk dk, (C 10) where f (x, k , k) = σ (x, k + k /2) σ (x, k − k /2). Then, assuming that the Wigner distribution around x can be expressed as a Taylor expansion, and after integrating over x, the following approximation of m 0 is obtained: where, due to the symmetry of f around k = 0, the exponent in the integral of the last expression can be replaced by a cosine. As opposed to the numerical implementation of (C 9), the implementation of (C 11) is straightforward. The first term in the expansion of (C 11) is exactly the formula used in SWAN, as given by (2.11). The high-order terms provide corrections (of O(µ); see the definition of µ in § 2.3) to the cross-correlation components that are stored in k. Ultimately, in the numerical examples, m 0 is evaluated up to second order using the following expression: A last and important step that is discussed here, which is necessary to the numerical implementation of (2.20), is the representation of the left-hand side of (D 7) in terms of the correlation function, Γ (x, x , t). One way to get to this representation involves a few algebraic steps. The other way is to see it directly through the convolution theorem. In order to present the second way, the multiplication termω(q, k) exp(iq · x) is replaced by the Fourier transform of ω around x, represented byω(q, k, x). Then, using the convolution theorem, the following equation is obtained: Appendix E. On the numerical model The steady-state numerical solution of (2.20) uses the following two-dimensional and equispaced grids: N x , N k , N x , N q (where x = x /2). These grids are constructed using the following spatial and spectral steps: x, k, x, q. The value of k is chosen according to the standard deviation, S (k) d , of the incoming wave spectrum as k = S (k) d /α, where α 1 serves as a resolution factor. Additionally, to ease the computation of the source term S QC , the value of k is selected such that k = q/2 (this selection prevents the need to perform interpolation in the calculation of the integral in (2.21)).
Next, the choice of q is explained. This choice stems from the fact that for any realistic sea state, the correlation function around a certain point, x, will effectively have a compact support in |x | < L c /2, where L c is the correlation length. Consequently, and as implied by (D 8), instead of an integral operation, the source term in (2.20), S QC , can be calculated as a discrete convolution between ω and W (and their derivatives) over the grid N q . This is done without introducing any discretization error if q 4π/L c (for details that also include the additional treatment required to compute the discrete version of ω, see Smit et al. (2015a)). If q is chosen such that q = 4π/L c , then the applied value of L c , which is taken into account in the numerical model, can be found from the definition of k as L c = 2πα/S (k) d , which is consistent with the expected order that should characterize the standard deviation of the envelope of the correlation function (O(1/S (k) d )). The choice of x is argued in a similar way to the selection of q, but now knowledge about the boundaries of W over N k at a certain point x (beyond which W equals zero) is required. As for the correlation function, also here, the boundaries are not easily predicted in advance, because the support of W over N k can change significantly over N x . For accuracy, the necessary x is evaluated in accordance with the boundaries introduced by N k . A more economical selection of x, which also introduces an acceptable error, is described by Smit et al. (2015b). Note that by introducing x, the summation in S QC becomes limited to the region [−q max , q max ], where q max = π/ x. Finally, the derivation of x is considered. Its value should be selected according to the characteristic variation length of W over x, and according to the adopted scheme for treating the spatial derivatives. Here, the second-order upwind scheme (see Hirsch 2007) is used. The local error introduced by this scheme is of O[( xµ/L) 3 ], where L is the characteristic wavelength (L/µ represents the characteristic length of wave interference at the considered location). Therefore, x should be chosen to be small enough so that the global error due to such a magnitude of local error would be acceptable.