## Appendix A. The Weyl operator and its asymptotic form

The purpose of this appendix is to provide the definition of the pseudo-differential operators employed in this study. To this end, the operator $\unicode[STIX]{x1D714}(\boldsymbol{x},-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}})$ will serve as a representative (the following also applies to the operators $\unicode[STIX]{x1D70E}(\boldsymbol{x},-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}})$ and $a(\boldsymbol{x},-\text{i}\unicode[STIX]{x1D735}_{\boldsymbol{x}})$ that were introduced in § 2). Additionally, for convenience, the expressions in this appendix (and in appendix B) are presented using the slow scale coordinates $\boldsymbol{x}_{m}=\unicode[STIX]{x1D716}\boldsymbol{x}$ and $t_{m}=\unicode[STIX]{x1D716}t$. However, in order to avoid cumbersome formulations, the subscript $m$ indicating these slow scale coordinates will be removed, keeping in mind that for the purposes of this appendix (and appendix B), $\boldsymbol{x}$ and $t$ are now serving as the slow scale coordinates. An additional notation is $D_{j}$, which will be used here and in the following appendices to represent the operator $-\text{i}\unicode[STIX]{x1D735}_{j}$.

The definition of the pseudo-differential $\unicode[STIX]{x1D714}(\boldsymbol{x},\unicode[STIX]{x1D716}D_{\boldsymbol{x}})$ is based on its association with a ‘phase-space’ symbol (a function which is defined in $(\boldsymbol{x},\boldsymbol{k})$ space). Here, it is assumed that such a ‘phase-space’ symbol can be defined locally (in this case, it is the usual dispersion relation, equation (2.1)), which basically requires that the characteristic length scale of the medium variation is much larger than the considered wavelength (e.g. Dingemans Reference Dingemans1997), i.e. that $\unicode[STIX]{x1D716}\ll 1$.

Given a ‘phase-space’ symbol, the corresponding operator in the physical space can be defined through the association between $\boldsymbol{k}$ and $D_{\boldsymbol{x}}$. However, because $\boldsymbol{x}$ and $D_{\boldsymbol{x}}$ do not commute, one must follow an association rule for an arbitrary symbol. Here, the Weyl rule of association is adopted (e.g. Cohen Reference Cohen2012), which is defined through the following Fourier transform of $\unicode[STIX]{x1D714}(\boldsymbol{x},\boldsymbol{k})$:

Then, the Weyl operator is obtained by substituting the operator $\unicode[STIX]{x1D716}D_{\boldsymbol{x}}$ instead of $\boldsymbol{k}$, which provides the following expression:

