Numerical investigation of turbulence of surface gravity waves

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.


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 η (k) ∝ P 1/3 k −5/2 , (1.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.1). In terms of the scaling S η (k) ∼ k α (or its frequency counterpart), different values of α 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 α depends on the forcing condition and approaches (1.1) at high (or certain) forcing amplitudes. In contrast, insensitivity of α to forcing amplitude is reported in Issenmann & Falcon (2013); Aubourg & Mordant (2016); Herbert et al. (2010), but with spectral slopes steeper than (1.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 η (k) and P is more elusive, with Falcon et al. (2007); Issenmann & Falcon (2013) reporting a scaling S η ∼ P in disagreement with (1.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 broadscale 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.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. 1999Plant et al. , 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.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 η (k) ∼ k −5/2 consistent with (1.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 freedecay (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 finitesize 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.

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 φ(x, z, t) satisfying the Laplace's equation. Here x = (x, y) is the horizontal coordinates, z is the vertical coordinate and t is time. The surface velocity potential is defined as φ S (x, t) = φ(x, z, t)| z=η , where η(x, t) is the surface elevation. The evolutions of η and φ S satisfy the Euler equations in Zakharov form (Zakharov 1968): where φ z (x, t) = ∂φ/∂z| z=η is the surface vertical velocity, ∇ x = (∂/∂x, ∂/∂y) denotes the horizontal gradient. In (2.2), 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.1) and (2.2) 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π × 2π corresponding to a fundamental wavenumber k 0 = 1, with spatial resolutions of 512×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.1) and (2.2): with the dissipation coefficient γ k defined as where γ 0 , k d and ν 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 γ 0 = −50, k d = 115 and ν = 30. We note that this formulation provides a sharp transition to dissipation range above k ≈ 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 (ω, θ) = D(θ)G(ω), where θ is the directional angle with respect to the positive x direction. The spreading function D(θ) characterizes the angular dependence of the spectra and is chosen to be a cosine-squared function (Tanaka 2001;Onorato et al. 2002): The frequency spectra G(ω) is set to the form of a Gaussian function: with σ = 0.4, and B a parameter characterizing the initial effective steepness ≡ 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 ω p = √ 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(θ, k, t) = H(θ)F (k, t) added to the right hand side of equation (2.2). The angular cut-off function H(θ) takes a value of one for |θ| π/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 where f 0 is the parameter determining the forcing amplitude, ω 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π] 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 (2.8), 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.

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 η (k) with k and P at different nonlinearity levels, and investigate the mechanisms underlying the variations.

Spectral slopes
We first define the omnidirectional wavenumber spectrum S η (k) (we neglect its t dependence for simplicity in the definition) by where k = (k x , k y ) and k = |k|,η(k) is the spatial Fourier transform of η(x): To demonstrate the (quasi-) stationarity of the spectral evolution, we define two integral measures E in and Ψ in respectively for the spectra and compensated spectra (which assigns more weight on the high-wavenumber part), given by 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 Ψ in /E in and Ψ 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 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 = 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 ) are shown in figure 1(b). At high nonlinearity level of = 0.151, we observe a clear power-law spectrum which has slope α ≈ −5/2 consistent with WTT solution (1.1) in an approximate range of [10,65]. We note that the power-law range ends at k ≈ 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 , the power-law spectrum becomes shorter and steeper, reaching α ≈ −3.4 at = 0.068. This steepening of the spectra is in contrast to the results from Majda-McLaughlin-Tabak (MMT) turbulence with ω = 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 = 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 α ≈ −5/2 are both conducted at a value of close to 0.15 (respectively 0.15 and 0.14).
The results from the forcing cases are shown in figure 2. The evolution of Ψ in in the narrow-and broadband cases with different forcing amplitudes f 0 (resulting in different values of 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 α 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  (d) stationary spectra in narrow-band forcing case with f 0 = 8 × 10 −7 , = 0.059 (-· -); f 0 = 4 × 10 −6 , = 0.114 (· · ·); f 0 = 1 × 10 −5 , = 0.148 ( ). The theoretical power-law scaling k −5/2 (---) and boundaries of the forcing ranges ( ) are indicated in (b) and (d) for reference. broadband case. It can also be noticed that, for the narrow-band case, the spectrum at low nonlinearity level ( = 0.059) exhibits discrete super-harmonic peaks, which will be explained through bound waves in §3.3.
To obtain a complete view of spectral slopes for both the free-decay and forcing cases, we plot the values of α 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 α approaches WTT value −5/2 at sufficiently high nonlinearity levels of ≈ 0.15. With the decrease of , 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 c ≈ 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.

Energy flux
In this section, we investigate the scaling between the spectral level of S η (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): where φ S (k) is the spatial Fourier transform of φ S (x). We note that (3.5) 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 ∼ 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 → 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 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. 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.

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(Plant et al. , 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 periodicdomain 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 η (k, ω), defined as whereη(k, ω) is the spatiotemporal Fourier transform of η(x, t): 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 ∼ 1500T p at the stationary state.

Generation mechanisms
The wavenumber-frequency spectra S η (k, ω) 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 subfigures. 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 (ω p , k p ) of carrier waves, in the form of (ω b , k b ) = (nω p , nk p ), n = 2, 3, 4, ... (3.8) The phase velocity for all bound wave components, described by (3.8), can be computed as ω b /k b = ω 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 η (ω, k). These boundwave 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 (ω, k) in the main (carrier) branch of the linear dispersion relation, satisfying (ω b , k b ) = (mω, mk), m = 2, 3, 4, ... (3.9) The second mechanism corresponds to bound waves generated by the (non-resonant) interactions between the harmonics of peak mode (ω p , k p ) and an arbitrary free wave (ω, k) in the main branch, satisfying The curves described by (3.9) and (3.10) 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 (3.9) and (3.10) are also separately observed in experiments of Herbert et al. (2010); Cobelli et al. (2011) andCampagne 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 (3.9) and (3.10) for explanation in the latter case.

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 . In (a), the peak modes are marked by • computed from (3.8) with n = 2, 3, 4, 5. In (bd), the lines corresponding to (3.9) with m = 2, 3 are indicated by · · ·, and the lines corresponding to (3.10) with l = −1, 1, 2 are indicated by ---. 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)) where δω = 0.6ω 0 is a parameter characterizing the width of the free-wave band, with ω 0 = 1 the fundamental frequency in the domain. We have tested that this choice of δω 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 (3.11) directly toη(k, ω) and obtain the free-wave and bound-wave components asη f (k, ω) = f (k, ω)η(k, ω) and η b (k, ω) =η(k, ω) −η f (k, ω). Then we compute the free-wave spectra S f η (k, ω), S f η (k) and bound-wave spectra S b η (k, ω), S b η (k) likewise using (3.6) and (3.1). We note that it is also possible to directly apply the filter (3.11) to S η (k, ω) for the separation, but our operation (w.r.tη) is desirable due to the study that will be discussed in §3.4. 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 ≈ 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 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 (3.12) but with S η (k, ω) replacing S b η (k, ω). 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 boundwave 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 ( c ≈ 0.11) below which E b /E t increases with the decrease of . (In general, we expect that the critical level c to increase with the decrease of the forcing bandwidth.) The critical level c is consistent with the transition in figure 3 below which a much more rapid steepening of the spectral slope is observed for the narrowband 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 ∼ P 1/3 shown in figure 4, since bound waves account for energy which does not follow the WTT cascade pathway.

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) (3.13) ω 1 + ω 2 = ω 3 + ω, (3.14) where ω 2 i = 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 (3.13) and (3.14) 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 (3.14) as with ∆ω 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 where k 3 ≡ −k + k 1 + k 2 , · denotes the time average (for stationary states), and * denotes the complex conjugate.η f (k, t) corresponds to the free-wave modes, and is computed by the spatial Fourier transform of η f (x, t) with the latter obtained after applying the filter (3.11) to the wavenumber-frequency spectra. By definition, the function T (k, k 1 , k 2 ) continuously varies between 0 and 1, with 1 corresponding to exact resonance among k, k 1 , k 2 and k 3 . In (3.16), we useη f instead ofη 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 k 1 = (40, 40) and k 2 = (20, −40) so that T (k) with k = (k x , k y ) can be shown by two-dimensional contour plots. Figure 9 shows such plots of T (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 (3.13) and (3.14), and discrete resonant points (red crosses) with 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η f in (3.16)). 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 k 1 and 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 ω = 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 measuresL b and I respectively for the broadening width and overall interaction strength. ForL b , we define it as a characteristic width (Pan & Yue 2017): where |Ω| is the normalized frequency mismatch given by |Ω| = |Ω|/k −1/2 with Ω ≡ ω 1 + ω 2 − ω 3 − ω k and the denominator k −1/2 estimating the frequency discreteness at k associated with the wavenumber spacing. In equation (3.17), the summation is for all the grid points of k. Thus,L b measures the nonlinear broadening around the exact solutions by the first moment of T (k).  figure 10. We can observe a general increase ofL b and I for increasing , 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 L b and I are only possible to resolve by usingη 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 ofL 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 than the other two cases, but not largerL b or I).

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.
(ii) For 0 < k x < 60, we need Ω d = ω 1 + ω 2 − ω 3 − ω = 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.
In summary, the resonant manifold shown in figure 9 corresponds to the second case above, which is the only possibility to have solutions in Ω d = 0 for a wave field traveling to the positive x direction.