2. Numerical model

In this study, we are focusing exclusively on wave turbulence arising from nonlinear interactions of random dispersive surface gravity waves. The numerical model employed here is the fully nonlinear model of Chalikov et al. (Reference Chalikov, Babanin and Sanina2014), specifically developed for long-term simulations of three-dimensional (3-D) fully nonlinear periodic potential gravity waves in deep water. This model is ideal for examining the long-term evolution of multimode gravity wave fields in deep water. Given that nonlinear wave interactions and spectral evolution can occur over tens or even hundreds of wave periods, we need a model that is precise enough to accurately reproduce the relatively slow spectral evolution while allowing for the consideration of various wave conditions with multiple realisations (Chalikov et al. Reference Chalikov, Babanin and Sanina2014).

The model involves a transformation of the 3-D equations for potential flow with periodic boundary conditions to a non-stationary surface-following non-orthogonal coordinate system. The model considers the non-dimensionalised form of the 3-D Laplace (2.1) for velocity potential, along with its corresponding boundary conditions at the free surface:

(2.1)$$\begin{gather} {\phi _{xx}} + {\varUpsilon ^2}{\phi _{yy}} + {\phi _{zz}} = 0, \end{gather}$$

(2.2)$$\begin{gather}\eta _{t}+\eta _{x}\varphi _{x}+\varUpsilon ^2\eta _{y}\varphi _{y}-\phi _{z}=0, \end{gather}$$

(2.3)$$\begin{gather}\varphi _{t}+\tfrac{1}{2}(\varphi _{x}^{2}+\varUpsilon ^2\varphi _{y}^{2}+\phi _{z}^{2})+\eta+p=0. \end{gather}$$

Here, $(x,y,z)$ is the Cartesian system, $t$ is time and ${{\eta }}(x,y,t)$ represents the single-valued interface (free surface), $\phi$ is the 3-D velocity potential, ${{\varphi }}$ is a value of $\phi$ at the surface and ${{p}}$ is the external pressure created by the flow above the surface and normalised using the density of water. Since the model focuses on gravity waves, the surface tension term is omitted (Chalikov et al. Reference Chalikov, Babanin and Sanina2014). Equations (2.1)–(2.3) are given in Cartesian coordinates; however, they are transformed into surface following a coordinate system.

In the model of Chalikov et al. (Reference Chalikov, Babanin and Sanina2014), the velocity potential is expressed as a combination of linear and nonlinear components. The linear component is computed using an analytical solution, whereas the nonlinear aspects are determined iteratively. Typically, the linear component of the velocity potential is larger by two orders of magnitude than the nonlinear component. This approach of separating the components enables the use of larger time steps while minimising the requirement for extensive iterations, thereby improving computational efficiency (Chalikov et al. Reference Chalikov, Babanin and Sanina2014). The following scales are used for the non-dimensional equations of (2.1)–(2.3): length ($L$), time ($L^{1/2}g^{1/2}$) and velocity potential ($L^{3/2}g^{1/2}$). The pressure is normalised by water density to obtain the pressure scale ($Lg$). Here $\varUpsilon$ is the ratio $L/L_y$, where $L$ and $L_y$ are length scales in the $x$ and $y$ directions, respectively. This ratio was introduced since the principal equations are solved in a square domain. This approach allows for faster computations for the considered principal equations (Chalikov et al. Reference Chalikov, Babanin and Sanina2014; Chalikov Reference Chalikov2020).