and which can be simplified using the commutator value, $[\text{i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{x},\unicode[STIX]{x1D716}i\boldsymbol{p}\boldsymbol{\cdot }D_{\boldsymbol{x}}]=-\text{i}\unicode[STIX]{x1D716}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{p}$, to obtain

An important step is to define the asymptotic form of the Weyl operator, which will be used quite often to understand and interpret the leading-order results of its operation on a certain variable. As shown below, this asymptotic form depends on a Taylor expansion of the dispersion relation, and therefore (at least conceptually) should be defined around $\boldsymbol{k}_{0}\neq 0$, since derivatives of the dispersion relation at $\boldsymbol{k}=0$ are singular. In order to obtain the asymptotic form of the Weyl operator, the Fourier transform of the dispersion relation around $\boldsymbol{k}_{0}$ is replaced by the Taylor expansion of the dispersion relation around that point:

which, after Fourier transform with respect to $\overline{\boldsymbol{k}}$, reduces to

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 $\boldsymbol{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 $\unicode[STIX]{x1D719}=A\exp [\text{i}(\boldsymbol{k}_{0}\boldsymbol{\cdot }\boldsymbol{x}-\unicode[STIX]{x1D714}_{0}t)/\unicode[STIX]{x1D716}]$, and consequently the operation $\unicode[STIX]{x1D714}\unicode[STIX]{x1D719}$ 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 $\unicode[STIX]{x1D719}=A(\boldsymbol{x})\exp [(S-\text{i}\unicode[STIX]{x1D714}_{0}t)/\unicode[STIX]{x1D716}]$. 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_{\boldsymbol{x}}$, operates only on the exponent of $\unicode[STIX]{x1D719}$. This term is the dispersion relation, equation (2.1). Terms of $O(\unicode[STIX]{x1D716})$ 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_{\boldsymbol{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_{\boldsymbol{x}}$, operates only on the exponent of $\unicode[STIX]{x1D719}$. 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, $\unicode[STIX]{x1D713}$.

## 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 Reference Dingemans1997). 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, $\boldsymbol{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 $\unicode[STIX]{x1D714}_{0}$. In this case, the surface variables of the considered wave obey to the following form:

where $B_{0}=\text{i}A_{0}(\unicode[STIX]{x1D714}_{0}-\boldsymbol{U}\boldsymbol{\cdot }(D_{\boldsymbol{x}}S))/g$, as follows from the linear relation between $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D702}$ (e.g. Dingemans Reference Dingemans1997):

The wavenumber, $D_{\boldsymbol{x}}S$, is constant, and the constant amplitude, $A_{0}$, assumed to be a ‘proper’ (e.g. Lapidoth Reference Lapidoth2017) Gaussian random variable, namely $\langle A_{0}^{2}\rangle =0$. The corresponding form of $\unicode[STIX]{x1D713}$ 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}+\text{i}aA_{0}$ and $D_{0}=ga^{-1}B_{0}^{\ast }+\text{i}aA_{0}^{\ast }$, and, following (A 7), $a=\sqrt{\unicode[STIX]{x1D70E}(D_{x}S)}$. Using these starting points, it is now aimed to show that (B 1) produces the correct magnitude of $D_{\boldsymbol{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 $\unicode[STIX]{x1D714}_{0}=\unicode[STIX]{x1D714}(D_{\boldsymbol{x}}S)$. Substituting this result into the definition of $B_{0}$, which becomes $B_{0}=\text{i}A_{0}\unicode[STIX]{x1D70E}(D_{\boldsymbol{x}}S)/g$, provides the necessary result that $D_{0}=0$ (this is necessary, since the second term of $\unicode[STIX]{x1D713}$ is not a solution of (B 1)), whereas the value of $C_{0}$ is then given by $C_{0}=2\text{i}A_{0}\sqrt{\unicode[STIX]{x1D70E}(D_{\boldsymbol{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 Reference Holmes2012), and assuming that the bathymetry and the current are not time dependent, the surface velocity potential and elevation, $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D702}$, can be defined to leading order in $\unicode[STIX]{x1D716}$ by (B 2), where now $A_{0}$ and $S$ are functions of $\boldsymbol{x}$. The corresponding leading-order definition of $\unicode[STIX]{x1D713}$ is given by (B 4), but now, following (A 8), $a=\sqrt{\unicode[STIX]{x1D70E}(\boldsymbol{x},D_{\boldsymbol{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 $\unicode[STIX]{x1D714}_{0}=\unicode[STIX]{x1D714}(\boldsymbol{x},D_{\boldsymbol{x}}S)$, which can be written as

Using this result, $B_{0}$ is obtained to be $B_{0}=\text{i}A_{0}\unicode[STIX]{x1D70E}(\boldsymbol{x},D_{\boldsymbol{x}}S)/g$, leading to the same necessary result as before, that $D_{0}=0$, whereas $C_{0}$ is given by $C_{0}=2\text{i}A_{0}\sqrt{\unicode[STIX]{x1D70E}(\boldsymbol{x},D_{\boldsymbol{x}}S)}$. The remaining unknown, $A_{0}$, is found through the following $O(\unicode[STIX]{x1D716})$ 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 Reference Dingemans1997), one should check that $\unicode[STIX]{x1D70C}\langle |\unicode[STIX]{x1D713}|^{2}\rangle =\unicode[STIX]{x1D70C}\langle |C_{0}|^{2}\rangle$ indeed defines the mean action density. To see this, one may calculate $\langle |\unicode[STIX]{x1D713}|^{2}\rangle$ directly, using the definition of $\unicode[STIX]{x1D713}$, 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, $(\unicode[STIX]{x1D70E}\unicode[STIX]{x1D719})_{0}$, is defined as $(\unicode[STIX]{x1D70E}\unicode[STIX]{x1D719})_{0}=\unicode[STIX]{x1D70E}(\boldsymbol{x},D_{\boldsymbol{x}}S)\unicode[STIX]{x1D719}$. 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(\boldsymbol{x})$ and $\boldsymbol{U}(\boldsymbol{x})$ (e.g. § 6.3 in Van Groesen and Molenaar (Reference Van Groesen and Molenaar2007)). Therefore, dividing the expression in the square brackets by the local intrinsic frequency, $\unicode[STIX]{x1D70E}$, leads to the definition of the mean action density (Bretherton & Garrett Reference Bretherton and Garrett1968).

## 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(\boldsymbol{x},\boldsymbol{k},t)$, the following example is discussed (note that here the formulation returns to be written in terms of the original spatial coordinates, $\boldsymbol{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:

where $A_{n}$ and $B_{n}$ are complex random amplitudes. The amplitude $B_{n}$ is given by $B_{n}=\text{i}A_{n}\unicode[STIX]{x1D70E}_{n}/g$, as follows from the linear relation between $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D702}$ (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, $\unicode[STIX]{x1D713}$, which for this example is given by

where $\unicode[STIX]{x1D6FF}(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 $\unicode[STIX]{x1D6FF}[k-(k_{1}+k_{2})/2]$. Following the definition of $k_{3}$, there are two such components, $2\unicode[STIX]{x1D70E}_{3}\langle |A_{3}|^{2}\rangle /g$ and $2\sqrt{\unicode[STIX]{x1D70E}_{1}}\sqrt{\unicode[STIX]{x1D70E}_{2}}\{\langle A_{1}A_{2}^{\ast }\rangle \exp [\text{i}(k_{1}-k_{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 $\unicode[STIX]{x1D70E}_{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, $\unicode[STIX]{x1D713}$, and the mean energy density, $m_{0}$, as given in (2.4), recalling that $a(\boldsymbol{x},D_{\boldsymbol{x}})$ is a Weyl operator associated with the square root of the intrinsic angular frequency, $\unicode[STIX]{x1D70E}^{1/2}(\boldsymbol{x},\boldsymbol{k})$. Following the definition of the Weyl operator, the expression in (2.4) can be written as

By substituting the Fourier transform of the correlation function, $\unicode[STIX]{x1D6E4}=\langle \unicode[STIX]{x1D713}(\boldsymbol{x}+\boldsymbol{p}_{1},t)\unicode[STIX]{x1D713}^{\ast }(\boldsymbol{x}-\boldsymbol{p}_{2},t)\rangle$, the following expression is obtained:

Integrating the above expression with respect to $\boldsymbol{q}_{1}$ and $\boldsymbol{q}_{2}$ leads to

By assuming that $\unicode[STIX]{x1D70E}$ is independent of $\boldsymbol{x}$, equation (C 6) reduces to the following expression:

where the change of variables $\boldsymbol{k}_{1}=\boldsymbol{k}+\boldsymbol{k}^{\prime }/2$ and $\boldsymbol{k}_{2}=\boldsymbol{k}-\boldsymbol{k}^{\prime }/2$ was applied.

Otherwise, equation (C 6) can be integrated once more. Now the integration is performed with respect to $\boldsymbol{p}_{1}$ and $\boldsymbol{p}_{2}$, leading to the following result:

Ultimately, if terms of $O(\unicode[STIX]{x1D716})$ 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:

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, $\boldsymbol{k}$, and around each spatial location, $\boldsymbol{x}$. In addition, one should distinguish, at each location $\boldsymbol{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, $\hat{\unicode[STIX]{x1D6E4}}(\boldsymbol{k}^{\prime },\boldsymbol{k},t)\exp (\text{i}\boldsymbol{k}^{\prime }\boldsymbol{\cdot }\boldsymbol{x})$, are replaced by the Fourier transform of the Wigner distribution around $\boldsymbol{x}$:

where $f(\boldsymbol{x},\boldsymbol{k}^{\prime },\boldsymbol{k})=\sqrt{\unicode[STIX]{x1D70E}(\boldsymbol{x},\boldsymbol{ k}+\boldsymbol{ k}^{\prime }/2)}\sqrt{\unicode[STIX]{x1D70E}(\boldsymbol{x},\boldsymbol{ k}-\boldsymbol{ k}^{\prime }/2)}$. Then, assuming that the Wigner distribution around $\boldsymbol{x}$ can be expressed as a Taylor expansion, and after integrating over $\overline{\boldsymbol{x}}$, the following approximation of $m_{0}$ is obtained:

where, due to the symmetry of $f$ around $\boldsymbol{k}^{\prime }=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(\unicode[STIX]{x1D707})$; see the definition of $\unicode[STIX]{x1D707}$ in § 2.3) to the cross-correlation components that are stored in $\boldsymbol{k}$. Ultimately, in the numerical examples, $m_{0}$ is evaluated up to second order using the following expression:

## Appendix D. The evolution equation for the Wigner distribution

This appendix presents a more detailed derivation of the transport equation for the Wigner distribution, equation (2.20), based on Weyl’s rule of association. The starting point of the derivation is the definition of the Weyl operator, equation (A 3). As a first step, the definition of the Weyl operator is used to define the operator that operates on the correlation function in (2.15) as

Then, the operator can be organized such that the exponential functions being generated due to the disentanglement of the exponential operators (for details, see § 2.4 in Cohen (Reference Cohen2012)) will cancel each other, leading to the following expression:

Following the above expression and the operator correspondences, $f(\boldsymbol{x}^{\prime })\leftrightarrow f(D_{\boldsymbol{k}})$ and $f(D_{\boldsymbol{x}^{\prime }})\leftrightarrow f(\boldsymbol{k})$, the corresponding operator that operates on the Wigner distribution in (2.16) is defined as

where the exponential operators that depend on $D_{\boldsymbol{k}}$ and $D_{\boldsymbol{x}}$ (the third exponential and the fourth exponential in the integral on the right-hand side of (D 3)) can be written as external operators outside of the integral, resulting in the following formulation:

which is exactly the exponential form of the operator as presented in (2.17). For the next steps, it will be convenient to write (D 4) as

which when operates on the Wigner distribution, leading to the following equation:

Finally, according to the assumption of small $\unicode[STIX]{x1D707}$ (see details in § 2.3), the exponential operator is defined through Taylor series. By approximating the exponential operator to first order in $\unicode[STIX]{x1D707}$, equation (D 6) becomes

which is exactly the operator shown in (2.19).

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, $\unicode[STIX]{x1D6E4}(\boldsymbol{x},\boldsymbol{x}^{\prime },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 $\hat{\unicode[STIX]{x1D714}}(\boldsymbol{q},\boldsymbol{k})\exp (\text{i}\boldsymbol{q}\boldsymbol{\cdot }\boldsymbol{x})$ is replaced by the Fourier transform of $\unicode[STIX]{x1D714}$ around $\boldsymbol{x}$, represented by $\hat{\unicode[STIX]{x1D714}}(\boldsymbol{q},\boldsymbol{k},\boldsymbol{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_{\boldsymbol{x}}$, $N_{\boldsymbol{k}}$, $N_{\overline{\boldsymbol{x}}}$, $N_{\boldsymbol{q}}$ (where $\overline{\boldsymbol{x}}=\boldsymbol{x}^{\prime }/2$). These grids are constructed using the following spatial and spectral steps: $\unicode[STIX]{x0394}x$, $\unicode[STIX]{x0394}k$, $\unicode[STIX]{x0394}\overline{x}$, $\unicode[STIX]{x0394}q$. The value of $\unicode[STIX]{x0394}k$ is chosen according to the standard deviation, $S_{d}^{(k)}$, of the incoming wave spectrum as $\unicode[STIX]{x0394}k=S_{d}^{(k)}/\unicode[STIX]{x1D6FC}$, where $\unicode[STIX]{x1D6FC}\geqslant 1$ serves as a resolution factor. Additionally, to ease the computation of the source term $S_{QC}$, the value of $\unicode[STIX]{x0394}k$ is selected such that $\unicode[STIX]{x0394}k=\unicode[STIX]{x0394}q/2$ (this selection prevents the need to perform interpolation in the calculation of the integral in (2.21)).

Next, the choice of $\unicode[STIX]{x0394}q$ is explained. This choice stems from the fact that for any realistic sea state, the correlation function around a certain point, $\boldsymbol{x}$, will effectively have a compact support in $|\boldsymbol{x}^{\prime }|<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 $\unicode[STIX]{x0394}\hat{\unicode[STIX]{x1D714}}$ and $W$ (and their derivatives) over the grid $N_{\boldsymbol{q}}$. This is done without introducing any discretization error if $\unicode[STIX]{x0394}q\leqslant 4\unicode[STIX]{x03C0}/L_{c}$ (for details that also include the additional treatment required to compute the discrete version of $\unicode[STIX]{x0394}\hat{\unicode[STIX]{x1D714}}$, see Smit *et al.* (Reference Smit, Janssen and Herbers2015*a*)). If $\unicode[STIX]{x0394}q$ is chosen such that $\unicode[STIX]{x0394}q=4\unicode[STIX]{x03C0}/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 $\unicode[STIX]{x0394}k$ as $L_{c}=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FC}/S_{d}^{(k)}$, which is consistent with the expected order that should characterize the standard deviation of the envelope of the correlation function ($O(1/S_{d}^{(k)})$).

The choice of $\unicode[STIX]{x0394}\overline{x}$ is argued in a similar way to the selection of $\unicode[STIX]{x0394}q$, but now knowledge about the boundaries of $W$ over $N_{\boldsymbol{k}}$ at a certain point $\boldsymbol{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_{\boldsymbol{k}}$ can change significantly over $N_{\boldsymbol{x}}$. For accuracy, the necessary $\unicode[STIX]{x0394}\overline{x}$ is evaluated in accordance with the boundaries introduced by $N_{\boldsymbol{k}}$. A more economical selection of $\unicode[STIX]{x0394}\overline{x}$, which also introduces an acceptable error, is described by Smit *et al.* (Reference Smit, Janssen and Herbers2015*b*). Note that by introducing $\unicode[STIX]{x0394}\overline{x}$, the summation in $S_{QC}$ becomes limited to the region $[-q_{max},q_{max}]$, where $q_{max}=\unicode[STIX]{x03C0}/\unicode[STIX]{x0394}\overline{x}$.

Finally, the derivation of $\unicode[STIX]{x0394}x$ is considered. Its value should be selected according to the characteristic variation length of $W$ over $\boldsymbol{x}$, and according to the adopted scheme for treating the spatial derivatives. Here, the second-order upwind scheme (see Hirsch Reference Hirsch2007) is used. The local error introduced by this scheme is of $O[(\unicode[STIX]{x0394}x\unicode[STIX]{x1D707}/L)^{3}]$, where $L$ is the characteristic wavelength ($L/\unicode[STIX]{x1D707}$ represents the characteristic length of wave interference at the considered location). Therefore, $\unicode[STIX]{x0394}x$ should be chosen to be small enough so that the global error due to such a magnitude of local error would be acceptable.

## References

Ardhuin, F., Gille, S. T., Menemenlis, D., Rocha, C. B., Rascle, N., Chapron, B., Gula, J. & Molemaker, J.2017 Small-scale open ocean currents have large effects on wind wave heights. J. Geophys. Res. 122 (6), 4500–4517.10.1002/2016JC012413

Ardhuin, F., O’reilly, W., Herbers, T. & Jessen, P.2003 Swell transformation across the continental shelf. Part I. Attenuation and directional broadening. J. Phys. Oceanogr. 33 (9), 1921–1939.

Belibassakis, K., Gerostathis, T. P. & Athanassoulis, G.2011 A coupled-mode model for water wave scattering by horizontal, non-homogeneous current in general bottom topography. Appl. Ocean Res. 33 (4), 384–397.

Besieris, I. M.1985 Wave-kinetic method, phase-space path integrals, and stochastic wave propagation. J. Opt. Soc. Am. A 2 (12), 2092–2099.

Besieris, I. M. & Tappert, F. D.1976 Stochastic wave-kinetic theory in the Liouville approximation. J. Math. Phys. 17 (5), 734–743.10.1063/1.522971

Booij, N., Ris, R. C. & Holthuijsen, L. H.1999 A third-generation wave model for coastal regions: 1. Model description and validation. J. Geophys. Res. 104 (C4), 7649–7666.

Bowen, A. J.1969 Rip currents: 1. Theoretical investigations. J. Geophys. Res. 74 (23), 5467–5478.10.1029/JC074i023p05467

Bretherton, F. P. & Garrett, C. J. R.1968 Wavetrains in inhomogeneous moving media. Proc. R. Soc. Lond. A 302 (1471), 529–554.

Chawla, A., Özkan-Haller, H. T. & Kirby, J. T.1998 Spectral model for wave transformation and breaking over irregular bathymetry. ASCE J. Waterway Port Coastal Ocean Engng 124 (4), 189–198.10.1061/(ASCE)0733-950X(1998)124:4(189)

Chen, Q., Dalrymple, R. A., Kirby, J. T., Kennedy, A. B. & Haller, M. C.1999 Boussinesq modeling of a rip current system. J. Geophys. Res. 104 (C9), 20617–20637.10.1029/1999JC900154

Cohen, L.2012 The Weyl Operator and Its Generalization. Springer Science & Business Media.

Craik, A. D. & Leibovich, S.1976 A rational model for langmuir circulations. J. Fluid Mech. 73 (3), 401–426.

Deigaard, R. et al. 1992 Mechanics of Coastal Sediment Transport, vol. 3. World Scientific Publishing Company.

Dingemans, M. W.1997 Water Wave Propagation Over Uneven Bottoms, vol. 13. World Scientific.

Dyhr-Nielsen, M. & Sørensen, T.1970 Some sand transport phenomena on coasts with bars. Coast. Engng Proc. 12, 855–865.

Fedele, F., Brennan, J., De León, S. P., Dudley, J. & Dias, F.2016 Real world ocean rogue waves explained without the modulational instability. Sci. Rep. 6, 27715.10.1038/srep27715

Group, T. W.1988 The WAM model: a third generation ocean wave prediction model. J. Phys. Oceanogr. 18 (12), 1775–1810.

Hirsch, C.2007 Numerical Computation of Internal and External Flows: The Fundamentals of Computational Fluid Dynamics. Elsevier.

Hlawatsch, F. & Flandrin, P.1997 The interference structure of the wigner distribution and related time-frequency signal representations. In The Wigner Distribution Theory and Applications in Signal Processing, pp. 59–133. Elsevier.

Holmes, M. H.2012 Introduction to Perturbation Methods, vol. 20. Springer Science & Business Media.

Janssen, T., Herbers, T. & Battjes, J.2008 Evolution of ocean wave statistics in shallow water: refraction and diffraction over seafloor topography. J. Geophys. Res. 113 (C3).

Janssen, T. T. & Herbers, T.2009 Nonlinear wave statistics in a focal zone. J. Phys. Oceanogr. 39 (8), 1948–1964.

Kirby, J. T. & Dalrymple, R. A.1986 An approximate model for nonlinear dispersion in monochromatic wave propagation models. Coast. Engng 9 (6), 545–561.

Krasitskii, V. P.1994 On reduced equations in the hamiltonian theory of weakly nonlinear surface waves. J. Fluid Mech. 272, 1–20.

Lapidoth, A.2017 A Foundation in Digital Communication. Cambridge University Press.10.1017/9781316822708

Longuet-Higgins, M. S.1970 Longshore currents generated by obliquely incident sea waves. Part 1. J. Geophys. Res. 75 (33), 6778–6789.

Mapp, G. R., Welch, C. S. & Munday, J. C.1985 Wave refraction by warm core rings. J. Geophys. Res. 90 (C4), 7153–7162.10.1029/JC090iC04p07153

McWilliams, J. C.2016 Submesoscale currents in the ocean. Proc. R. Soc. Lond. A 472 (2189), 20160117.

McWilliams, J. C.2018 Surface wave effects on submesoscale fronts and filaments. J. Fluid Mech. 843, 479–517.

Metzger, J. J., Fleischmann, R. & Geisel, T.2014 Statistics of extreme waves in random media. Phys. Rev. Lett. 112 (20), 203903.10.1103/PhysRevLett.112.203903

Papoulis, A. & Pillai, S. U.2002 Probability, Random Variables, and Stochastic Processes. Tata McGraw-Hill Education.

Poje, A. C., Özgökmen, T. M., Lipphardt, B. L., Haus, B. K., Ryan, E. H., Haza, A. C., Jacobs, G. A., Reniers, A., Olascoaga, M. J., Novelli, G. et al. 2014 Submesoscale dispersion in the vicinity of the deepwater horizon spill. Proc. Natl Acad. Sci. USA 111 (35), 12693–12698.

Reniers, A. & Battjes, J.1997 A laboratory study of longshore currents over barred and non-barred beaches. Coast. Engng 30 (1–2), 1–21.

Rice, S. O.1945 Mathematical analysis of random noise. Bell Syst. Tech. J. 24 (1), 46–156.

Roelvink, D., Reniers, A., Van Dongeren, A., de Vries, J. v. T., McCall, R. & Lescinski, J.2009 Modelling storm impacts on beaches, dunes and barrier islands. Coast. Engng 56 (11–12), 1133–1152.

Ruessink, B., Miles, J., Feddersen, F., Guza, R. & Elgar, S.2001 Modeling the alongshore current on barred beaches. J. Geophys. Res. 106 (C10), 22451–22463.

Smit, P. & Janssen, T.2013 The evolution of inhomogeneous wave statistics through a variable medium. J. Phys. Oceanogr. 43 (8), 1741–1758.

Smit, P., Janssen, T. & Herbers, T.2015a Stochastic modeling of coherent wave fields over variable depth. J. Phys. Oceanogr. 45 (4), 1139–1154.

Smit, P., Janssen, T. & Herbers, T.2015b Stochastic modeling of inhomogeneous ocean waves. Ocean Model. 96, 26–35.

Soong, T. T.1973 Random Differential Equations in Science and Engineering. Elsevier.

Stive, M. J. & De Vriend, H. J.1994 Shear stresses and mean flow in shoaling and breaking waves. Coast. Engng Proc. 24, 594–608.

Tolman, H. L.1991 A third-generation model for wind waves on slowly varying, unsteady, and inhomogeneous depths and currents. J. Phys. Oceanogr. 21 (6), 782–797.

Van Groesen, E. & Molenaar, J.2007 Continuum Modeling in the Physical Sciences, vol. 13. Siam.

Van Rijn, L. C.1993 Principles of Sediment Transport in Rivers, Estuaries and Coastal Seas, vol. 1006. Aqua Publications.

Vellinga, P.1982 Beach and dune erosion during storm surges. Coast. Engng 6 (4), 361–387.

Vincent, C. L. & Briggs, M. J.1989 Refraction–diffraction of irregular waves over a mound. J. Waterway Port Coastal Ocean Engng 115 (2), 269–284.

Yoon, S. B. & Liu, P. L.-F.1989 Interactions of currents and weakly nonlinear water waves in shallow water. J. Fluid Mech. 205, 397–419.

Zijlema, M. & van der Westhuysen, A. J.2005 On convergence behaviour and numerical accuracy in stationary swan simulations of nearshore wind wave spectra. Coastal Engng 52 (3), 237–256.