# Numerical investigation of turbulence of surface gravity waves

Zhou Zhang1  Yulin Pan1<EMAIL_ADDRESS>1Department of Naval Architecture
and Marine Engineering, University of Michigan, Ann Arbor, MI 48109, USA

###### Abstract

In this paper, we numerically study the wave turbulence of surface gravity
waves in the framework of Euler equations of the free surface. The purpose is
to understand the variation of the scaling of the spectra with wavenumber $k$
and energy flux $P$ at different nonlinearity levels under different
forcing/free-decay conditions. For all conditions (free decay, narrow- and
broadband forcing) we consider, we find that the spectral forms approach wave
turbulence theory (WTT) solution $S_{\eta}\sim k^{-5/2}$ and $S_{\eta}\sim
P^{1/3}$ at high nonlinearity levels. With the decrease of nonlinearity level,
the spectra for all cases become steeper, with the narrow-band forcing case
exhibiting the most rapid deviation from WTT. To interpret these spectral
variations, we further investigate two hypothetical and disputable mechanisms
about bound waves and finite-size effect. Through a tri-coherence analysis, we
find that the finite-size effect is present in all cases, which is responsible
for the overall steepening of the spectra and reduced capacity of energy flux
at lower nonlinearity levels. The fraction of bound waves in the domain
generally decreases with the decrease of nonlinearity level, except for the
narrow-band case, which exhibits a transition at some critical nonlinearity
level below which a rapid increase is observed. This increase serves as the
main reason for the fastest deviation from WTT with the decrease of
nonlinearity in the narrow-band forcing case.

###### keywords:

Authors should not enter keywords on the manuscript, as these must be chosen
by the author during the online submission process and will then be added
during the typesetting process (see
http://journals.cambridge.org/data/relatedlink/jfm-keywords.pdf for the full
list)

##  1 Introduction

The normal state of the ocean surface is characterized by a large number of
waves at difference scales subject to nonlinear interactions in the presence
of wind forcing and viscous dissipation. In such a state referred to as wave
turbulence, a continuous surface elevation spectrum is usually developed with
a power-law form in the inertial range, as a result of the energy cascade
through scales. In the framework of weak turbulence theory (WTT), the wave
spectra can be analytically computed based on the assumptions of weak
nonlinearity, infinite domain and phase stochasticity, leading to the so-
called Kolmogorov-Zakharov (KZ) spectra (Zakharov & Filonenko, 1967). For
surface gravity waves, the omnidirectional KZ wavenumber spectra has the form

$S_{\eta}(k)\propto P^{1/3}k^{-5/2},$ (1)

where $P$ is the energy flux to small scales and $k$ is the wavenumber.

Since the assumptions in WTT are difficult to satisfy in finite facilities,
experimental attempts to verify WTT often show deviations from (1). In terms
of the scaling $S_{\eta}(k)\sim k^{\alpha}$ (or its frequency counterpart),
different values of $\alpha$ are observed in experiments under different
conditions (e.g., forcing with different amplitudes and bandwidths), and these
findings are sometimes in disagreement with one another. For example, Falcon
et al. (2007); Nazarenko et al. (2010); Deike et al. (2015); Denissenko et al.
(2007) show that the spectral slope $\alpha$ depends on the forcing condition
and approaches (1) at high (or certain) forcing amplitudes. In contrast,
insensitivity of $\alpha$ to forcing amplitude is reported in Issenmann &
Falcon (2013); Aubourg & Mordant (2016); Herbert et al. (2010), but with
spectral slopes steeper than (1). In addition, Cobelli et al. (2011) shows
that the observed spectra also depend strongly on the forcing bandwidth. The
situation for the scaling between $S_{\eta}(k)$ and $P$ is more elusive, with
Falcon et al. (2007); Issenmann & Falcon (2013) reporting a scaling
$S_{\eta}\sim P$ in disagreement with (1). However, their measurement of $P$
is based on the mean power injected by the wave maker, which may lead to
inconsistency with the concept of energy flux due to the broad-scale
dissipation (Deike et al., 2014; Pan & Yue, 2015). Using the same measurement
of $P$, Cobelli et al. (2011) further suggests that $P$ and $P^{1/3}$ scalings
can be realized for respectively broadband and narrow-band forcings.

The inconsistencies in experiments (with WTT and with one another) are usually
attributed to factors including finite-size effect, bound waves and coherent
structures. First, finte-size effect (e.g. Pushkarev & Zakharov, 2000; Lvov et
al., 2006; Nazarenko, 2006) occurs at low nonlinerity level when the nonlinear
broadening is not sufficient to overcome the discreteness of $k$ caused by the
finite size of the facility. This is in contrast to the continous $k$
configuration in deriving (1). Second, bound waves can be considered as wave
components not satisfying the dispersion relation, generated from harmonics
(i.e., non-resonant interactions) or distortion of the carrier wave (Herbert
et al., 2010; Plant et al., 1999, 2004). It is found in Michel et al. (2018);
Campagne et al. (2019) that bound waves are dominant at high frequencies and
are likely responsible for the deviation of the measurements from WTT (and its
dependence on forcing amplitudes). Third, coherent structures (such as rogue
waves and wave breaking) occurring at high nonlinearity levels are not
described by WTT, thus may lead to spectra different from (1). Candidate
theories to model such spectra include the Phillips spectra (Phillips, 1958)
and the Kuznetsov spectra (Kuznetsov, 2004), which are however out of the
scope of the present paper. Finally, all experiments involve other
complexities such as reflection from boundaries, broad-scale dissipation and
surface tension effect which inevitably affect the spectra to some extent.

In numerical simulations, we are able to have a better control of the wave
field by precisely specifying the forcing/dissipation and implementation of
periodic boundary conditions. This offers us a clean environment to study wave
turbulence at various conditions free of the complexities that are present in
experiments. Existing work include Dyachenko et al. (2004); Lvov et al. (2006)
for forcing turbulence and Onorato et al. (2002); Yokoyama (2004) for free-
decay turbulence of gravity waves in the context of Euler equations. While all
these work report a scaling $S_{\eta}(k)\sim k^{-5/2}$ consistent with (1),
the simulation (in each of them) is conducted at a single nonlinearity level,
and therefore is not capable of resolving/understanding the sensitivity of
spectra to various conditions and their scaling with $P$. On the other hand,
free-decay turbulence is studied for capillary waves in Pan & Yue (2014),
which reveals steepened spectra with the decrease of nonlinearity level.
However, the finding cannot be naively applied to gravity waves, because the
spectral behavior at low nonlinearity critically depends on the discrete
resonant manifold (Hrabski & Pan, 2020) which has not been characterized for
surface gravity waves.

In this work, we conduct a numerical study on the spectral properties of
gravity wave turbulence at different forcing (in terms of bandwidths and
amplitudes) and free-decay (with relatively broadband initial data)
conditions. The purpose is to elucidate the mechanisms underlying the spectral
variation under different conditions. In particular, we will focus on the
hypothetical mechanisms of finite-size effect and bound waves, and will leave
the study of coherent structures to our future work which directly simulates
the two-phase Navier-Stokes equations (since only some of the coherent
structures can be simulated in the framework of Euler equations). We also
remark that while the presented findings are clear in our numerical
simulations, applying them to explain the aforementioned experimental
observations will necessarily require the considerations of further
complexities in experiments.

The outline and some main findings of the paper are as follows. In §2, we
present the numerical setup of the Euler equations under forcing and free-
decay conditions. In §3, we show the numerical results including the scaling
of the wave spectra with $k$ and $P$ at different nonlinearity levels. It is
found that the WTT solution is approached at high nonlinearity levels for all
conditions (free-decay, narrow- and broadband forcing). The spectra deviate
from WTT as the nonlinearity level decreases with the largest deviation rate
observed in the narrow-band forcing case. Mechanisms leading to the variation
of spectra with nonlinearity levels are discussed in terms of bound waves and
finite-size effect. Through a tri-coherence analysis we find that finite-size
effect is present at low nonlinearities for all cases, responsible for the
overall steepening of the spectra and reduced energy flux capacity. The
fraction of bound waves generally decreases with the decrease of nonlinearity,
but exhibits a sharp transition and explains the rapid deviation of spectra
from WTT in the narrow-band forcing case. The conclusions are provided in §4.