Two damping mechanisms have been included in the model, namely, conventional viscous damping and wave breaking. The accurate modelling of wave breaking has received particular attention, as it is a crucial process that affects wave statistics and intermittency (Falcon, Roux & Laroche Reference Falcon, Roux and Laroche2010b; Babanin Reference Babanin2013; Slunyaev & Kokorina Reference Slunyaev and Kokorina2019). Since the numerical model assumes the free surface to be a single-valued function, wave breaking cannot be solved explicitly. Instead, wave breaking is handled as a two-step calculation. First, the prediction of the wave breaking onset, and then the calculation of parameterisation of breaking dissipation (energy loss) (Seiffert, Ducrozet & Bonnefoy Reference Seiffert, Ducrozet and Bonnefoy2017; Chalikov Reference Chalikov2020). Prediction of the wave breaking onset refers to predicting the time and location of the start of the breaking process so that numerical instability in models can be prevented. More specifically, the wave breaking onset is defined with the local steepness of the simulated surface reaches the threshold value of 1.1 (Chalikov et al. Reference Chalikov, Babanin and Sanina2014). Once the breaking onset criterion is fulfilled, it triggers the calculation of energy loss due to breaking. This energy loss is parameterised by introducing the conventional diffusion terms defined with the second horizontal derivatives of the surface height in the evolutionary equations of surface elevation and velocity potential (Chalikov et al. Reference Chalikov, Babanin and Sanina2014; Chalikov Reference Chalikov2016, Reference Chalikov2020). The details of the implemented numerical scheme for dissipation due to wave breaking including performance benchmarking can be found in Chalikov (Reference Chalikov2016).

For the initial wave conditions, we have considered the conventional Joint North Sea Wave Project (JONSWAP) spectrum form (Hasselmann et al. Reference Hasselmann1973; Holthuijsen Reference Holthuijsen2007):

(2.4)\begin{equation} S(\omega,\theta)= F(\omega)G(\theta), \end{equation}

where $F(\omega )$ represents the distribution of energy across frequencies, which is subsequently transformed into an energy distribution over wavenumbers. Meanwhile, $G(\theta )$ describes the directional spreading, hence, the anisotropy of the spectrum (see Holthuijsen Reference Holthuijsen2007; Hasselmann et al. Reference Hasselmann1973)

(2.5)\begin{equation} F(\omega)=\frac{\alpha g^2}{\omega^5}\exp{\left[-\frac{5}{4} \left(\frac{\omega}{\omega_{peak}}\right)^{{-}4}\right]} \gamma^{\exp{\left[-\frac{(\omega-\omega_{peak})^2}{2 \sigma^2{\omega_{peak}}^2}\right]}}, \end{equation}

where $\alpha$ is the Phillips parameter that controls the spectrum energy and the significant wave height (see the following), $\gamma$ is the peak enhancement factor, $\omega _{peak}$ is the peak angular frequency and $\sigma$ is the dimensionless ‘width’ of the spectrum that varies with respect to frequency as in

(2.6) \begin{equation} \sigma = \left\{\begin{array}{@{}cc@{}} {\sigma_a} = 0.07, & \omega \le {\omega_{peak}},\\ {\sigma _b} = 0.09, & \omega >{\omega_{peak}}. \end{array} \right\} \end{equation}

The directional spreading function $G(\theta )$ has a form similar to that given in Donelan, Hui & Hamilton (Reference Donelan, Hui and Hamilton1985):

(2.7)\begin{equation} G(\theta)= A_N\cos^N(\theta), \end{equation}

where $\theta$ is the propagation angle and $N$ is the directional spreading exponent. The higher value of $N$ corresponds to the narrower spectrum.

In the simulations, initial wave fields are constructed in the wavenumber domain, translating the initial spectra from their frequency form, $S(\omega,\theta )$, into the non-dimensional wavenumber form, $S(k,\theta )$. One of the main advantages of using the non-dimensional wavenumber spectra in the numerical model is the increased resolution towards higher frequencies, due to the $k = \omega ^2$ relationship derived from the dispersion equation. In other words, to account for $\delta \omega$ towards the higher frequency end of the spectrum, we need to move $\delta k^2$ on the wavenumber axis, where the resolution is constant, as non-dimensional wavenumbers are represented as $k = 1, 2, 3, \ldots, k_{{max}}$. Considering higher resolution at high-frequency modes enhances accuracy and provides a more accurate representation of nonlinear interactions at the spectrum's tail, which is crucial for wave turbulence studies. These nonlinear interactions, involving both bound and free modes, are known to lead to the emergence of coherent structures and deviations from the predictions of WWT. However, the trade-off of this methodology is that it becomes numerically challenging to consider longer tails on the frequency axis, as it requires two orders of magnitude more wave modes. Another limitation of the model is that it does not consider surface tension and, thus, the simulations do not account for capillary effects (Falcon et al. Reference Falcon, Fauve and Laroche2007, Reference Falcon, Roux and Laroche2010b). Small-amplitude capillary bursts are known to drive intermittency in wave turbulence (Deike et al. Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015; Falcon et al. Reference Falcon, Michel, Prabhudesai, Cazaubiel, Berhanu, Mordant, Aumaître and Bonnefoy2020). Nevertheless, the focus of this study is not to investigate the origin of intermittency but to explore its relationship with wave characteristics such as steepness and directionality.

