1 Introduction
The studies of nonlinear waves initiated by Stokes in the middle of the nineteenth century (Stokes Reference Stokes1847) focused on finite-amplitude effects for a single wave train. In this case, the mechanism of nonlinear interaction between waves which appears at third order in the steepness is not readily apparent. Work of Tick (Reference Tick1959), Phillips (Reference Phillips1960) and Longuet-Higgins (Reference Longuet-Higgins1962) investigating the interaction between several wave trains first shed light on the importance of this resonant energy exchange mechanism, and elucidated how waves of small steepness change their amplitudes and frequencies in the presence of other such waves.
The mutual effects on the phase speed of two interacting wave trains were first studied explicitly by Longuet-Higgins & Phillips (Reference Longuet-Higgins and Phillips1962), by means of classical asymptotic expansion techniques. A unified method for treating arbitrarily many modes, or indeed a continuum of modes, is provided by the Zakharov equation (Zakharov Reference Zakharov1968), which was employed by Hogan, Gruman & Stiassnie (Reference Hogan, Gruman and Stiassnie1988) to extend Longuet-Higgins & Phillips’ analysis to capillary-gravity waves. A thorough comparison of the two methods was made by Zhang & Chen (Reference Zhang and Chen1999) for collinear waves in deep water. Madsen & Fuhrman generalized the theory to include the effects of finite depth and an ambient current for bi-directional waves (Reference Madsen and Fuhrman2006), and subsequently for multi-directional, irregular waves (Reference Madsen and Fuhrman2012).
In what follows, an investigation into the nonlinear corrections on the frequencies is undertaken using the Zakharov equation. General forms for the corrected frequency are thus available for continuous one-dimensional (1D) and two-dimensional (2D) energy density spectra, giving an analogue to the Stokes correction for any wave of a spectral sea-state. These are applied to a number of calculated examples in deep water.
The theoretical background is developed in § 2. Results for some simple cases, as well as for Pierson–Moskowitz (PM) spectra in one and two dimensions are presented in § 3. The main conclusion on the importance of these results to wave forecasting is indicated in § 4. Some comments on the robustness of the results are also given in appendix A, which deals with deterministic and random quartets of waves.
2 Theoretical background
For a wave field with a discretized spectrum, Zakharov’s equation with the complex amplitudes
$B_{n}=B(\boldsymbol{k}_{n},t)$
reads
$$\begin{eqnarray}\text{i}\frac{\text{d}B_{n}}{\text{d}t}=\mathop{\sum }_{p,q,r=1}^{N}T_{npqr}\unicode[STIX]{x1D6FF}_{np}^{qr}\text{e}^{\text{i}\unicode[STIX]{x1D6E5}_{npqr}t}B_{p}^{\ast }B_{q}B_{r},\quad n=1,2,\ldots ,N,\end{eqnarray}$$
where
$T_{npqr}=T(\boldsymbol{k}_{n},\boldsymbol{k}_{p},\boldsymbol{k}_{q},\boldsymbol{k}_{r})$
are the kernels defined in Mei, Stiassnie & Yue (Reference Mei, Stiassnie and Yue2018),
$\unicode[STIX]{x1D6FF}_{np}^{qr}$
is a Kronecker delta function such that
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}_{np}^{qr}=\left\{\begin{array}{@{}ll@{}}1\quad & \text{for }\boldsymbol{k}_{n}+\boldsymbol{k}_{p}=\boldsymbol{k}_{q}+\boldsymbol{k}_{r},\\ 0\quad & \text{otherwise,}\end{array}\right. & & \displaystyle\end{eqnarray}$$
and
$\unicode[STIX]{x1D6E5}_{npqr}=\unicode[STIX]{x1D714}_{n}+\unicode[STIX]{x1D714}_{p}-\unicode[STIX]{x1D714}_{q}-\unicode[STIX]{x1D714}_{r}.$
The asterisk denotes a complex conjugate. Here
$\unicode[STIX]{x1D714}_{i}$
are the frequencies (in
$\text{rad}~\text{s}^{-1}$
),
$g$
is the constant acceleration of gravity (taken to be
$9.8~\text{m}~\text{s}^{-2}$
in computations),
$t$
is time (in s) and
$\boldsymbol{k}=(k_{x},k_{y})$
is the wavenumber.
The linear dispersion relation for gravity waves in water of constant depth
$h$
is
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{n}^{2}=g|\boldsymbol{k}_{n}|\tanh (|\boldsymbol{k}_{n}|h).\end{eqnarray}$$
The relation of the complex amplitudes
$B_{n}$
to the free surface elevation
$\unicode[STIX]{x1D702}$
is given, to leading order, by
$$\begin{eqnarray}\unicode[STIX]{x1D702}=\frac{1}{2\unicode[STIX]{x03C0}}\mathop{\sum }_{n}\left(\frac{\unicode[STIX]{x1D714}_{n}}{2g}\right)^{1/2}(B_{n}\text{e}^{\text{i}(\boldsymbol{k}_{n}\boldsymbol{\cdot }\boldsymbol{x}-\unicode[STIX]{x1D714}_{n}t)}+\text{*}).\end{eqnarray}$$
Substituting
$B_{n}=|B_{n}|\exp (\text{i}\arg B_{n})$
into (2.1) and separating the real and imaginary parts:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}|B_{n}|}{\text{d}t}=\mathop{\sum }_{p,q,r}T_{npqr}\unicode[STIX]{x1D6FF}_{np}^{qr}|B_{p}||B_{q}||B_{r}|\sin (\unicode[STIX]{x1D703}_{npqr}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}(\arg B_{n})}{\text{d}t}=-|B_{n}|^{-1}\mathop{\sum }_{p,q,r}T_{npqr}\unicode[STIX]{x1D6FF}_{np}^{qr}|B_{p}||B_{q}||B_{r}|\cos (\unicode[STIX]{x1D703}_{npqr}), & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\unicode[STIX]{x1D703}_{npqr}=\unicode[STIX]{x1D6E5}_{npqr}t-\arg B_{n}-\arg B_{p}+\arg B_{q}+\arg B_{r}.\end{eqnarray}$$
This same substitution in (2.4) yields
$$\begin{eqnarray}\unicode[STIX]{x1D702}=\frac{1}{2\unicode[STIX]{x03C0}}\mathop{\sum }_{n}\left(\frac{\unicode[STIX]{x1D714}_{n}}{2g}\right)^{1/2}(|B_{n}|\text{e}^{\text{i}(\boldsymbol{k}_{n}\boldsymbol{\cdot }\boldsymbol{x}-\unicode[STIX]{x1D714}_{n}t+\arg B_{n})}+\text{*}).\end{eqnarray}$$
From (2.6) it is clear that the wave amplitudes
$a_{n}$
and the frequencies
$\unicode[STIX]{x1D6FA}_{n}$
are
$$\begin{eqnarray}\displaystyle & \displaystyle a_{n}=\frac{1}{\unicode[STIX]{x03C0}}\left(\frac{\unicode[STIX]{x1D714}_{n}}{2g}\right)^{1/2}|B_{n}|, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{n}=\unicode[STIX]{x1D714}_{n}-\frac{\text{d}}{\text{d}t}(\arg B_{n}). & \displaystyle\end{eqnarray}$$
$e_{np}=1$
for
$n=p$
and
$e_{np}=2$
for
$n\neq p,$
and noting that
$\sin (\unicode[STIX]{x1D703}_{npnp})=0$
and
$\cos (\unicode[STIX]{x1D703}_{npnp})=1,$
we find 
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}|B_{n}|}{\text{d}t}=\mathop{\sum }_{p}\mathop{\sum }_{q\neq n}\mathop{\sum }_{r\neq n}T_{npqr}\unicode[STIX]{x1D6FF}_{np}^{qr}|B_{p}||B_{q}||B_{r}|\sin (\unicode[STIX]{x1D703}_{npqr}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}(\arg B_{n})}{\text{d}t}=-\mathop{\sum }_{p}\!e_{np}T_{npnp}|B_{p}|^{2}\!-|B_{n}|^{-1}\mathop{\sum }_{p}\mathop{\sum }_{q\neq n}\mathop{\sum }_{r\neq n}T_{npqr}\unicode[STIX]{x1D6FF}_{np}^{qr}|B_{p}||B_{q}||B_{r}|\cos (\unicode[STIX]{x1D703}_{npqr}).\qquad & \displaystyle\end{eqnarray}$$
In cases without resonant quartets, or without nearly resonant quartets,
$\sin (\unicode[STIX]{x1D703}_{npqr})$
and
$\cos (\unicode[STIX]{x1D703}_{npqr})$
oscillate on a fast time scale, and the integration over time of (2.8) and (2.9) yields at leading order
$$\begin{eqnarray}\displaystyle & \displaystyle |B_{n}(t)|=|B_{n}(0)|, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \arg (B_{n})=-\!\left(\mathop{\sum }_{p}e_{np}T_{npnp}|B_{p}|^{2}\right)t+\arg (B_{n}(0)). & \displaystyle\end{eqnarray}$$
Substituting (2.11) into (2.7b ) gives the discrete nonlinear dispersion relation
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{n}=\unicode[STIX]{x1D714}_{n}+\mathop{\sum }_{p}e_{np}T_{npnp}|B_{p}|^{2}.\end{eqnarray}$$
For a continuous wavenumber energy spectrum
$\unicode[STIX]{x1D6F9}(\boldsymbol{k})$
it is assumed that the number of modes
$N$
tends to infinity. In the limit they become densely distributed over the relevant domain in the wavenumber plane. This passage to a continuum of modes is accompanied by the assumption that the phases, i.e.
$\arg B_{n}$
, to leading order are uncorrelated and uniformly distributed over
$(-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0}]$
; they are denoted below by
$\unicode[STIX]{x1D700}(\boldsymbol{k})$
.
The limit of a continuous wavenumber energy spectrum
$\unicode[STIX]{x1D6F9}(\boldsymbol{k})$
is approached by considering a square grid of wavenumbers with spacing
$\text{d}\boldsymbol{k},$
so that
$$\begin{eqnarray}a_{n}^{2}/2=\unicode[STIX]{x1D6F9}_{n}=\unicode[STIX]{x1D6F9}(\boldsymbol{k})\,\text{d}\boldsymbol{k}.\end{eqnarray}$$
Using (2.7a,b
) and (2.13) and taking the limit
$\text{d}\boldsymbol{k}\rightarrow 0,$
(2.6) is rewritten as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}(\boldsymbol{x},t)=\int \cos (\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}-\unicode[STIX]{x1D6FA}t+\unicode[STIX]{x1D700}(\boldsymbol{k}))\sqrt{2\unicode[STIX]{x1D6F9}(\boldsymbol{k})\,\text{d}\boldsymbol{ k}} & & \displaystyle\end{eqnarray}$$
(see Kinsman Reference Kinsman1984, chap. 8), and (2.12) is rewritten as
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}(\boldsymbol{k})=\unicode[STIX]{x1D714}(\boldsymbol{k})+4\unicode[STIX]{x03C0}^{2}g\int e_{\boldsymbol{k},\boldsymbol{k}_{1}}T(\boldsymbol{k},\boldsymbol{k}_{1},\boldsymbol{k},\boldsymbol{k}_{1})\frac{\unicode[STIX]{x1D6F9}(\boldsymbol{k}_{1})}{\unicode[STIX]{x1D714}(\boldsymbol{k}_{1})}\,\text{d}\boldsymbol{k}_{1}, & & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle e_{\boldsymbol{k},\boldsymbol{k}_{1}}=\left\{\begin{array}{@{}ll@{}}1\quad & \text{for }\boldsymbol{k}=\boldsymbol{k}_{1},\\ 2\quad & \text{otherwise.}\end{array}\right. & & \displaystyle\end{eqnarray}$$
Equation (2.14b
) is the continuous nonlinear dispersion relation. The kernel
$T(\boldsymbol{k}_{i},\boldsymbol{k}_{j},\boldsymbol{k}_{i},\boldsymbol{k}_{j})$
for finite depth is given in (3.9a,b) and (4.11) of Stiassnie & Gramstad (Reference Stiassnie and Gramstad2009). The applicability of (2.14a
), (2.14b
) depends on the validity of the transition from (2.8), (2.9) to (2.10), (2.11). For continuous spectra the above transition must be based on the assumption that the ensemble averages of
$\sin (\unicode[STIX]{x1D703}_{npqr})$
and
$\cos (\unicode[STIX]{x1D703}_{npqr})$
remain zero, which is true only as long as the phases
$\arg B_{n}$
maintain their random independence. Note that this is exactly the condition used to ensure the applicability of the kinetic equation, and implies a separation between the time scale of nonlinear interaction and that of phase mixing, more details of which may be found in L’vov & Nazarenko (Reference L’vov and Nazarenko2010, § 2) or Zakharov et al. (Reference Zakharov, L’vov and Falkovich1992, § 2.1.4–5). In particular, this means that the spectrum must be sufficiently broad for this approach to be strictly applicable.
Although the above derivation is valid for finite depth, taking the deep-water limit simplifies the kernels dramatically. In this case, substituting
$\unicode[STIX]{x1D714}_{n}^{2}=g|\boldsymbol{k}_{n}|$
into the kernel yields
$$\begin{eqnarray}\displaystyle & & \displaystyle T(\boldsymbol{k}_{i},\boldsymbol{k}_{j},\boldsymbol{k}_{i},\boldsymbol{k}_{j})=-\frac{1}{16\unicode[STIX]{x03C0}^{2}(|\boldsymbol{k}_{i}||\boldsymbol{k}_{j}|)^{1/2}}\left[3(|\boldsymbol{k}_{i}||\boldsymbol{k}_{j}|)^{2}\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\,(\boldsymbol{k}_{i}\boldsymbol{\cdot }\boldsymbol{k}_{j})(\boldsymbol{k}_{i}\boldsymbol{\cdot }\boldsymbol{k}_{j}-4(|\boldsymbol{k}_{i}|+|\boldsymbol{k}_{j}|)(|\boldsymbol{k}_{i}||\boldsymbol{k}_{j}|)^{1/2})\nonumber\\ \displaystyle & & \displaystyle \quad +\left.\frac{2(\unicode[STIX]{x1D714}_{i}-\unicode[STIX]{x1D714}_{j})^{2}(\boldsymbol{k}_{i}\boldsymbol{\cdot }\boldsymbol{k}_{j}+|\boldsymbol{k}_{i}||\boldsymbol{k}_{j}|)^{2}}{g|\boldsymbol{k}_{i}-\boldsymbol{k}_{j}|-(\unicode[STIX]{x1D714}_{i}-\unicode[STIX]{x1D714}_{j})^{2}}+\frac{2(\unicode[STIX]{x1D714}_{i}+\unicode[STIX]{x1D714}_{j})^{2}(\boldsymbol{k}_{i}\boldsymbol{\cdot }\boldsymbol{k}_{j}-|\boldsymbol{k}_{i}||\boldsymbol{k}_{j}|)^{2}}{g|\boldsymbol{k}_{i}+\boldsymbol{k}_{j}|-(\unicode[STIX]{x1D714}_{i}+\unicode[STIX]{x1D714}_{j})^{2}}\right]\end{eqnarray}$$
(see Leblanc Reference Leblanc2009). The examples to be treated subsequently are given for deep water only. Note that an early derivation of (2.12) and (2.14b ) was outlined in a conference in 1988 (see Stiassnie Reference Stiassnie and Osborne1991), as well as by Dyachenko & Zakharov (Reference Dyachenko and Zakharov1994) and Zakharov (Reference Zakharov1999). Appendix A provides some further justification of the assumptions underlying (2.12) and (2.14b ), via direct numerical computations for a quartet of waves.
3 Results
3.1 Simple cases, based on (2.12)
The simplest example that shows the corrections to frequency due to nonlinearity is the third-order Stokes wave (Stokes Reference Stokes1847, p. 450), with the well-known frequency correction
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D714}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D714}k^{2}a^{2}, & & \displaystyle\end{eqnarray}$$
where
$a$
is the wave amplitude,
$k$
the wavenumber and
$\unicode[STIX]{x1D714}$
the linear frequency. The collinear bi-modal case, where the spectrum consists of two wave trains with wavenumbers
$k_{a}$
and
$k_{b},$
likewise has a simple form:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{a}=\unicode[STIX]{x1D714}_{a}\left(1+\frac{1}{2}\unicode[STIX]{x1D716}_{a}^{2}+\frac{\unicode[STIX]{x1D714}_{b}}{\unicode[STIX]{x1D714}_{a}}\cdot \frac{k_{a}^{2}}{k_{b}^{2}}\cdot \unicode[STIX]{x1D716}_{b}^{2}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{b}=\unicode[STIX]{x1D714}_{b}\left(1+\frac{1}{2}\unicode[STIX]{x1D716}_{b}^{2}+\frac{\unicode[STIX]{x1D714}_{a}}{\unicode[STIX]{x1D714}_{b}}\cdot \frac{k_{b}}{k_{a}}\cdot \unicode[STIX]{x1D716}_{a}^{2}\right). & \displaystyle\end{eqnarray}$$
$k_{b}>k_{a},$
and
$\unicode[STIX]{x1D716}_{a}=k_{a}a_{a},\unicode[STIX]{x1D716}_{b}=k_{b}a_{b}$
are the slopes of the two waves involved. It is readily observed that the longer wave has a substantial effect on the shorter wave, but not vice versa. Note that when
$\boldsymbol{k}=(k_{x},0)$
and
$\boldsymbol{k}_{1}=(k_{1x},0),$
the Zakharov kernel simplifies to 
$$\begin{eqnarray}T(\boldsymbol{k},\boldsymbol{k}_{1},\boldsymbol{k},\boldsymbol{k}_{1})=\left\{\begin{array}{@{}ll@{}}\displaystyle \frac{k_{x}k_{1x}^{2}}{4\unicode[STIX]{x03C0}^{2}}\quad & \text{for }k_{1x}<k_{x},\\ \displaystyle \frac{k_{x}^{2}k_{1x}}{4\unicode[STIX]{x03C0}^{2}}\quad & \text{for }k_{1x}\geqslant k_{x},\end{array}\right.\end{eqnarray}$$
(see equation (4.18) in Zakharov Reference Zakharov, Banner and Grimshaw1992). Expressions similar to (3.2a ) and (3.2b ), albeit with more modes, have been obtained by Qi et al. (Reference Qi, Wu, Liu, Kim and Yue2018), and used to improve the performance of their nonlinear phase-resolved reconstruction method.
3.2 Pierson–Moskowitz spectra, based on (2.14b )
A more realistic case, with energy distributed over a continuum of modes, consists of taking a unidirectional PM spectrum, which is given in terms of wavenumber
$k$
by
$$\begin{eqnarray}\unicode[STIX]{x1D6F9}(k)=\frac{0.004}{k^{3}}\exp (-0.554g^{2}U^{-4}k^{-2}),\quad \text{for }k>0,\end{eqnarray}$$
where
$U$
is the wind velocity at ten metres above the sea surface, blowing in the
$x$
-direction. In the unidirectional setting we assume that all wavevectors
$\boldsymbol{k}$
have no component in the
$y$
-direction, and are written as scalar wavenumbers, so that
$\boldsymbol{k}_{i}=(k_{i},0)$
is denoted by
$k_{i}$
. Substituting (3.4) with (3.3) into (2.14b
) and integrating over
$k_{1}$
gives
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}(k)=\unicode[STIX]{x1D714}(k)+2k\left(\int _{0}^{k}k_{1}\unicode[STIX]{x1D714}(k_{1})\unicode[STIX]{x1D6F9}(k_{1})\,\text{d}k_{1}+k\int _{k}^{\infty }\unicode[STIX]{x1D714}(k_{1})\unicode[STIX]{x1D6F9}(k_{1})\,\text{d}k_{1}\right)\!,\end{eqnarray}$$
which can be evaluated analytically to yield
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}=\unicode[STIX]{x1D714}(1+0.004[\unicode[STIX]{x1D705}^{-0.25}\unicode[STIX]{x1D6E4}(0.25,\unicode[STIX]{x1D705})+\unicode[STIX]{x1D705}^{-0.75}\unicode[STIX]{x1D6FE}(0.75,\unicode[STIX]{x1D705})]),\end{eqnarray}$$
where
$\unicode[STIX]{x1D705}=0.554g^{2}U^{-4}k^{-2},$
and
$\unicode[STIX]{x1D6E4},\unicode[STIX]{x1D6FE}$
are incomplete gamma functions; see Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik1980, p. 941). Equation (3.6) enables one to obtain a value for the relative correction of the phase speed due to the presence of other waves in the spectrum. Denoting the linear/nonlinear phase speed by
$c/C$
respectively,
$c=\unicode[STIX]{x1D714}/k$
and
$C=\unicode[STIX]{x1D6FA}/k$
give
$$\begin{eqnarray}\frac{C-c}{c}=0.00405(\unicode[STIX]{x1D705}^{-0.25}\unicode[STIX]{x1D6E4}(0.25,\unicode[STIX]{x1D705})+\unicode[STIX]{x1D705}^{-0.75}\unicode[STIX]{x1D6FE}(0.75,\unicode[STIX]{x1D705})).\end{eqnarray}$$
The dashed line in figure 1 shows
$(C-c)/c$
as a function of the dimensionless wavenumber
$k/k_{p},$
where
$k_{p}=0.6657(g/U^{2})$
is the peak of the PM spectrum (3.4). For
$k=k_{p}$
the nonlinear correction for the phase speed
$(C-c)/c$
is only 0.4 %, but this grows to 12.8 % for
$k=100k_{p}.$
For 2D sea states, the directional spreading of the spectrum may be controlled by
$$\begin{eqnarray}D(\unicode[STIX]{x1D703}_{1})=\left\{\begin{array}{@{}ll@{}}A_{1}\cos ^{m}(\unicode[STIX]{x1D703}_{1})\quad & \text{for }|\unicode[STIX]{x1D703}_{1}|\leqslant \unicode[STIX]{x03C0}/2,\\ 0\quad & \text{for }|\unicode[STIX]{x1D703}_{1}|>\unicode[STIX]{x03C0}/2,\end{array}\right.\end{eqnarray}$$
where
$A_{1}=\unicode[STIX]{x1D6E4}(1+m/2,0)/(\sqrt{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D6E4}(0.5+m/2,0))$
is a normalization coefficient (see Holthuijsen (Reference Holthuijsen2007, p. 164), and note that
$\unicode[STIX]{x1D6E4}(x,0)$
denotes the complete gamma function). The directional spectrum is then given by
$D(\unicode[STIX]{x1D703}_{1})\unicode[STIX]{x1D6F9}(k_{1}),$
where
$\unicode[STIX]{x1D6F9}(k_{1})$
is the unidirectional PM spectrum given in (3.4).
Plot of relative corrections to phase speed for a unidirectional and 2D PM spectrum with directional spreading coefficient
$m=25,$
on a logarithmic scale. Given is
$(C-c)/c$
for a unidirectional spectrum (black, dashed line) and values for angles
$\unicode[STIX]{x1D703}$
from
$0^{\circ }$
(the principal direction of wave propagation of the 2D spectrum) to
$67.5^{\circ }.$
Substituting
$D(\unicode[STIX]{x1D703}_{1})\unicode[STIX]{x1D6F9}(k_{1})$
into (2.14b
) yields
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}(\boldsymbol{k})=\unicode[STIX]{x1D714}(\boldsymbol{k})+4\unicode[STIX]{x03C0}^{2}\int _{-\unicode[STIX]{x03C0}/2}^{\unicode[STIX]{x03C0}/2}\int _{0}^{\infty }e_{\boldsymbol{ k},\boldsymbol{k}_{1}}T(\boldsymbol{k},\boldsymbol{k}_{1},\boldsymbol{k},\boldsymbol{k}_{1})\frac{\unicode[STIX]{x1D714}(k_{1})}{k_{1}}\unicode[STIX]{x1D6F9}(k_{1})D(\unicode[STIX]{x1D703}_{1})\,\text{d}k_{1}\,\text{d}\unicode[STIX]{x1D703}_{1},\end{eqnarray}$$
where
$\boldsymbol{k}_{1}=(k_{1}\cos (\unicode[STIX]{x1D703}_{1}),k_{1}\sin (\unicode[STIX]{x1D703}_{1})),$
and
$\boldsymbol{k}=(k\cos (\unicode[STIX]{x1D703}),k\sin (\unicode[STIX]{x1D703})).$
This allows for a numerical calculation of the frequency correction for waves
$\boldsymbol{k}$
travelling at an angle to the principal direction of wave propagation (at
$\unicode[STIX]{x1D703}_{1}=0$
). The relative corrections to the phase speeds derived from the corrected frequency are given in figure 1, where
$\unicode[STIX]{x1D6FA}(\boldsymbol{k})=\unicode[STIX]{x1D6FA}(k,\unicode[STIX]{x1D703}).$
With most of the energy concentrated about
$\unicode[STIX]{x1D703}_{1}=0$
(given
$m=25$
), the corrections are largest for small angles (nearly collinear waves, or
$\unicode[STIX]{x1D703}$
close to zero). For a wave travelling at a large angle to the principal direction of wave propagation, the corrections are less significant.
4 Relevance to wave forecasting, based on (2.14b )
According to equation (6.4.9) in Holthuijsen (Reference Holthuijsen2007) the balance of spectral energy given by a slowly varying wavenumber spectrum
$\unicode[STIX]{x1D6F9}(\boldsymbol{k},t,\boldsymbol{x})$
in deep water is given by
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6F9}}{\unicode[STIX]{x2202}t}+\boldsymbol{c}_{g}(\boldsymbol{k})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{x}}\unicode[STIX]{x1D6F9}=S(\boldsymbol{k},t,\boldsymbol{x}).\end{eqnarray}$$
The source term
$$\begin{eqnarray}S=S_{in}+S_{nl}+S_{d}\end{eqnarray}$$
represents wave generation by the wind, nonlinear wave interaction (by a stochastic generalization of Zakharov’s equation) and dissipation by wave breaking, respectively. In (4.1),
$\boldsymbol{c}_{g}(\boldsymbol{k})$
is the linear group velocity:
$$\begin{eqnarray}\boldsymbol{c}_{g}(\boldsymbol{k})=\frac{\text{d}\unicode[STIX]{x1D714}}{\text{d}|\boldsymbol{k}|}\frac{\boldsymbol{k}}{|\boldsymbol{k}|}=\frac{1}{2}\sqrt{\frac{g}{|\boldsymbol{k}|}}\frac{\boldsymbol{k}}{|\boldsymbol{k}|}.\end{eqnarray}$$
A more accurate formulation of the convective left-hand side of (4.1) would be to replace the linear group velocity
$\boldsymbol{c}_{g}$
by the nonlinear counterpart:
$$\begin{eqnarray}\boldsymbol{C}_{g}=\frac{\text{d}\unicode[STIX]{x1D6FA}}{\text{d}|\boldsymbol{k}|}\boldsymbol{\cdot }\frac{\boldsymbol{k}}{|\boldsymbol{k}|}.\end{eqnarray}$$
We call the magnitude of the linear/nonlinear group velocity
$\boldsymbol{c}_{g}/\boldsymbol{C}_{g}$
the linear/nonlinear group speed and refer to it by
$c_{g}=|\boldsymbol{c}_{g}|$
or
$C_{g}=|\boldsymbol{C}_{g}|$
, respectively. In the case of a 1D PM spectrum, an analytical expression can be found by differentiating (3.6) with respect to
$k$
to yield
$\text{d}\unicode[STIX]{x1D6FA}/\text{d}k,$
and thus
$$\begin{eqnarray}\frac{C_{g}-c_{g}}{c_{g}}=0.0081\unicode[STIX]{x1D705}^{-0.25}\unicode[STIX]{x1D6E4}(0.25,\unicode[STIX]{x1D705})+0.0162\unicode[STIX]{x1D705}^{-0.75}\unicode[STIX]{x1D6FE}(0.75,\unicode[STIX]{x1D705}).\end{eqnarray}$$
Plot of relative corrections to group speed for a unidirectional and 2D PM spectrum with directional spreading coefficient
$m=25,$
on a logarithmic scale. Given is the correction
$(C_{g}-c_{g})/c_{g}$
for a unidirectional spectrum (black, dashed line) and values for angles
$\unicode[STIX]{x1D703}$
from
$0^{\circ }$
(the principal direction of wave propagation of the 2D spectrum) to
$67.5^{\circ }.$
From the black dashed line in figure 2, which is calculated from (4.5), the nonlinear correction to the group speed is about 1.5 % for
$k=k_{p}$
, and it grows to 7.7 % for
$k=10k_{p},$
and to 26.8 % for
$k=100k_{p}.$
For a directional PM spectrum, with directional spreading given in (3.8), it is possible to give the analogous relative corrections to the group speed using (4.4) together with (2.14b
). Figure 2 gives this relative correction to the group speed for a directional spreading coefficient
$m=25,$
for waves
$\boldsymbol{k}$
travelling at various angles to the principal direction of wave propagation. We believe that these changes are significant enough to be included in modern wave-forecasting software.
Acknowledgement
M.S. is grateful for the support of the Israel Science Foundation, grant no. 261/17.
Appendix A. Comparisons for four waves, based on (2.12)
If only a single quartet of waves is considered, the Zakharov equations written for these waves
$\boldsymbol{k}_{a},\boldsymbol{k}_{b},\boldsymbol{k}_{c}$
and
$\boldsymbol{k}_{d}$
become rather simple, and can be solved either analytically (Stiassnie & Shemer Reference Stiassnie and Shemer2005) in terms of elliptic functions, or by numerical integration. This allows for a direct comparison of the corrected frequency (2.12) with (2.7b
). For the latter,
$\arg B_{n}$
are calculated from (2.8) to (2.9). Since (2.12) is a leading-order approximation valid in the absence of resonance, or for random waves, these are the scenarios of interest. For quartets close to resonance, the large energy exchanges introduce further fluctuations in the frequencies, albeit on the slow Zakharov equation time scale
$T_{2}=\unicode[STIX]{x1D716}^{2}t.$
The quartet
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \boldsymbol{k}_{a}=[0.9806,-0.1961],\quad \boldsymbol{k}_{b}=[0.9806,0.1961],\\ \displaystyle \boldsymbol{k}_{c}=[1.2903-\unicode[STIX]{x1D707},0.2747+\unicode[STIX]{x1D707}],\quad \boldsymbol{k}_{d}=[0.6709+\unicode[STIX]{x1D707},-0.2747-\unicode[STIX]{x1D707}]\end{array}\right\} & & \displaystyle\end{eqnarray}$$
is used for computations, where
$\unicode[STIX]{x1D707}$
allows for a move out of resonance. Figure 3 shows numerical results for this quartet, where
$\unicode[STIX]{x1D714}_{a}=\unicode[STIX]{x1D714}_{b}=3.1305,\unicode[STIX]{x1D714}_{c}=3.7057,\unicode[STIX]{x1D714}_{d}=2.4189$
are the linear frequencies in
$\text{rad}~\text{s}^{-1}$
, given as dotted lines. The corresponding values of steepness are
$\unicode[STIX]{x1D700}_{a}=0.15,\unicode[STIX]{x1D700}_{b}=0.12,\unicode[STIX]{x1D700}_{c}=0.08,\unicode[STIX]{x1D700}_{d}=0.03,$
and
$\unicode[STIX]{x1D707}=-0.15,$
so that
$\min \unicode[STIX]{x1D6E5}_{abcd}/\unicode[STIX]{x1D714}_{i}=0.0547>\max \unicode[STIX]{x1D700}_{i}^{2}$
. The initial values are then given by
$|B_{i}(0)|=\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{i}\sqrt{2g/\unicode[STIX]{x1D714}_{i}}/|k_{i}|$
(see Stiassnie & Shemer Reference Stiassnie and Shemer2005) while the initial arguments are taken to be zero.
Frequencies
$\unicode[STIX]{x1D6FA}_{i}$
versus time
$t$
(in seconds) for each of the four modes
$\boldsymbol{k}_{a},\boldsymbol{k}_{b},\boldsymbol{k}_{c},\boldsymbol{k}_{d}$
in the non-resonant quartet. Solid lines are computed from (2.7b
) by solving the Zakharov equation. Dashed lines are computed from (2.12) using initial amplitudes
$|B_{i}(0)|=\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{i}\sqrt{2g/\unicode[STIX]{x1D714}_{i}}/|\boldsymbol{k}_{i}|.$
Dotted lines are the uncorrected, linear frequencies
$\unicode[STIX]{x1D714}_{i}$
.
The solid lines show
$\unicode[STIX]{x1D6FA}_{i},i\in \{a,b,c,d\}$
obtained from solutions to the Zakharov equations (2.8), (2.9), whereby
$\unicode[STIX]{x1D6FA}_{i}$
is determined from (2.7b
). The dashed lines show the leading-order Stokes correction given in (2.12), which assumes constant magnitudes
$|B_{i}|$
(see (2.10)). The leading-order Stokes correction is in all cases in very good agreement with solutions obtained from the Zakharov equation, and presents a substantial correction to the linear frequencies.
Frequencies
$\unicode[STIX]{x1D6FA}_{i}$
versus time
$t$
(in seconds) for each of the four modes
$\boldsymbol{k}_{a},\boldsymbol{k}_{b},\boldsymbol{k}_{c},\boldsymbol{k}_{d}$
in the resonant quartet. Solid lines are computed from (2.7b
) after solving the Zakharov equation and averaging over 2000 realizations with random phases and amplitudes. Dashed lines are computed from (2.12) using initial amplitudes
$|B_{i}(0)|=\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{i}\sqrt{2g/\unicode[STIX]{x1D714}_{i}}/|\boldsymbol{k}_{i}|.$
Dotted lines are the uncorrected, linear frequencies
$\unicode[STIX]{x1D714}_{i}.$
It is also possible to consider the case of a resonant quartet with random phases and amplitudes. When
$\unicode[STIX]{x1D707}=0$
in the above quartet, this is very close to exact resonance, with
$\max \unicode[STIX]{x1D6E5}_{abcd}/\unicode[STIX]{x1D714}\approx 10^{-6}<\min \unicode[STIX]{x1D716}_{i}^{2},$
with the same values of
$\unicode[STIX]{x1D700}_{i}$
as above. Owing to this resonance, any single solution to the Zakharov equation shows large oscillations in the amplitudes
$|B_{i}(t)|$
with time, and the leading order corrected frequencies of (2.12), based on the initial values
$|B_{i}(0)|,$
are not in good agreement with those obtained from solving the full Zakharov equation. However, the averaged results shown in figure 4 do agree very well with this leading-order correction.
In figure 4, the phases are chosen randomly and uniformly distributed over
$(0,2\unicode[STIX]{x03C0}]$
; the amplitudes are chosen from a Rayleigh distribution with mean
$\unicode[STIX]{x1D707}_{i}$
equal to the initial amplitude used for the deterministic case, i.e.
$\unicode[STIX]{x1D707}_{i}=|B_{i}(0)|.$
An average is taken over 2000 realizations for the resonant quartet. The observed agreement demonstrates the applicability of (2.12) for a simple, discretized spectrum containing four modes.