##  2 Mathematical formulation

We consider gravity waves on a two-dimensional free surface of an
incompressible, inviscid and irrotational fluid. The flow can be described by
a velocity potential $\phi(\bm{x},z,t)$ satisfying the Laplace’s equation.
Here $\bm{x}=(x,y)$ is the horizontal coordinates, $z$ is the vertical
coordinate and $t$ is time. The surface velocity potential is defined as
$\phi^{S}(\bm{x},t)=\phi(\bm{x},z,t)|_{z=\eta}$, where $\eta(\bm{x},t)$ is the
surface elevation. The evolutions of $\eta$ and $\phi^{S}$ satisfy the Euler
equations in Zakharov form (Zakharov, 1968):

$\eta_{t}+\nabla_{\bm{x}}\eta\cdot\nabla_{\bm{x}}\phi^{S}-(1+\nabla_{\bm{x}}\eta\cdot\nabla_{\bm{x}}\eta)\phi_{z}=0$
(2)
$\phi_{t}^{S}+\eta+\frac{1}{2}\nabla_{\bm{x}}\phi^{S}\cdot\nabla_{\bm{x}}\phi^{S}-\frac{1}{2}(1+\nabla_{\bm{x}}\eta\cdot\nabla_{\bm{x}}\eta)\phi_{z}^{2}=0$
(3)

where $\phi_{z}(\bm{x},t)=\partial\phi/\partial z|_{z=\eta}$ is the surface
vertical velocity, $\nabla_{\bm{x}}=(\partial/\partial x,\partial/\partial y)$
denotes the horizontal gradient. In (3), we assume that the mass and time
units are properly chosen such that density and gravitational acceleration
both take values of unity (Dommermuth & Yue, 1987).

To integrate (2) and (3) in time, we use the higher-order spectral (HOS)
method (Dommermuth & Yue, 1987; West et al., 1987). We use a nonlinearity
order $M=3$ which includes cubic nonlinear terms corresponding to four-wave
interactions in spectral space, consistent with the dominant energy
transferring processes for gravity waves (e.g. Hammack & Henderson, 1993; Mei
et al., 2005). All simulations are conducted in doubly periodic square domain
of size $2\pi\times 2\pi$ corresponding to a fundamental wavenumber $k_{0}=1$,
with spatial resolutions of $512\times 512$ (that are sufficient to capture
the phenomena of interest in this work). To account for dissipation at high
wavenumbers, we add two artificial terms respectively on the right hand sides
of (2) and (3):

$D_{\eta}(k)=\gamma_{k}\eta,$ (4) $D_{\phi^{S}}(k)=\gamma_{k}\phi^{S},$ (5)

with the dissipation coefficient $\gamma_{k}$ defined as

$\gamma_{k}=\gamma_{0}(k/k_{d})^{\nu},$ (6)

where $\gamma_{0}$, $k_{d}$ and $\nu$ are parameters characterizing the
dissipation. This formulation is equivalent to the low-pass filter operation
used in Xiao (2013), developed through phenomenological matching with
measurement of dissipation in experiments. In current work, we use values of
$\gamma_{0}=-50$, $k_{d}=115$ and $\nu=30$. We note that this formulation
provides a sharp transition to dissipation range above $k\approx k_{d}$, which
is essential for us to use dissipation rate to measure energy flux in our
numerical studies (i.e., free of the broad-scale dissipation effect).

In this work, we conduct simulations for a free-decay case, and two cases with
external forcing of broad and narrow bandwidth. For the free-decay case, we
use as initial condition a wavenumber spectra converted from a directional
frequency spectra $S_{D}(\omega,\theta)=D(\theta)G(\omega)$, where $\theta$ is
the directional angle with respect to the positive $x$ direction. The
spreading function $D(\theta)$ characterizes the angular dependence of the
spectra and is chosen to be a cosine-squared function (Tanaka, 2001; Onorato
et al., 2002):