It is conventionally accepted that the most important ‘aggregated’ parameter that controls the nonlinear behaviour of waves is wave steepness (Holthuijsen Reference Holthuijsen2007; Babanin et al. Reference Babanin, Chalikov, Young and Savelyev2010). The formal relation for wave steepness, $\varepsilon$, and wave spectrum is given by the following expressions:

(2.8)\begin{equation} \varepsilon=(H_s/2)|\boldsymbol{k}_{peak}| =2\sqrt{E}|\boldsymbol{k}_{peak}|,\quad E = \mathop \int \nolimits_0^\infty S( {\omega } )\, {\rm d}\omega, \end{equation}

where $H_s$ is the significant wave height and $E$ is the zeroth-order moment of the wave spectrum; for qualitative trend analysis we can use an approximate expression proposed in Onorato et al. (Reference Onorato2009).

3. Numerical results

In this study, our objective is to investigate the effect of wave steepness and wave directionality (as defined in (2.7)) on intermittency through numerical simulations. Previous studies (Fadaeiazar et al. Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018; Skvortsov et al. Reference Skvortsov, Kirezci, Sgarioto and Babanin2022) have demonstrated that nonlinearity induces deviations from Gaussian behaviour in the free surface and also has the potential to lead to wave breaking, thereby introducing intermittency. As a result, parameters influencing the nonlinear properties of the surface also play a role in the intermittency of surface wave turbulence.

It is known that an increase in the directional spreading of wave energy weakens the magnitude of nonlinearity effects in the wave turbulence (Onorato et al. Reference Onorato2009; Waseda, Kinoshita & Tamura Reference Waseda, Kinoshita and Tamura2009), which also leads to changes in intermittent behaviour (Deike et al. Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015; Fadaeiazar et al. Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018). Based on the laboratory experiments, Fadaeiazar et al. (Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018) reported that the narrow directional distribution corresponds to a greater intermittency.

In our simulations, the selection of these parameters was informed by insights from previous research (Fadaeiazar et al. Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018; Alberello et al. Reference Alberello, Onorato, Frascoli and Toffoli2019; Violante-Carvalho et al. Reference Violante-Carvalho, Skvortsov, Babanin, Pereira, Pinho and Esperança2021; Skvortsov et al. Reference Skvortsov, Kirezci, Sgarioto and Babanin2022) and aimed to encompass a diverse range of conditions within the realistic limits of oceanic scenarios. In simulations, 1024 grid points on the $x$ direction and 512 grid points on the $y$ direction were used. The temporal and spatial resolutions of simulations were selected as $0.04$ $\lambda _{peak}$ and $0.005$ $T_{peak}$, respectively. The simulations utilise 256 modes, with the largest wavenumber set to 6.56 times the peak wavenumber. This set-up provides significantly high resolution at the tail, albeit resulting in a shorter tail. Since this study does not extend into the capillary range, the shorter tail is considered acceptable. Although we explored using longer spectral tails, maintaining the same high resolution proved challenging due to numerical instabilities, especially in the highly nonlinear cases that are crucial for this study. The duration of simulations in both models was chosen as $200 T_{{peak}}$ to consider longer-term nonlinear interactions. To ensure robust statistical analysis of the $\zeta _p$ exponents, we conducted 10 realisations for each scenario with initial phase randomisation. This intensive approach to data collection is made feasible by the computational efficiency of our model, allowing us to perform multiple cases with 10 realisations each. This capability sets our chosen model apart, as it enables thorough exploration and reliable statistical conclusions which would be impractical with other models:

(3.1a–c)\begin{equation} \lambda_{peak} =\frac{2{\rm \pi}}{k_{peak}}, \quad T_{peak} = \frac{2{\rm \pi}}{\omega_{peak}}, \quad \omega_{peak}=\sqrt{g {k_{peak}}}, \end{equation}

and $\lambda _{peak}$ and $T_{peak}$ are the peak wavelength and the peak wave period, respectively, and $g$ is the acceleration of gravity. In simulations, we have considered four distinct frequency spectra: the Pierson–Moskowitz (PM) spectrum, which represents fully developed wind–sea states, and the conventional JONSWAP spectrum, with shape parameters $\gamma =3.30$ and $\alpha =0.0081$, used to represent developing sea conditions (table 1). In addition, we have considered two additional spectra labelled as A and B, which exhibit higher steepness conditions (table 1). These four frequency spectra are employed in conjunction with five different directional spreading conditions, where the directional spreading parameter $N$ in (2.7) ranges from 1 to 800 (as detailed in table 1). It should be noted that the spectral shapes and steepness values used here refer to the initial spectra at $t=0$. As the spectra evolve over time, both the steepness values and spectral shapes change. These changes are more pronounced in spectra with larger steepness. The limit of high $N$ was explored to investigate the transition from two-dimensional to unidirectional behaviour of gravity wave turbulence which according to Zakharov (Reference Zakharov1967), Zakharov, Dias & Pushkarev (Reference Zakharov, Dias and Pushkarev2004), Majda, McLaughlin & Tabak (Reference Majda, McLaughlin and Tabak1997) and others may incur a new wave phenomenology. We have also included the unidirectional simulations results presented in Skvortsov et al. (Reference Skvortsov, Kirezci, Sgarioto and Babanin2022) here as a reference to quantify the intermittent behaviour.

Table 1. Values of $\zeta _p/\zeta _2$ at $p=12$.

$^*$Unidirectional results are obtained from Skvortsov et al. (Reference Skvortsov, Kirezci, Sgarioto and Babanin2022) and from additional simulations conducted using the one-dimensional model developed by Chalikov & Sheinin (Reference Chalikov and Sheinin1996).

As a preliminary step, we validated the predictions of WWT for the wave energy density spectrum of surface waves simulated, and we expected that the simulation should saturate to a Kolmogorov-type velocity spectrum. The results of the simulations for various scenarios are presented in figure 1. It is observed that the energy density spectrum $E(\omega )$ did resemble the power-law asymptotes with $E(\omega ) \propto \omega ^{-5}$. As we begin with an initial JONSWAP sea state calculated from wavenumber spectra with random phase and no forcing conditions considered, wave breaking dissipation becomes dominant as the spectra evolves over time, leading to an $\omega ^{-5}$ relation observed in the spectral tail. However, we observed a slight deviation from this $\omega ^{-5}$ relation at higher frequencies. This deviation can be attributed to the flux of energy resulting from nonlinear interactions occurring in higher modes (Chalikov et al. Reference Chalikov, Babanin and Sanina2014). The decay observed at the very end of the tail can be attributed to a damping mechanism similar to viscous damping introduced in other numerical solutions (Chalikov et al. Reference Chalikov, Babanin and Sanina2014).

Figure 1. (a) One-dimensional normalised energy density spectra of the cases with spreading parameter $N=10$. The dashed line is the reference slope of $\omega ^{-5}$. (b) Directional spreading as a function of angle $\theta$.

In line with the frozen turbulence hypothesis of Taylor (Reference Taylor1922), we implemented the framework of Nazarenko (Reference Nazarenko2011), Falcon et al. (Reference Falcon, Roux and Audit2010a), Fadaeiazar et al. (Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018) and Rusaouen et al. (Reference Rusaouen, Chabaud, Salort and Roche2017) and we performed calculations of structure-functions ($S_p$) in the time domain (for a given position in space), namely, $S_p(\tau )$, where $\tau$ is the small time separation (for the rationale of this approach, see Nazarenko Reference Nazarenko2011). In the time domain, $S_p$ and surface elevation increments, $\delta \eta (\tau )$, can be rewritten as

(3.2)\begin{equation} S_p=\langle|\delta \eta(\tau)|^p\rangle, \quad \delta \eta(\tau)=\eta(t+\tau)-\eta(t). \end{equation}