$D(\theta)=\left\\{\begin{array}[]{lc}\frac{2}{\pi}\cos^{2}\theta,&|\theta|\leq\pi/2\\\
0,&|\theta|>\pi/2\end{array}.\right.$ (7)

The frequency spectra $G(\omega)$ is set to the form of a Gaussian function:

$G(\omega)=\frac{B}{\sqrt{2\pi\sigma^{2}}}\exp\left[-\frac{(\omega-\omega_{p})^{2}}{2\sigma^{2}}\right]$
(8)

with $\sigma=0.4$, and $B$ a parameter characterizing the initial effective
steepness $\epsilon\equiv k_{p}H_{s}/2$ (with $H_{s}$ the significant wave
height and $k_{p}$ the peak wavenumber) as a measure of the nonlinearity
level. We use $\omega_{p}=\sqrt{10}$ corresponding to $k_{p}=10$ in the
initial wavenumber spectra.

For the forcing cases, the initial condition is a quiescent water surface and
the waves are generated by forcing with different bandwidths and amplitudes.
Numerically the forcing is modeled by an artificial pressure term
$Q(\theta,k,t)=H(\theta)F(k,t)$ added to the right hand side of equation (3).
The angular cut-off function $H(\theta)$ takes a value of one for
$|\theta|\leq\pi/4$ and zero otherwise. (We have tested that the spreading
angles in all cases do not critically affect the results presented in this
paper, but using a relatively narrower angle in the forcing case favorably
stabilizes the simulation.) The function $F(k,t)$ takes the form (e.g.
Dyachenko et al., 2004; Pan, 2020)

$F(k,t)=\left\\{\begin{array}[]{lc}f_{k}\exp[-Ct+i(\omega_{k}t+R)],&t\leq
T_{c}\\\ f_{k}\exp[-CT_{c}+i(\omega_{k}t+R)],&t>T_{c}\end{array},\right.$ (9)

with

$f_{k}=\left\\{\begin{array}[]{lc}f_{0}\frac{(k-k_{1})(k_{2}-k)}{(k_{1}-k_{2})^{2}},&k_{1}\leq
k\leq k_{2}\\\ 0,&\mathrm{otherwise}\end{array},\right.$ (10)

where $f_{0}$ is the parameter determining the forcing amplitude, $\omega_{k}$
is the angular frequency for wavenumber $k$ calculated from the dispersion
relation, $k_{1}$ and $k_{2}$ are the lower and upper bounds of the forcing
range, and $R$ is a random number uniformly distributed in $[0,2\pi]$ that are
different for each wave mode. In the broadband and narrow-band cases, we use
$[k_{1},k_{2}]=[1,19]$ and $[k_{1},k_{2}]=[9,11]$ respectively, both
corresponding to a peak mode of $k_{p}=(k_{1}+k_{2})/2=10$. As described by
(9), the forcing level decays exponentially in time with rate $C$ before
$t=T_{c}$, and then remains constant. This provides a fast convergence to
stationary turbulence state where the forcing balances the dissipation. We use
values of $T_{c}=500T_{p}$ and $C=\ln 5/T_{c}$ (where $T_{p}$ is the peak
period corresponding to $k_{p}$) which lead to favorable convergence rate in
our study.

##  3 Results

In this section, we present the results from simulations of free-decay
turbulence and forcing turbulence with broad and narrow bandwidths. The former
is conducted at different effective steepness of the initial conditions, and
the latter at different forcing amplitudes, in order to cover a sufficient
range of nonlinearity levels. In the following, we will focus on the scaling
of $S_{\eta}(k)$ with $k$ and $P$ at different nonlinearity levels, and
investigate the mechanisms underlying the variations.

###  3.1 Spectral slopes

We first define the omnidirectional wavenumber spectrum $S_{\eta}(k)$ (we
neglect its $t$ dependence for simplicity in the definition) by

$S_{\eta}(k)=\int_{0}^{2\pi}|\tilde{\eta}(\bm{k})|^{2}kd\theta,$ (11)

where $\bm{k}=(k_{x},k_{y})$ and $k=|\bm{k}|$, $\tilde{\eta}(\bm{k})$ is the
spatial Fourier transform of $\eta(\bm{x})$:

$\tilde{\eta}(\bm{k})=\iint_{[0,2\pi]\times[0,2\pi]}\eta(\bm{x})e^{-i\bm{k}\cdot\bm{x}}d\bm{x}.$
(12)

To demonstrate the (quasi-) stationarity of the spectral evolution, we define
two integral measures $E_{in}$ and $\Psi_{in}$ respectively for the spectra
and compensated spectra (which assigns more weight on the high-wavenumber
part), given by

$E_{in}(t)\equiv\int_{k_{c}}^{k_{d}}S_{\eta}(k,t)dk,$ (13)
$\Psi_{in}(t)\equiv\int_{k_{c}}^{k_{d}}k^{5/2}S_{\eta}(k,t)dk,$ (14)

where $k_{c}=19$ locates within the inertial range as is shown later (also
corresponding to the upper bound of the broadband forcing). We check the
evolutions of $\Psi_{in}/E_{in}$ and $\Psi_{in}$ in free-decay and forcing
cases to characterize their stationary state, where the denominator in the
first quantity is used to account for the slow decay of energy level with time
in the free-decay case.

The results obtained from free-decay cases are presented in figure 1. Figure
1(a) shows evolution of $\Psi_{in}/E_{in}$ with different nonlinearity levels,
measured by the value of $\epsilon\equiv k_{p}H_{s}/2$ evaluated at
$t=500T_{p}$. It can be seen that quasi-stationary states are established for
all cases after $t=500T_{p}$. A typical spectral evolution for the case with
$\epsilon=0.151$ is also shown as an inset to demonstrate the convergence of
the spectrum to a power-law state. The spectra at quasi-stationary states for
different nonlinearity levels (i.e., values of $\epsilon$) are shown in figure
1(b). At high nonlinearity level of $\epsilon=0.151$, we observe a clear
power-law spectrum which has slope $\alpha\approx-5/2$ consistent with WTT
solution (1) in an approximate range of $[10,65]$. We note that the power-law
range ends at $k\approx 65$ which is smaller than $k_{d}=115$, similar to
other numerical simulations with a sharp dissipation cutoff (e.g. Pan & Yue,
2014; Dyachenko et al., 2004), probably due to the interaction of the spectra
with the dissipation range. With the decrease of nonlinearity level
$\epsilon$, the power-law spectrum becomes shorter and steeper, reaching
$\alpha\approx-3.4$ at $\epsilon=0.068$. This steepening of the spectra is in
contrast to the results from Majda-McLaughlin-Tabak (MMT) turbulence with
$\omega=k^{2}$ and quartet resonance (Hrabski & Pan, 2020), mainly because the
latter forms a continuous resonant system (Faou et al., 2016) at low
nonlinearity (we will elaborate this more in §3.4). We also remark that
$\epsilon=0.151$ is about the highest nonlinearity we can reach in the current
HOS context. Previous free-decay simulations (Onorato et al., 2002; Yokoyama,
2004) which result in $\alpha\approx-5/2$ are both conducted at a value of
$\epsilon$ close to 0.15 (respectively 0.15 and 0.14).

Figure 1: (a) Time evolution of $\Psi_{in}/E_{in}$ and (b) stationary spectra
in free-decay cases with $\epsilon=0.068$ (\- $\cdot$ -), $0.106$ ($\cdot\
\cdot\ \cdot$) and $0.151$ ( ) evaluated at $t=500T_{p}$ from different runs.
Inset of (a): wave spectra obtained at $t=0$ (– – –), $100T_{p}$ ( ),
$200T_{p}$ ( ) and $500T_{p}$ ( ) with $\epsilon=0.151$. The theoretical
power-law scaling $k^{-5/2}$ (– – –) is indicated in (b) for reference.

The results from the forcing cases are shown in figure 2. The evolution of
$\Psi_{in}$ in the narrow- and broadband cases with different forcing
amplitudes $f_{0}$ (resulting in different values of $\epsilon$ at stationary
state) are plotted in figure 2(a) and (c), all showing stationary states in
$[1000T_{p},1500T_{p}]$. The stationary power-law spectra for respectively the
two cases are plotted in (b) and (d). For both cases, we observe that the
spectral slopes $\alpha$ approach $-5/2$ at sufficiently high
forcing/nonlinearity (consistent with previous work Dyachenko et al. (2004)).
With the decrease of nonlinearity, the spectra for both cases become steeper,
which are also reported in experiments in Falcon et al. (2007); Nazarenko et
al. (2010); Deike et al. (2015). However, we observe that the spectra at low
nonlinearity clearly show a dependence on the bandwidth of the forcing, with
the one from narrow-band case much steeper (even not showing a power law) than
the one from broadband case. It can also be noticed that, for the narrow-band
case, the spectrum at low nonlinearity level ($\epsilon=0.059$) exhibits
discrete super-harmonic peaks, which will be explained through bound waves in
§3.3.

Figure 2: (a) Time evolution of $\Psi_{in}$ and (b) stationary spectra in
broadband forcing case with $f_{0}=3\times 10^{-7}$, $\epsilon=0.061$ (\-
$\cdot$ -); $f_{0}=1.6\times 10^{-6}$, $\epsilon=0.118$ ($\cdot\ \cdot\
\cdot$); $f_{0}=3.2\times 10^{-6}$, $\epsilon=0.145$ ( ). (c) Time evolution
of $\Psi_{in}$ and (d) stationary spectra in narrow-band forcing case with
$f_{0}=8\times 10^{-7}$, $\epsilon=0.059$ (\- $\cdot$ -); $f_{0}=4\times
10^{-6}$, $\epsilon=0.114$ ($\cdot\ \cdot\ \cdot$); $f_{0}=1\times 10^{-5}$,
$\epsilon=0.148$ ( ). The theoretical power-law scaling $k^{-5/2}$ (– – –) and
boundaries of the forcing ranges ( ) are indicated in (b) and (d) for
reference.

To obtain a complete view of spectral slopes for both the free-decay and
forcing cases, we plot the values of $\alpha$ in all cases as functions of the
nonlinearity level in figure 3. This plot is limited above by the stability of
HOS and below by the existence of a power-law spectra. In practice, we
consider a power-law spectrum not existing if the power-law range is less than
0.5 decade or if the spectrum is dominated by discrete peaks such as the one
with low nonlinearity level in figure 2(d). We see that for all cases, the
spectral slope $\alpha$ approaches WTT value $-5/2$ at sufficiently high
nonlinearity levels of $\epsilon\approx 0.15$. With the decrease of
$\epsilon$, all spectra become steeper but with different steepening rates
especially at relatively low nonlinearity level. It is also clear that the
narrow-band forcing case shows a transition at $\epsilon_{c}\approx 0.11$
below which a very rapid steepening is observed. The mechanisms underlying
these behaviors will be further analyzed in §3.3 and §3.4.

Figure 3: The spectral slope $\alpha$ as functions of effective steepness
$\epsilon$ for the free-decay ($\circ$), broadband forcing ($\Box$) and
narrow-band forcing ($\triangle$) cases, compared to the prediction of WTT
with $\alpha=-5/2$ (– – –).

###  3.2 Energy flux

In this section, we investigate the scaling between the spectral level of
$S_{\eta}(k)$ and the energy flux $P$. For the evaluation of $P$, we can
directly compute the energy transfer due to nonlinear terms in the primitive
equation (Hrabski & Pan, 2020) or use energy dissipation rate as a measure of
$P$ (e.g. Pan & Yue, 2014). For dissipation localized at high wavenumbers
(such as our cases), the two approaches are equivalent (as all energy flux
through the inertial range is dissipated at high wavenumbers in a stationary
state). Therefore, we use dissipation rate as a measure of $P$, which takes
the form (Pushkarev & Zakharov, 1996; Pan & Yue, 2014):

$P=\iint_{\bm{k}}\gamma_{k}(|\tilde{\eta}({\bm{k})}|^{2}+k|\widetilde{\phi^{S}}({\bm{k}})|^{2})d\bm{k},$
(15)

where $\widetilde{\phi^{S}}(\bm{k})$ is the spatial Fourier transform of
$\phi^{S}(\bm{x})$. We note that (15) is slightly less accurate in the free-
decay case since the spectra slowly evolves in the quasi-stationary state.
However, the evolution rate of the spectra is much smaller than the energy
flux so that the associated error is negligible.

The scaling between $E_{in}$ and $P^{1/3}$ is plotted in figure 4 for the
free-decay case and the two forcing cases. The variable $P^{1/3}$ is used for
the horizontal axis so that WTT scaling $E_{in}\sim P^{1/3}$ becomes a
straight line in the figure. We see that, at high nonlinearity level, the
scaling in all cases approach the WTT scaling, although the range of
consistency is longer in the free-decay and the broadband forcing cases. With
the decrease of nonlinearity, the scaling deviates from the WTT prediction
with a smaller value of $P$ for given $E_{in}$. This indicates a reduced
capacity of energy cascade with the reduction of nonlinearity. As nonlinearity
level is further decreased, all curves approach states with $P\rightarrow 0$
and finite $E_{in}$, suggesting the formation of “frozen turbulence”, which is
previously (only) introduced for capillary waves (Pushkarev & Zakharov, 2000;
Pan & Yue, 2014). This result is remarkable because of the existence of exact
resonances for gravity waves on a discrete grid of $\bm{k}$ (rational torus),
in contrast to capillary waves. Before this study it is not clear whether the
exact quartet resonances of gravity waves lead to a cascade as in MMT
turbulence (Hrabski & Pan, 2020). The results here suggest that the energy
cascade should not be expected in spite of the existence of energy transfer
within a small number of resonant quartets.

Figure 4: Plots of the energy density in the inertial range $E_{in}$ as a
function of $P^{1/3}$ for free-decay ($\circ$), broad-band forcing ($\Box$)
and narrow-band forcing ($\triangle$) cases, comparing to the prediction of
WTT $E_{in}\sim P^{1/3}$ (– – –).

In the following, we will study hypothetical physical mechanisms about bound
waves and finite-size effect, which may lead to the spectral behaviors
presented in §3.1 and §3.2.

###  3.3 Bound waves

In this section, we study the effect of bound waves, which is argued as a
major factor influencing spectral behavior in previous experiments (Michel et
al., 2018). Here we generally define bound waves as wave components that do
not satisfy the linear dispersion relation, no matter if they are generated by
carrier-wave distortion (Plant et al., 1999, 2004), super-harmonics (Herbert
et al., 2010; Cobelli et al., 2011; Michel et al., 2018) or other mechanisms
(Campagne et al., 2019; Longuet-Higgins, 1992) proposed in previous
literature. In fact, as we will show in §3.3.1, all bound waves in a periodic-
domain simulation can be interpreted through non-resonant nonlinear
interactions. This definition distinguishes bound waves from free waves which
satisfy the linear dispersion relation (or lie in its vicinity). To separate
bound waves from the wave field, it is necessary to conduct a spatiotemporal
analysis and obtain wavenumber-frequency spectra $S_{\eta}(k,\omega)$, defined
as

$S_{\eta}(k,\omega)=\int_{0}^{2\pi}|\tilde{\eta}(\bm{k},\omega)|^{2}kd\theta,$
(16)

where $\tilde{\eta}(\bm{k},\omega)$ is the spatiotemporal Fourier transform of
$\eta(\bm{x},t)$:

$\tilde{\eta}(\bm{k},\omega)=\iiint_{[0,T_{L}]\times[0,2\pi]\times[0,2\pi]}\eta(\bm{x},t)h_{T}(t)e^{-i(\bm{k}\cdot\bm{x}-\omega
t)}d\bm{x}dt,$ (17)

with $h_{T}(t)$ the Tukey window (Bloomfield, 2004) of length $T_{L}=20T_{p}$,
the time duration of the collected data within $1480T_{p}\sim 1500T_{p}$ at
the stationary state.

####  3.3.1 Generation mechanisms

The wavenumber-frequency spectra $S_{\eta}(k,\omega)$ are plotted in figure 5
from (a) to (d) for respectively narrow-band forcing with low and high
nonlinearity levels, and broadband forcing with low and high nonlinearity
levels. The results for free-decay cases are not shown since they are somewhat
similar to the broadband forcing case. In all cases, we (seem to) observe
significant amount of energy away from the linear dispersion relation (red
curve in each sub-figure), indicating the presence of bound waves in the
simulations. In addition, the plot for narrow-band forcing at low nonlinearity
(figure 5(a)) shows discrete peaks of bound wave components, a visually
distinct pattern from other sub-figures. In retrospect to results in §3.1,
these discrete peaks correspond to the peaks in wavenumber spectra with low
nonlinearity in figure 2(d). As marked by red dots in figure 5(a), these peaks
can be quantified as super-harmonics of the peak mode $(\omega_{p},k_{p})$ of
carrier waves, in the form of

$(\omega_{b},k_{b})=(n\omega_{p},nk_{p}),\quad n=2,3,4,...$ (18)

The phase velocity for all bound wave components, described by (18), can be
computed as $\omega_{b}/k_{b}=\omega_{p}/k_{p}$, which is the same as the
phase velocity of the peak mode of carrier waves. Therefore, the bound waves
can be considered as non-dispersive, i.e., they are “bound” to the carrier
wave as the latter travels. This behavior agrees with some of the classical
views of bound waves described in Lake & Yuen (1978); Plant (2003).

For cases in 5(b)-(d), however, the bound-wave patterns are dramatically
different from that in (a). It is clear from the plots that these bound wave
components are not necessarily non-dispersive, which is consistent with other
general view of bound waves (e.g. Phillips, 1981). In these cases, we observe
several bound-wave branches (in contrast to the main branch of linear
dispersion relation) on the spectra $S_{\eta}(\omega,k)$. These bound-wave
branches can be explained through two different mechanisms discussed below.

The first mechanism corresponds to the super-harmonics of all free waves,
i.e., modes $(\omega,k)$ in the main (carrier) branch of the linear dispersion
relation, satisfying

$(\omega_{b},k_{b})=(m\omega,mk),\quad m=2,3,4,...$ (19)

The second mechanism corresponds to bound waves generated by the (non-
resonant) interactions between the harmonics of peak mode $(\omega_{p},k_{p})$
and an arbitrary free wave $(\omega,k)$ in the main branch, satisfying

$(\omega_{b},k_{b})=(\omega+l\omega_{p},k+lk_{p}),\quad l=\pm 1,\pm 2,\pm
3...$ (20)

The curves described by (19) and (20) are marked in 5(b)-(d) respectively by
black solid lines and brown dashed lines. The two mechanisms co-explain the
super-branches of bound waves on the right of the main branch, as they both
contribute to energy in each of the visible super-branches. The second
mechanism with $l=-1$ further explains the sub-branches of bound waves on the
left of the main branch. We remark that the two mechanisms (19) and (20) are
also separately observed in experiments of Herbert et al. (2010); Cobelli et
al. (2011) and Campagne et al. (2019). Our results here provide a more
comprehensive view of bound waves: we need to consider the property of carrier
waves to distinguish cases in 5(a) and (b)-(d), and consider combined
mechanisms (19) and (20) for explanation in the latter case.

Figure 5: Normalized wavenumber-frequency spectra
$S_{\eta}(\omega,k)/S_{\eta}(\omega_{p},k_{p})$ in log scale for narrow-band
forcing cases with (a) $\epsilon=0.059$, (b) $\epsilon=0.148$ and broadband
forcing cases with (c) $\epsilon=0.061$, (d) $\epsilon=0.145$. The linear
dispersion relation is marked by  . In (a), the peak modes are marked by
$\bullet$ computed from (18) with $n=2,3,4,5$. In (b-d), the lines
corresponding to (19) with $m=2,3$ are indicated by $\cdot\ \cdot\ \cdot$, and
the lines corresponding to (20) with $l=-1,1,2$ are indicated by – – –.

####  3.3.2 Effects on wave spectra

Our next goal is to separate bound waves and free waves to elucidate their
relative importance to the wave turbulence of gravity waves. The algorithm for
the separation is illustrated in figure 6, which shows the definition of a
free-wave finite band (by dashed line) in the vicinity of the linear
dispersion relation. Specifically, this finite band is generated through a
filter (similar to Campagne et al. (2019))

$f(\bm{k},\omega)\equiv\left\\{\begin{array}[]{lc}1,&|\omega-\omega_{k}|\leq\delta\omega\\\
0,&|\omega-\omega_{k}|>\delta\omega\end{array}\right.$ (21)

where $\delta\omega=0.6\omega_{0}$ is a parameter characterizing the width of
the free-wave band, with $\omega_{0}=1$ the fundamental frequency in the
domain. We have tested that this choice of $\delta\omega$ lies in a stationary
range in terms of the results discussed below, i.e., reducing or increasing it
results in respectively insufficient free-wave energy or contamination by
bound-wave branches. For the separation, we apply the filter (21) directly to
$\tilde{\eta}(\bm{k},\omega)$ and obtain the free-wave and bound-wave
components as
$\tilde{\eta}_{f}(\bm{k},\omega)=f(\bm{k},\omega)\tilde{\eta}(\bm{k},\omega)$
and
$\tilde{\eta}_{b}(\bm{k},\omega)=\tilde{\eta}(\bm{k},\omega)-\tilde{\eta}_{f}(\bm{k},\omega)$.
Then we compute the free-wave spectra $S_{\eta}^{f}(k,\omega)$,
$S_{\eta}^{f}(k)$ and bound-wave spectra $S_{\eta}^{b}(k,\omega)$,
$S_{\eta}^{b}(k)$ likewise using (16) and (11). We note that it is also
possible to directly apply the filter (21) to $S_{\eta}(k,\omega)$ for the
separation, but our operation (w.r.t $\tilde{\eta}$) is desirable due to the
study that will be discussed in §3.4.

Figure 6: A typical illustration for the separation of free and bound waves in
a wavenumber-frequency spectrum. The region of
$[\omega_{k}-\delta\omega,\omega_{k}+\delta\omega]$ defined in (21) is marked
by – – –. The region bounded by \- $\cdot$ - corresponds to the inertial
ranges $[\omega_{c},\omega_{d}]\times[k_{c},k_{d}]$ used in (22).

The free-wave and bound-wave wavenumber spectra are shown in figure 7 for all
three (free-decay, braodband and narrow-band forcing) cases at high and low
nonlinearity levels. In most cases, we find that the power-law range of the
spectra are dominated by free-wave components with their energy at least one
order of magnitude higher than those of the bound-wave components. The only
exception is the case with narrow-band forcing at low nonlinearity level shown
in figure 7(e), where bound waves dominate most of the range above $k\approx
17$ exhibiting discrete peaks.

The general behavior of bound waves in figure 7 suggest that they are not the
major factor causing steepening of the spectra except in the narrow-band
forcing case. This point will be made more clear after further quantifying the
fraction of bound waves in the wave field. For this purpose, we define the
bound-wave energy in the inertial range as

$E_{b}\equiv\int_{\omega_{c}}^{\omega_{d}}\int_{k_{c}}^{k_{d}}S_{\eta}^{b}(k,\omega)dkd\omega$
(22)

where $[k_{c},k_{d}]$ and $[\omega_{c}=k_{c}^{1/2},\omega_{d}=k_{d}^{1/2}]$
are a box describing the limits of inertial range as shown in figure 6. We
quantify the fraction of bound waves using the ratio of $E_{b}/E_{t}$, where
$E_{t}$ is the total wave energy in the inertial range computed in a similar
way as (22) but with $S_{\eta}(k,\omega)$ replacing $S_{\eta}^{b}(k,\omega)$.

Figure 7: $S_{\eta}^{f}(k)$ (– – –), $S_{\eta}^{b}(k)$ (\- $\cdot$ -) and
$S_{\eta}(k)$ ( ) in free-decay cases with (a) $\epsilon=0.068$ and (b)
$\epsilon=0.151$; broadband forcing cases with (c) $\epsilon=0.071$ and (d)
$\epsilon=0.145$; narrow-band forcing cases with (e) $\epsilon=0.059$ and (f)
$\epsilon=0.148$. The boundaries of the power-law ranges are indicated by – –
–, except for (e) where discrete peaks are observed.

Figure 8 summarizes the quantity $E_{b}/E_{t}$ at different nonlinearity
levels for all three cases. For the broadband forcing and free-decay cases, we
see that the fraction of bound-wave energy consistently decreases with the
decrease of nonlinearity level. This provides a clear evidence that the
steepening of the spectra with decrease of nonlinearity is not caused by bound
waves in these cases. For narrow-band forcing case, the fraction of bound-wave
energy is consistently larger than those in the other two cases (at all
nonlinearity levels). In addition, with the decrease of nonlinearity, a
transition occurs at a critical nonlinearity level ($\epsilon_{c}\approx
0.11$) below which $E_{b}/E_{t}$ increases with the decrease of $\epsilon$.
(In general, we expect that the critical level $\epsilon_{c}$ to increase with
the decrease of the forcing bandwidth.) The critical level $\epsilon_{c}$ is
consistent with the transition in figure 3 below which a much more rapid
steepening of the spectral slope is observed for the narrow-band forcing case.
We thus conclude that the presence of bound waves explains the largest
steepening rate of the spectra at low nonlinearity levels in the narrow-band
forcing case. Moreover, the significant fraction of bound waves in the narrow-
band forcing case also provide an explanation of the rapid deviation from
$E_{in}\sim P^{1/3}$ shown in figure 4, since bound waves account for energy
which does not follow the WTT cascade pathway.

Figure 8: The proportions of bound wave energy $E_{b}/E_{t}$ as functions of
$\epsilon$ for the free-decay ($\circ$), broadband forcing ($\Box$) and
narrow-band forcing ($\triangle$) cases.

###  3.4 Finite-size effect

In order to further understand the steepening of spectra and reduced energy
flux capacity at low nonlinearity levels, we investigate another hypothetical
mechanism of finite-size effect. In general, finite-size effect arises due to
the violation of the assumption of the infinite domain in WTT (Pushkarev &
Zakharov, 2000; Lvov et al., 2006; Nazarenko, 2006). In the framework of WTT
kinetic equation, the energy cascade is enabled by modes on the continuous
resonant manifold satisfying (for gravity waves)

$\bm{k}_{1}+\bm{k}_{2}=\bm{k}_{3}+\bm{k},$ (23)
$\omega_{1}+\omega_{2}=\omega_{3}+\omega,$ (24)

where $\omega_{i}^{2}=k_{i}\ (i=1,2,3)$. In a finite domain, discreteness of
wavenumber and frequency (imposed by the domain boundary) reduces the manifold
defined by (23) and (24) by limiting the resonances to discrete points. This
discreteness effect can be compensated by the nonlinear broadening which
allows the occurrence of quasi-resonances, characterized by a modification of
(24) as

$|\omega_{1}+\omega_{2}-\omega_{3}-\omega|<\Delta\omega,$ (25)

with $\Delta\omega$ a nonlinear broadening parameter depending on the
nonlinearity level.

One way to quantify the finite-size effect (in particular to measure the
nonlinear broadening and interaction strength) of gravity waves is through a
tricoherence analysis (see e.g. Pan & Yue, 2017, for bi-coherence study for
capillary waves). We first define a tricoherence function

$T(\bm{k},\bm{k}_{1},\bm{k}_{2})\equiv\frac{|\langle\tilde{\eta}_{f}(\bm{k}_{1},t)\tilde{\eta}_{f}(\bm{k}_{2},t)\tilde{\eta}_{f}^{*}(\bm{k}_{3},t)\tilde{\eta}_{f}^{*}(\bm{k},t)\rangle|}{\langle|\tilde{\eta}_{f}(\bm{k}_{1},t)||\tilde{\eta}_{f}(\bm{k}_{2},t)||\tilde{\eta}_{f}^{*}(\bm{k}_{3},t)||\tilde{\eta}_{f}^{*}(\bm{k},t)|\rangle}$
(26)

where $\bm{k}_{3}\equiv-\bm{k}+\bm{k}_{1}+\bm{k}_{2}$, $\langle\cdot\rangle$
denotes the time average (for stationary states), and $*$ denotes the complex
conjugate. $\tilde{\eta}_{f}(\bm{k},t)$ corresponds to the free-wave modes,
and is computed by the spatial Fourier transform of $\eta_{f}(\bm{x},t)$ with
the latter obtained after applying the filter (21) to the wavenumber-frequency
spectra. By definition, the function $T(\bm{k},\bm{k}_{1},\bm{k}_{2})$
continuously varies between 0 and 1, with 1 corresponding to exact resonance
among $\bm{k}$, $\bm{k}_{1}$, $\bm{k}_{2}$ and $\bm{k}_{3}$. In (26), we use
$\tilde{\eta}_{f}$ instead of $\tilde{\eta}$ for the evaluation to remove the
contamination by bound waves since they result in noisy values of $T$ not
lying in the vicinity of the resonant manifold (as found in previous work Pan
& Yue (2017)). This operation is also compatible with findings in §3.3 that
bound waves do not contribute to the general steepening of the spectra (except
in the narrow-band forcing case).

To facilitate the visualization of $T$, we fix $\bm{k}_{1}=(40,40)$ and
$\bm{k}_{2}=(20,-40)$ so that $T(\bm{k})$ with $\bm{k}=(k_{x},k_{y})$ can be
shown by two-dimensional contour plots. Figure 9 shows such plots of
$T(\bm{k})$ in the broadband forcing case at three different nonlinearity
levels. The results for the narrow-band forcing and free-decay cases are not
included due to their similarity to the presented results (the quantification
of broadening and interaction strength for all cases will be presented later).
Also shown in figure 9 is the continuous resonant manifold (red lines)
satisfying (23) and (24), and discrete resonant points (red crosses) with
$\bm{k}$ only taking integer values. In particular, to compute the resonant
manifold we need to consider the wave traveling direction and complex
conjugate relation for real functions, with details presented in Appendix A.

From figure 9 we can see that all significant values of $T$ are concentrated
close to the resonant manifold (as a result of using $\tilde{\eta}_{f}$ in
(26)). It is also clear that the nonlinear broadening is visibly wider at
higher nonlinearity level compared to that at lower nonlinearity level.
Furthermore, the discrete resonant points seem to be sparse on the continuous
manifold, with all four points in each sub-figure of figure 9 as trivial
quartet solutions (either with repetition or symmetry with given vectors of
$\bm{k}_{1}$ and $\bm{k}_{2}$). This sparsity can be further demonstrated by a
numerical analysis which identifies no non-trivial solution with stretching of
the current resonant manifold (see Appendix B). This problem has also been
considered analytically and numerically in Kartashova (2006); Lvov et al.
(2006), which reach a consistent conclusion on the rareness of the discrete
solutions except for two special types of collinear quartets and tridents
(which however do not exist in our case with waves propagating only to the
positive $x$ direction). Moreover, the sparsity of discrete solutions in this
case is in contrast to the situation of $\omega=k^{2}$ MMT dispersion relation
studied in Hrabski & Pan (2020). The difference in the structure of the
discrete resonant solutions is the key for the gravity-wave spectra and MMT
spectra to show completely different behaviors at low nonlinearity, with the
former deviating from WTT spectral slope but the latter approaching WTT
spectral slope.

To further quantify the effect of nonlinearity level to quartet resonances, we
define two measures $\hat{L}_{b}$ and $I$ respectively for the broadening
width and overall interaction strength. For $\hat{L}_{b}$, we define it as a
characteristic width (Pan & Yue, 2017):

$\hat{L}_{b}\equiv\frac{\sum\limits_{\bm{k}}|\hat{\Omega}|T(\bm{k})}{\sum\limits_{\bm{k}}T(\bm{k})},$
(27)

where $|\hat{\Omega}|$ is the normalized frequency mismatch given by
$|\hat{\Omega}|=|\Omega|/k^{-1/2}$ with
$\Omega\equiv\omega_{1}+\omega_{2}-\omega_{3}-\omega_{k}$ and the denominator
$k^{-1/2}$ estimating the frequency discreteness at $k$ associated with the
wavenumber spacing. In equation (27), the summation is for all the grid points
of $\bm{k}$. Thus, $\hat{L}_{b}$ measures the nonlinear broadening around the
exact solutions by the first moment of $T(\bm{k})$.

We further define $I$ as

$I\equiv\sum\limits_{\bm{k}}T(\bm{k}),$ (28)

which measures the overall interaction strength by summing up contributions
from both resonant and quasi-resonant interactions.

The characteristic width $\hat{L}_{b}$ and the interaction strength $I$ are
plotted as functions of $\epsilon$ for all cases in figure 10. We can observe
a general increase of $\hat{L}_{b}$ and $I$ for increasing $\epsilon$,
indicating a wider nonlinear broadening and a larger interaction strength with
more quasi-resonances alleviating the finite-size effect. We remark that this
clear behavior of $\hat{L}_{b}$ and $I$ are only possible to resolve by using
$\tilde{\eta}_{f}$ to evaluate the tri-coherence. The reduction of nonlinear
broadening and interaction strength with decreasing nonlinearity is consistent
with steepening of spectra and reduction of energy flux capacity discussed in
§3.1 and §3.2. Therefore, we conclude that the finite-size effect is a major
contributor to the spectral behaviors at low nonlinearity level (especially
for cases with sufficient spectral bandwidth). Finally, the values of
$\hat{L}_{b}$ and $I$ in the narrow-band forcing case are slightly smaller
than the other two cases for all nonlinearity levels, probably because of the
higher fraction of bound waves (which leads to a larger $\epsilon$ than the
other two cases, but not larger $\hat{L}_{b}$ or $I$).

Figure 9: The tricoherence $T(\bm{k},\bm{k}_{1},\bm{k}_{2})$ with
$\bm{k}_{1}=(40,40)$, $\bm{k}_{2}=(20,-40)$ for broad-band forcing cases with
(a) $\epsilon=0.071$, (b) $\epsilon=0.118$ and (c) $\epsilon=0.145$. The
continuous resonant manifold is plotted by  . The discrete resonant solutions
with $\bm{k}\in(\mathbb{Z},\mathbb{Z})$ are marked by $\times$. The data for
generating these figures are collected from $t=1460T_{p}$ to $t=1500T_{p}$ in
the stationary state.

Figure 10: (a) The characteristic length of nonlinear broadening $\hat{L}_{b}$
and (b) the intensity of nonlinear interactions $I$ as functions of $\epsilon$
for the free-decay ($\circ$), broadband forcing ($\Box$) and narrow-band
forcing ($\triangle$) cases.

##  4 Conclusions

We conduct numerical simulation of Euler equations to study the surface
gravity wave turbulence in three representative conditions, namely free-decay,
narrow-band forcing and broadband forcing turbulence. In all cases, We find
that the scaling of the wave spectra with wavenumber and energy flux both
approach the WTT solution at sufficiently high nonlinearity level. With the
decrease of nonlinearity level, steeper spectra and reduced energy flux
capacity can be observed indicating deviation from WTT, with the largest
deviation rate found in the narrow-band forcing case. Two hypothetical
mechanisms on bound waves and finite-size effect to explain these spectral
variations are investigated. For bound waves, we elucidate their generation
mechanisms through a spatiotemporal analysis (which generalizes the previous
study on this topic) and find that their fraction generally decreases with the
decrease of nonlinearity level except for the narrow-band forcing case. This
suggests that bound waves only account for the rapid deviation from WTT in the
narrow-band forcing case (but not for the other two cases). For finite-size
effect, we perform a tri-coherence analysis and find that both the nonlinear
broadening and interaction strength decrease with the decrease of nonlinearity
level, which accounts for the deviation from WTT at low nonlinearity level in
all cases of our simulation. We finally remark that cautions have to be taken
in applying these numerical findings to experiments due to the additional
complexity inevitably involved in the latter.

## Declaration of Interests

The authors report no conflict of interest.

##  Appendix A Computation of resonant manifold in figure 9

We start by considering a (quasi-)stationary wave field described by

$\eta(\bm{x},t)=\sum_{i}A_{i}e^{i(\bm{\kappa}_{i}\cdot\bm{x}-\omega_{i}t)}+A_{i}^{*}e^{-i(\bm{\kappa}_{i}\cdot\bm{x}-\omega_{i}t)},$
(29)

where $A_{i}=|A_{i}|e^{i\phi_{i}}\in\mathbb{C}$ and
$\omega_{i}\in\mathbb{R}^{+}$. Since the wave fields in our simulations
generally only contain waves traveling to the positive $x$ direction (due to
the setting of initial condition and forcing), for each mode in (29) we have
$\kappa_{ix}\cdot\omega_{i}>0$ so that
$\bm{\kappa}_{i}=(\kappa_{ix},\kappa_{iy})\in(\mathbb{R}^{+},\mathbb{R})$.

Considering $\tilde{\eta}(\bm{k},t)$ with $\bm{k}\in(\mathbb{R},\mathbb{R})$
as the spatial Fourier transform of $\eta(\bm{x},t)$, we can see that

$\tilde{\eta}(\bm{k}_{i},t)=\left\\{\begin{array}[]{lc}A_{i}e^{-i\omega_{i}t},&k_{ix}>0\\\
A_{i}^{*}e^{i\omega_{i}t},&k_{ix}<0\end{array},\right.$ (30)

Substituting (30) to (26), we obtain

$T=|\langle e^{i(-\Omega_{d}t+\Phi_{d})}\rangle|,$ (31)

where $\Phi_{d}=\phi_{1}+\phi_{2}-\phi_{3}-\phi$, and

$\Omega_{d}=\operatorname{sgn}(k_{1x})\omega_{1}+\operatorname{sgn}(k_{2x})\omega_{2}-\operatorname{sgn}(k_{3x})\omega_{3}-\operatorname{sgn}(k_{x})\omega,$
(32)

with

$\operatorname{sgn}(x)=\left\\{\begin{array}[]{lc}1,&x>0\\\
-1,&x<0\end{array},\right.$ (33)

In our case, we have $k_{1x}=40>0$, $k_{2x}=20>0$ and
$k_{3x}=k_{1x}+k_{2x}-k_{x}=60-k_{x}$, so that (32) is reduced to
$\Omega_{d}=\omega_{1}+\omega_{2}-\operatorname{sgn}(60-k_{x})\omega_{3}-\operatorname{sgn}(k_{x})\omega$.
We are interested in the case of $\Omega_{d}=0$, which results in $T=1$
indicating the resonant manifold. Whether $\Omega_{d}=0$ can be realized needs
to be discussed in the following three situations (summarized in figure
11(a)):

  1. 1.

For $k_{x}<0$, we need $\Omega_{d}=\omega_{1}+\omega_{2}-\omega_{3}+\omega=0$
which has no solution (for $\omega\in\mathbb{R}^{+}$) on the $(k_{x},k_{y})$
plane.

  2. 2.

For $0<k_{x}<60$, we need
$\Omega_{d}=\omega_{1}+\omega_{2}-\omega_{3}-\omega=0$ which has solutions
shown as an ellipse bounded by $k_{x}=0$ and $k_{x}=60$ in figure 11(b). This
solutions corresponds to the resonant manifold shown in figure 9.

  3. 3.

For $k_{x}>60$, we have $\Omega_{d}=\omega_{1}+\omega_{2}+\omega_{3}-\omega=0$
which has no solution on the $(k_{x},k_{y})$ plane.

In summary, the resonant manifold shown in figure 9 corresponds to the second
case above, which is the only possibility to have solutions in $\Omega_{d}=0$
for a wave field traveling to the positive $x$ direction.

Figure 11: (a) Sketches of the three situation in discussion. Vectors
$\bm{k}_{1}$ and $\bm{k}_{2}$ are indicated in the figure, $\bm{k}$ and
$\bm{k}_{3}$ are denoted by red and yellow arrows. The two vertical lines
marked the boundaries of situation (i), (ii) and (iii) in the discussion. (b)
The solution in situation (ii) as denoted by red lines in the figure (ellipse
bounded by the two vertical line).

##  Appendix B Sparsity of the discrete resonant solutions

We count the number of discrete resonant solutions as in figure 9, but with
stretched $\bm{k}_{1}$ and $\bm{k}_{2}$. The purpose is to demonstrate the
sparsity of discrete resonances with the expansion of the resonant manifold.
In particular, we define $\bm{k}_{1}=\lambda(40,40)$ and
$\bm{k}_{2}=\lambda(20,-40)$ where $\lambda$ is a factor of stretching, and
compute the number of exact resonances $N$ as a function of $\lambda$. For
$\lambda$ varying from 1 to 10, our numerical searching algorithm gives $N=4$
always, indicating that no non-trivial solution can be identified. While this
numerical study is presented for $\bm{k}_{1}$ and $\bm{k}_{2}$ in particular
directions, we have also tested other directions to confirm that the exact
resonances are indeed sparse for gravity waves.

## References

  * Aubourg & Mordant (2016) Aubourg, Quentin & Mordant, Nicolas 2016 Investigation of resonances in gravity-capillary wave turbulence. Physical Review Fluids 1 (2), 023701.
  * Bloomfield (2004) Bloomfield, Peter 2004 Fourier analysis of time series: an introduction. John Wiley & Sons.
  * Campagne et al. (2019) Campagne, Antoine, Hassaini, Roumaissa, Redor, Ivan, Valran, Thomas, Viboud, Samuel, Sommeria, Joël & Mordant, Nicolas 2019 Identifying four-wave-resonant interactions in a surface gravity wave turbulence experiment. Physical Review Fluids 4 (7), 074801.
  * Cobelli et al. (2011) Cobelli, P, Przadka, Adam, Petitjeans, Philippe, Lagubeau, Guillaume, Pagneux, Vincent & Maurel, Agnes 2011 Different regimes for water wave turbulence. Physical Review Letters 107 (21), 214503.
  * Deike et al. (2014) Deike, Luc, Berhanu, Michael & Falcon, Eric 2014 Energy flux measurement from the dissipated energy in capillary wave turbulence. Physical Review E 89 (2), 023003.
  * Deike et al. (2015) Deike, Luc, Miquel, Benjamin, Gutiérrez, P, Jamin, Timothée, Semin, Benoit, Berhanu, Michael, Falcon, Eric & Bonnefoy, Félicien 2015 Role of the basin boundary conditions in gravity wave turbulence. Journal of Fluid Mechanics 781, 196–225.
  * Denissenko et al. (2007) Denissenko, Petr, Lukaschuk, Sergei & Nazarenko, Sergey 2007 Gravity wave turbulence in a laboratory flume. Physical review letters 99 (1), 014501.
  * Dommermuth & Yue (1987) Dommermuth, Douglas G & Yue, Dick KP 1987 A high-order spectral method for the study of nonlinear gravity waves. Journal of Fluid Mechanics 184, 267–288.
  * Dyachenko et al. (2004) Dyachenko, Alexander I, Korotkevich, Alexander O & Zakharov, Vladimir E 2004 Weak turbulent kolmogorov spectrum for surface gravity waves. Physical review letters 92 (13), 134501\.
  * Falcon et al. (2007) Falcon, Eric, Laroche, Claude & Fauve, Stéphan 2007 Observation of gravity-capillary wave turbulence. Physical review letters 98 (9), 094503.
  * Faou et al. (2016) Faou, Erwan, Germain, Pierre & Hani, Zaher 2016 The weakly nonlinear large-box limit of the 2d cubic nonlinear schrödinger equation. Journal of the American Mathematical Society 29 (4), 915–982.
  * Hammack & Henderson (1993) Hammack, Joseph L & Henderson, DM 1993 Resonant interactions among surface water waves. Annual review of fluid mechanics 25 (1), 55–97.
  * Herbert et al. (2010) Herbert, Eric, Mordant, Nicolas & Falcon, Eric 2010 Observation of the nonlinear dispersion relation and spatial statistics of wave turbulence on the surface of a fluid. Physical review letters 105 (14), 144502.
  * Hrabski & Pan (2020) Hrabski, Alexander & Pan, Yulin 2020 Effect of discrete resonant manifold structure on discrete wave turbulence. Physical Review E 102 (4), 041101.
  * Issenmann & Falcon (2013) Issenmann, Bruno & Falcon, Eric 2013 Gravity wave turbulence revealed by horizontal vibrations of the container. Physical Review E 87 (1), 011001.
  * Kartashova (2006) Kartashova, Elena 2006 Fast computation algorithm for discrete resonances among gravity waves. Journal of low temperature physics 145 (1), 287–295.
  * Kuznetsov (2004) Kuznetsov, Evgenii Aleksandrovich 2004 Turbulence spectra generated by singularities. Journal of Experimental and Theoretical Physics Letters 80 (2), 83–89.
  * Lake & Yuen (1978) Lake, Bruce M & Yuen, Henry C 1978 A new model for nonlinear wind waves. part 1. physical model and experimental evidence. Journal of fluid mechanics 88 (1), 33–62.
  * Longuet-Higgins (1992) Longuet-Higgins, Michael S 1992 Capillary rollers and bores. Journal of Fluid Mechanics 240, 659–679.
  * Lvov et al. (2006) Lvov, Yuri V, Nazarenko, Sergey & Pokorni, Boris 2006 Discreteness and its effect on water-wave turbulence. Physica D: Nonlinear Phenomena 218 (1), 24–35.
  * Mei et al. (2005) Mei, Chiang C, Stiassnie, Michael & Yue, Dick K-P 2005 Theory and applications of ocean surface waves: nonlinear aspects, , vol. 23. World scientific.
  * Michel et al. (2018) Michel, Guillaume, Semin, Benoît, Cazaubiel, Annette, Haudin, Florence, Humbert, Thomas, Lepot, Simon, Bonnefoy, Félicien, Berhanu, Michaël & Falcon, Éric 2018 Self-similar gravity wave spectra resulting from the modulation of bound waves. Physical Review Fluids 3 (5), 054801.
  * Nazarenko (2006) Nazarenko, Sergey 2006 Sandpile behaviour in discrete water-wave turbulence. Journal of Statistical Mechanics: Theory and Experiment 2006 (02), L02002.
  * Nazarenko et al. (2010) Nazarenko, Sergey, Lukaschuk, Sergei, McLelland, Stuart & Denissenko, Petr 2010 Statistics of surface gravity wave turbulence in the space and time domains. Journal of Fluid Mechanics 642, 395–420.
  * Onorato et al. (2002) Onorato, Miguel, Osborne, Alfred Richard, Serio, Marina, Resio, D, Pushkarev, A, Zakharov, Vladimir E & Brandini, C 2002 Freely decaying weak turbulence for sea surface gravity waves. Physical review letters 89 (14), 144501.
  * Pan (2020) Pan, Yulin 2020 High-order spectral method for the simulation of capillary waves with complete order consistency. Journal of Computational Physics 408, 109299.
  * Pan & Yue (2014) Pan, Yulin & Yue, Dick KP 2014 Direct numerical investigation of turbulence of capillary waves. Physical review letters 113 (9), 094501.
  * Pan & Yue (2015) Pan, Yulin & Yue, Dick KP 2015 Decaying capillary wave turbulence under broad-scale dissipation. Journal of Fluid Mechanics 780.
  * Pan & Yue (2017) Pan, Yulin & Yue, Dick KP 2017 Understanding discrete capillary-wave turbulence using a quasi-resonant kinetic equation. Journal of Fluid Mechanics 816.
  * Phillips (1981) Phillips, OM 1981 The dispersion of short wavelets in the presence of a dominant long wave. Journal of Fluid Mechanics 107, 465–485.
  * Phillips (1958) Phillips, Oliver M 1958 The equilibrium range in the spectrum of wind-generated waves. Journal of Fluid Mechanics 4 (4), 426–434.
  * Plant (2003) Plant, William J 2003 A new interpretation of sea-surface slope probability density functions. Journal of Geophysical Research: Oceans 108 (C9).
  * Plant et al. (2004) Plant, William J, Dahl, Peter H, Giovanangeli, Jean-Paul & Branger, Hubert 2004 Bound and free surface waves in a large wind-wave tank. Journal of Geophysical Research: Oceans 109 (C10).
  * Plant et al. (1999) Plant, William J, Keller, William C, Hesany, Vahid, Hara, Tetsu, Bock, Erik & Donelan, Mark A 1999 Bound waves and bragg scattering in a wind-wave tank. Journal of Geophysical Research: Oceans 104 (C2), 3243–3263.
  * Pushkarev & Zakharov (1996) Pushkarev, AN & Zakharov, VE 1996 Turbulence of capillary waves. Physical review letters 76 (18), 3320.
  * Pushkarev & Zakharov (2000) Pushkarev, AN & Zakharov, VE 2000 Turbulence of capillary waves—theory and numerical simulation. Physica D: Nonlinear Phenomena 135 (1-2), 98–116.
  * Tanaka (2001) Tanaka, Mitsuhiro 2001 Verification of hasselmann’s energy transfer among surface gravity waves by direct numerical simulations of primitive equations. Journal of Fluid Mechanics 444, 199\.
  * West et al. (1987) West, Bruce J, Brueckner, Keith A, Janda, Ralph S, Milder, D Michael & Milton, Robert L 1987 A new numerical method for surface hydrodynamics. Journal of Geophysical Research: Oceans 92 (C11), 11803–11824.
  * Xiao (2013) Xiao, Wenting 2013 Study of directional ocean wavefield evolution and rogue wave occurrence using large-scale phase-resolved nonlinear simulations. PhD thesis, Massachusetts Institute of Technology.
  * Yokoyama (2004) Yokoyama, Naoto 2004 Statistics of gravity waves obtained by direct numerical simulation. Journal of Fluid Mechanics 501, 169–178.
  * Zakharov & Filonenko (1967) Zakharov, VE & Filonenko, NN 1967 Energy spectrum for stochastic oscillations of the surface of a liquid. SPhD 11, 881.
  * Zakharov (1968) Zakharov, Vladimir E 1968 Stability of periodic waves of finite amplitude on the surface of a deep fluid. Journal of Applied Mechanics and Technical Physics 9 (2), 190–194.