The power exponents $n$ for the wave spectra obtained in the simulations are greater than three, so scaling of $S_2^{(\tau )} \propto \tau ^{\zeta _2}$ with $\zeta _2=n-1$ does not hold (Falcon et al. Reference Falcon, Roux and Audit2010a) and leads to the trivial scaling when the first-order differences of the surface elevation is used as given in (3.2). To remove this trivial scaling from $S_p(\tau )$, we have employed the ESS approach and followed the framework of Falcon et al. (Reference Falcon, Roux and Audit2010a) and applied the higher-degree difference statistics adapted to the exponent of the power spectrum $n$. In our study, we considered

(3.3)\begin{equation} S_2^{(3)} (\tau)=\langle|\delta^{(3)} \eta|^2\rangle, \quad \delta^{(3)}\eta=\eta(t+3\tau)-3\eta(t+2\tau)+3\eta(t+\tau)-\eta(t), \end{equation}

and then used $S_2^{(3)} (\tau )$ for estimation of $\zeta ^{q}_p$ from the log–log slope of expression $S_p(r) \propto [S_2(r)]^{\zeta _p}$. This statistical framework employing higher-order differences (similar to (3.3)) of the surface elevation is universal (namely, not restricted by any particular wave interaction mechanism) and can be applied for resonant (waves) and non-resonant (conventional turbulent flow) models of turbulence.

Figure 2 illustrates the structure-functions of third-order surface elevation differences derived from the JONSWAP spectra with $N=30$ directional spreading. These functions are presented in relation to both non-dimensional time lag and the second-order structure-function. With our specified spectral resolution and tail length, we were able to compute $S_p$ up to the $12$th order. Notably, the structure-functions exhibit a linear increase in a double logarithmic plane up to $\tau / T_{peak} \approx 0.20$ (figure 2a). The structure-functions continue to increase until $\tau / T_{peak} \approx 0.45$, after which they decline, reflecting the periodic nature of water waves, consistent with findings by Fadaeiazar et al. (Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018) and Alberello et al. (Reference Alberello, Onorato, Frascoli and Toffoli2019). We determined the relative scaling exponents using the ESS approach, applying least-squares fit to figure 2(b).

Figure 2. Sample structure-functions of the third-order differences of surface elevation (a) as a function of non-dimensional time lag up to the $12$th order and (b) as a function of the second-order structure-function. The plots are given for the case with JONSWAP spectrum and $N=30$.

In figure 3, the probability density function (p.d.f.) of the third-order differences of surface elevation for different time lags is presented. Deviations from Gaussian behaviour are observed at small time separations. As the time separation approaches $\tau / T_{peak} \approx 0.40\unicode{x2013}0.45$, it is noted that the p.d.f.s tend to converge towards a Gaussian shape, consistent with findings in previous studies (Falcon et al. Reference Falcon, Roux and Audit2010a; Chibbaro et al. Reference Chibbaro, De Lillo and Onorato2017; Fadaeiazar et al. Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018).

Figure 3. Example probability density function (p.d.f.) of the third-order differences of surface elevation for different time lags given for JONSWAP spectrum and $N=30$. Gaussian distribution with zero mean and unit standard deviation (black dashed line).

In figure 3 and table 1, the scaling relations of $\zeta _p/ \zeta _2$ are given as a function of the structure-function order, $p$. The nonlinear behaviour of the function $\zeta _p/ \zeta _2$ is linked to deviations from Gaussian statistics, serving as a signature of intermittency in wave turbulence (Denissenko et al. Reference Denissenko, Lukaschuk and Nazarenko2007; Falcon et al. Reference Falcon, Fauve and Laroche2007; Chibbaro et al. Reference Chibbaro, De Lillo and Onorato2017; Fadaeiazar et al. Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018). Stronger deviations from the straight line, $\zeta _p/\zeta _2=p/2$, correspond to higher intermittency. As shown in figure 3, sea states with a narrower directional spreading exhibit more pronounced deviations from the line of Gaussian statistics, $\zeta _p/\zeta _2=p/2$, indicating stronger intermittency, consistent with previously reported trends (Fadaeiazar et al. Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018).

We have also integrated the outcomes of unidirectional simulations, as detailed in Skvortsov et al. (Reference Skvortsov, Kirezci, Sgarioto and Babanin2022), to serve as a reference for quantifying the intermittent behaviour observed in directional simulations. As expected, unidirectional simulations consistently show higher intermittency compared with their directional counterparts (figure 4). This difference in the intermittent behaviour of unidirectional and directional simulations is strongly related to reducing nonlinearity with the directional distribution of energy. Nevertheless, other factors, including the absence of resonant interaction in unidirectional simulations, wave breaking at the high frequency (short waves), and resolution differences between models also contribute to the different intermittent behaviour. Nonetheless, our findings substantiate the occurrence of heightened intermittency in directionally narrow sea states. However, it is essential to note that exceptions to this trend can arise. For instance, the scenarios with the broadest directionality, $N=4$, display higher intermittency than the scenarios with narrower conditions of $N=10$ and $N=30$ (figure 4).

Figure 4. The effect of directional spreading on intermittency: ratio $\zeta _p/ \zeta _2$ vs the exponent of the structure-function as a function of the order $p$ for the different directional spreading parameter $N$, (2.7) and wave spectra for the selected scenarios (a) PM spectra, (b) JONSWAP spectra ($\gamma =3.3$), (c) spectra A and (d) spectra B. The solid lines correspond to the no-intermittency relation ($\zeta _p/ \zeta _2=p/2$) given in (1.3).

We speculate that the occurrence of this non-monotonic pattern may be related to a possible interlay of the reduced efficiency of the quasi-resonant wave–wave interactions, and emerging mechanism of the so-called directional energy focusing of random waves (Kharif & Pelinovsky Reference Kharif and Pelinovsky2003; Fochesato, Grilli & Dias Reference Fochesato, Grilli and Dias2007; Babanin et al. Reference Babanin, Waseda, Shugan and Hwung2011; Kirezci, Babanin & Chalikov Reference Kirezci, Babanin and Chalikov2021) resulting in a random local increase of wave amplitudes (superposition of waves). This is one of the possible mechanisms of rogue wave formation (Fochesato et al. Reference Fochesato, Grilli and Dias2007; Häfner, Gemmrich & Jochum Reference Häfner, Gemmrich and Jochum2021) and a common cause of wave breaking (Babanin Reference Babanin2011). Moreover, within the late stages of this superposition, a greater contribution of nonlinear wave interaction occurs due to increased local wave steepness. This, in turn, leads to more frequent occurrences of sharper crests and wave breaking hence results in higher intermittency. In addition, the scaling relations presented in (table 1) and in figure 3 also highlights that scenarios with a high steepness spectrum exhibit the most pronounced intermittency. The instances of intermittency are most pronounced especially in cases A and B. This observation aligns with established knowledge that the steepness of waves influences wave breaking, a phenomenon known to contribute to intermittent behaviour (Connaughton, Nazarenko & Newell Reference Connaughton, Nazarenko and Newell2003; Yokoyama Reference Yokoyama2004; Skvortsov et al. Reference Skvortsov, Kirezci, Sgarioto and Babanin2022).

Finally, we compare our intermittency results with those from previously reported studies (Falcon et al. Reference Falcon, Fauve and Laroche2007; Mininni & Pouquet Reference Mininni and Pouquet2009; Falcon et al. Reference Falcon, Roux and Audit2010a,Reference Falcon, Roux and Larocheb; Deike et al. Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015; Chibbaro & Josserand Reference Chibbaro and Josserand2016) (see figure 4) and also present our findings across various universal scaling models (see figure 5). According to Falcon et al. (Reference Falcon, Fauve and Laroche2007, Reference Falcon, Roux and Audit2010a,Reference Falcon, Roux and Larocheb), Deike et al. (Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015) and Mordant et al. (Reference Mordant, Delour, Léveque, Arnéodo and Pinton2002), the dependency of $\zeta _p$ can be described by a quadratic polynomial model

(3.4)\begin{equation} \zeta_p=c_1p - (c_2/2)p^2, \end{equation}

where factors $c_1$ and $c_2$ are shown to be forcing dependent and are in the range $1.5< c_1<3.0$, $0.1< c_2<0.5$ with respect to applied conditions in Falcon et al. (Reference Falcon, Fauve and Laroche2007, Reference Falcon, Roux and Audit2010a,Reference Falcon, Roux and Larocheb) and Deike et al. (Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015). The formal justification of this approximation comes from the seminal Kolmogorov model that incorporates the fluctuations of the dissipation energy (Kolmogorov Reference Kolmogorov1962). In our case, these fluctuations can be attributed to random wave breaking. The second term in (3.4) is the nonlinear correction to the linear scaling predicted by the Gaussian statistics and $c_2\neq 0$ indicates the occurrence of intermittency (Falcon et al. Reference Falcon, Fauve and Laroche2007; Deike et al. Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015).

Figure 5. Parabolic fits to the simulation data for the low-order range ($p \le 6$) structure-functions, see (3.4): 1 is the line of no-intermittency relation $\zeta _p/\zeta _2=p/2$; 2 is the shaded area of parabolic fit with the value and variations of coefficients $c_1$ and $c_2$ taken from Falcon et al. (Reference Falcon, Roux and Laroche2010b); 3 and 4 are the scatter points of averaged values of $\zeta _p/\zeta _2$ with the error bars for the unidirectional and directional simulations, respectively; 5 and 6 are respective parabolic model fits for the results of unidirectional and directional simulations, respectively; 7 and 8 are the experimental observation values taken from Fadaeiazar et al. (Reference Fadaeiazar, Alberello, Onorato, Leontini, Frascoli, Waseda and Toffoli2018) for JONSWAP ($\gamma =3.0$) spectra with $N=10$ and $N=1000$ directional spreading parameters, respectively.

Here, the parameters $c_1$ and $c_2$ are calculated using averaged $\zeta _p$ values. For the high-order range ($p \le 12$), we have estimated $c_1 = 2.38$ and $c_2 = 0.006$. If we restrict our consideration to the range $p \le 6$ (as reported in Falcon et al. Reference Falcon, Fauve and Laroche2007, Reference Falcon, Roux and Laroche2010b), then our estimation of $c_2$ aligns even more closely with the values given in Falcon et al. (Reference Falcon, Fauve and Laroche2007, Reference Falcon, Roux and Laroche2010b) and Deike et al. (Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015). Specifically, when $p \le 6$, $c_1 = 2.74$ and $c_2 = 0.11$, respectively. In addition, parameters $c_1$ and $c_2$ are related to the slope of the wave energy spectrum ($n$ in (1.4)) by $n = \zeta _2 + 1 = 2(c_1 - c_2) + 1$ (Falcon et al. Reference Falcon, Fauve and Laroche2007; Deike et al. Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015). Our estimation of $n$ for wave energy spectra falls within the range of $[-5,-7]$. Consequently, the values of $c_1$ and $c_2$ for both the low- and high-order range application satisfy the relation $n=\zeta _2+1$.

Translation of (3.4) to the ratio $\zeta _p/\zeta _2$ leads to

(3.5)\begin{equation} \zeta_p/\zeta_2 =C_1 p/2 - C_2(p/2)^2,\end{equation}

where $C_1 = c_1 (c_1 - c_2)$ and $C_2 = c_2(c_1 - c_2)$. Scaling $\zeta _p/\zeta _2 = p/2$, predicted by Kolmogorov–Zakharov scaling, implies that $C_1 = 1$ in (3.5). For $C_1 = 1$, we calculated $C_2=0.015$. A comparison of our numerical results and the results of Falcon et al. (Reference Falcon, Fauve and Laroche2007, Reference Falcon, Roux and Laroche2010b), and the fit given by (3.5) are depicted in figure 5. Within the significant uncertainties in estimation of parameters $C_1$ and $C_2$ we observed good agreement.

It is evident that (3.4) cannot be extended for $p \rightarrow \infty$ (where this expression becomes negative), and other models should be considered (She & Leveque Reference She and Leveque1994; Jou Reference Jou1997; St-Jean Reference St-Jean2005; Yakhot Reference Yakhot2006; Mininni & Pouquet Reference Mininni and Pouquet2009; Chibbaro & Josserand Reference Chibbaro and Josserand2016). We begin with the conventional She–Leveque (SL) model (also known as the hierarchical structure model) (She & Leveque Reference She and Leveque1994). The original SL model is formulated as a correction for Kolmogorov scaling of hydrodynamic turbulence $\zeta _p/\zeta _3=p/3$. Here, we used the one-parameter version of the SL model (Meyrand, Kiyani & Galtier Reference Meyrand, Kiyani and Galtier2015) which has been adapted for the scaling relations for surface wave turbulence. We follow the framework of Meyrand et al. (Reference Meyrand, Kiyani and Galtier2015) and use the following ansatz for $\zeta _p/\zeta _2$:

(3.6)\begin{equation} \zeta_p/\zeta_2 = C_1\frac{p}{2} + C_{0} - C_{0} \left ( 1- \frac{2}{3C_{0} } \right)^{p/2}, \end{equation}

where $C_1, C_{0}= const$. By imposing condition $\zeta _2 = 1$, we obtain that $C_1=1/3$ and $C_{0}$ values are evaluated from data ($C_{0} =11.37$).

In a similar way, we evaluated the fit of the multifractal (MF) model of intermittency (Benzi et al. Reference Benzi, Paladin, Parisi and Vulpiani1984). The MF model was initially proposed by Benzi et al. (Reference Benzi, Paladin, Parisi and Vulpiani1984) and can be reduced to

(3.7)\begin{equation} \zeta_p/\zeta_2 = \frac{p}{2} - \frac{1}{2}{\log_{2}}\left[C_{mf} + (1-C_{mf}) \left(\frac{1}{2}\right)^{1-p/2} \right], \end{equation}

where $C_{mf}$ is the ‘free’ parameter evaluated from data, calculated as $C_{mf} =0.92$. For high values of $p$, our numerical results are plotted in figure 6 which also presents the results for SF (3.6) and MF (3.7) models. We observe that both models provide a reasonable fit for the intermediate range of $p$, $p <8$, but for higher $p$ the SL model performs better.

Figure 6. She–Leveque (SL) and multifractal (MF) model fits to the simulation data for the high-order range ($p \le 12$) structure-functions, see ((3.6) and (3.7)): 1 is the line of no-intermittency relation $\zeta _p/\zeta _2=p/2$; 2 and 3 are the shaded areas for unidirectional and directional simulations, respectively; 4 and 5 are the lines for the averaged $\zeta _p/\zeta _2$ values for unidirectional and directional simulations, respectively; 6 and 7 are the SL fits for the results of unidirectional and directional simulations, respectively; 8 and 9 are the MF fits for the results of unidirectional and directional simulations, respectively.

Comparing our findings with other experimental studies such as Falcon et al. (Reference Falcon, Fauve and Laroche2007) and Deike et al. (Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015), we observe a lower level of intermittency. This discrepancy may arise from several factors, including the exclusion of capillary range effects in our study. It is important to recognise that gravity dominates in oceanic spectra, whereas capillary effects are more pronounced in laboratory settings due to the limited size of wave basins. Therefore, these differences can be expected. Furthermore, it is known that inverse cascades occur in gravity wave scales but not in pure capillary waves, which likely influences the intermittency characteristics observed in wave turbulence studies.

Another simulation set was conducted to consider a longer tail case ($k_{{max}} / k_{p} = 16$). The long tail set employs reduced spectral resolution to prevent numerical instabilities, focusing exclusively on JONSWAP and PM spectra, as spectra A and B encounter numerical instability with this configuration. Even though the results of that simulation set are not given here, we observed a higher level of intermittency similar to what was seen in a one-dimensional simulation dataset that also considered a longer tail. Including higher frequencies, as expected, enhanced the formation of coherent structures and resulted in increased anomalous scaling exponents, indicative of intermittent behaviour. However, the structure-functions exhibited a linear trend over a limited range of time lags up to the seventh order. Attempts to extend beyond this order led to spurious oscillations, likely due to low resolution. Furthermore, we did not observe the anticipated decline in structure-functions after $\tau / T_{peak} \approx 0.5$, which typically signifies the periodic nature of water waves.