Intermittency of gravity wave turbulence on the surface of an infinitely deep fluid: directional effects

Abstract

This study investigates the influence of surface wave characteristics, specifically wave steepness and directional spreading, on intermittency in deep-water gravity wave turbulence through long-term numerical simulations of three-dimensional potential fully nonlinear periodic gravity waves. We conducted this investigation by estimating the scaling exponent of the surface elevation under different sea state conditions. With our numerical methods, we were able to evaluate the scaling exponents of the structure-function up to 12th order. The observed increased intermittency in directionally narrower sea states and in higher steepness conditions aligns with known effects of quasi-resonant wave–wave interactions and wave breaking. Comparative analyses reveal that both the conventional She–Leveque model and the multifractal models, also used to represent intermittency in wave turbulence of a different nature, exhibit a strong correlation in this study. This observation underscores the universality of intermittency phenomena within wave turbulence.

1. Introduction

The nonlinear interactions of random dispersive surface gravity waves lead to an energy cascade similar to that observed in hydrodynamic turbulence when out of equilibrium (Mordant & Miquel 2017; Fadaeiazar et al. 2018). This phenomenon is known as wave turbulence. Unlike hydrodynamic turbulence, where the motion is a hydrodynamic vortex, wave turbulence involves the propagation of waves (Nazarenko 2011). The wave turbulence is also free from the closure problem encountered in hydrodynamic turbulence (Galtier 2024). Here, we are investigating a specific phenomenon: intermittency in gravity wave turbulence exhibited by nonlinear gravity waves. This phenomenon is commonly depicted through the utilisation of statistical measures applied to surface increments (free-surface elevation and velocity). Although several physical processes are known to affect intermittent behaviour in wave turbulence, their effects are yet to be quantified accurately. This study deals with the intermittent characteristics of directional nonlinear surface conditions.

The general theory of intermittency in turbulence remains an open question and an active area of research. There is no universal method to derive intermittency statistics from the underlying equations of motion of a dynamic system (Falcon, Fauve & Laroche 2007). Various phenomenological models employ different physical assumptions, including fluctuations of dissipation energy, wave breaking, bifurcation and others (Frisch 1995; Biven, Nazarenko & Newell 2001; Newell, Nazarenko & Biven 2001; Nazarenko 2011). Conventionally, the signature of intermittency is associated with the notable departure from the predictions of the Gaussian statistics for random field variables or the emergence of the large deviation statistics for the turbulence bursts (Frisch 1995). Mathematically, intermittency in turbulence can be defined as the ‘anomalous’ scaling of the structural function (Frisch 1995)

(1.1)\begin{equation} S_p(r) \propto r^{\zeta_p}, \end{equation}

where $S_p (r) = \langle |\delta \boldsymbol {v} |^p \rangle$, $\delta \boldsymbol {v} = \boldsymbol {v} (\boldsymbol {r}_1, t) - \boldsymbol {v}(\boldsymbol {r}_2, t)$, $r = |\boldsymbol {r}_1 - \boldsymbol {r}_2|$ and $\boldsymbol {v}(\boldsymbol {r}, t)$ is the random velocity field of a turbulent motion (or any other field variable). With this notation, the intermittency simply corresponds to the nonlinearity of function $\zeta _p \equiv \zeta (p)$ (Frisch 1995). For Kolmogorov theory of isotropic turbulence $S_3 \propto r$ and this leads to $\zeta _p = p/3$ (Frisch 1995); for wave turbulence there is a less-definitive conjecture for $\zeta _2$, see the following.

In this study, we utilise an ensemble of random nonlinear gravity waves on the surface of an infinitely deep fluid and numerically simulate its dynamics. This choice is motivated by two factors. First, the theory of weak wave turbulence (WWT) provides numerous analytical predictions, which can serve as a benchmark for comparison. Second, in this scenario, the dynamics of the field and, consequently, the $\zeta _p$ function can be assessed directly from the ‘first principles’, e.g. by numerical solution of the underlying equations of motion, without any additional closure assumptions (e.g. phase randomisation) and reduction to the kinetic equations as in the conventional framework of the WWT (Chalikov, Babanin & Sanina 2014; Chalikov 2016). From this perspective, the ‘first principles’ approach provides an important insight to the ‘universality’ of the functional form of $\zeta _p$.

As the order $p$ increases, both numerical and experimental uncertainties escalate, posing substantial challenges in the direct validation of the scaling law (1.1) and the robust estimation of the exponents $\zeta (p)$ for wave turbulence, see Benzi et al. (1993), She & Leveque (1994), Fadaeiazar et al. (2018) and others. To overcome these difficulties, the extended self-similarity (ESS) framework is usually employed (Benzi et al. 1993; Falcon et al. 2007; Deike et al. 2015; Chibbaro & Josserand 2016; Fadaeiazar et al. 2018; Alberello et al. 2019). The ESS implies that all exponents $\zeta _p$ can be defined relatively to a particular structure-function $S_p(r)$ for some $p = p_*$ for which $\zeta _{p_*}$ is known. For Kolmogorov turbulence, we have the exact result for $p =3$, $S_3 \propto r$, so $\zeta _3 = 1$ and $S_p(r) \propto [S_3(r)]^{\zeta _p}$ with $\zeta _p = (p/3) \zeta _3$ when intermittency is absent (Benzi et al. 1993; She & Leveque 1994; Frisch 1995). The relation $S_3 \propto r$ does not hold for wave turbulence, and $S_2$ (e.g. $p=2$) is used as a reference structure-function in this study:

(1.2)\begin{equation} S_p(r) \propto [S_2(r)]^{\zeta_p}. \end{equation}

In line with the above comments, the following scaling law for the non-intermittent regime (Gaussian statistics) can be obtained (Newell et al. 2001; Falcon, Roux & Audit 2010a; Chibbaro, De Lillo & Onorato 2017; Alberello et al. 2019):

(1.3)\begin{equation} \zeta_p = (p/2) \zeta_2. \end{equation}

The use of $S_2(r)$ as a reference structure-function is motivated by its relation to the slope $n$ of the energy spectra of wave turbulence that can be estimated by other means (Frisch 1995). More specifically, if the energy spectra is deduced in a power-law form

(1.4)\begin{equation} E(\omega) \propto \omega^{- n}, \end{equation}

where $\omega$ is the angular frequency of field harmonics, the exponent $n$ may provide a conjecture for $\zeta _2 =p/2$ (see Frisch 1995; Chibbaro et al. 2017). In this regard, the choice of $S_2(r)$ as a reference structure-function also helps mitigate the potential non-uniformity of scaling $S_2(r) \propto r^{\zeta _2}$ (and, consequently, higher order structure-functions), which manifests itself in the well-known discontinuity of the slope $n$ in (1.4) of the energy spectra (Frisch 1995; Fadaeiazar et al. 2018; Alberello et al. 2019). The exponent $n = 4$ was recovered in the theoretical studies of Zakharov et al. (2019), Zakharov (1967) and Zakharov & Filonenko (1967), and $n = 5$ is the so-called Phillips spectrum (Phillips 1977; Zakharov et al. 2019). It is assumed that the Phillips spectrum dominates at approximately $\omega \ge 5 \omega _p$, and $\omega _p$ is the circular frequency of the spectral peak (Zakharov et al. 2019). It is worth noting that the model with $n= 6$ has also been reported previously (Denissenko, Lukaschuk & Nazarenko 2007).

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. (2014), 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. 2014).

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. 2014). 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. (2014), 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. 2014). 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. 2014; Chalikov 2020).

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 2010b; Babanin 2013; Slunyaev & Kokorina 2019). 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 2017; Chalikov 2020). 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. 2014). 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. 2014; Chalikov 2016, 2020). The details of the implemented numerical scheme for dissipation due to wave breaking including performance benchmarking can be found in Chalikov (2016).

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

(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 2007; Hasselmann et al. 1973)

(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 (1985):

(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. 2007, 2010b). Small-amplitude capillary bursts are known to drive intermittency in wave turbulence (Deike et al. 2015; Falcon et al. 2020). 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 2007; Babanin et al. 2010). 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. (2009).

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. 2018; Skvortsov et al. 2022) 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. 2009; Waseda, Kinoshita & Tamura 2009), which also leads to changes in intermittent behaviour (Deike et al. 2015; Fadaeiazar et al. 2018). Based on the laboratory experiments, Fadaeiazar et al. (2018) 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. 2018; Alberello et al. 2019; Violante-Carvalho et al. 2021; Skvortsov et al. 2022) 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 (1967), Zakharov, Dias & Pushkarev (2004), Majda, McLaughlin & Tabak (1997) and others may incur a new wave phenomenology. We have also included the unidirectional simulations results presented in Skvortsov et al. (2022) 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. (2022) and from additional simulations conducted using the one-dimensional model developed by Chalikov & Sheinin (1996).

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. 2014). 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. 2014).

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 (1922), we implemented the framework of Nazarenko (2011), Falcon et al. (2010a), Fadaeiazar et al. (2018) and Rusaouen et al. (2017) 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 2011). 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. 2010a) 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. (2010a) 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. (2018) and Alberello et al. (2019). 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. 2010a; Chibbaro et al. 2017; Fadaeiazar et al. 2018).

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. 2007; Falcon et al. 2007; Chibbaro et al. 2017; Fadaeiazar et al. 2018). 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. 2018).

We have also integrated the outcomes of unidirectional simulations, as detailed in Skvortsov et al. (2022), 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 2003; Fochesato, Grilli & Dias 2007; Babanin et al. 2011; Kirezci, Babanin & Chalikov 2021) 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. 2007; Häfner, Gemmrich & Jochum 2021) and a common cause of wave breaking (Babanin 2011). 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 2003; Yokoyama 2004; Skvortsov et al. 2022).

Finally, we compare our intermittency results with those from previously reported studies (Falcon et al. 2007; Mininni & Pouquet 2009; Falcon et al. 2010a,b; Deike et al. 2015; Chibbaro & Josserand 2016) (see figure 4) and also present our findings across various universal scaling models (see figure 5). According to Falcon et al. (2007, 2010a,b), Deike et al. (2015) and Mordant et al. (2002), 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. (2007, 2010a,b) and Deike et al. (2015). The formal justification of this approximation comes from the seminal Kolmogorov model that incorporates the fluctuations of the dissipation energy (Kolmogorov 1962). 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. 2007; Deike et al. 2015).

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. (2010b); 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. (2018) 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. 2007, 2010b), then our estimation of $c_2$ aligns even more closely with the values given in Falcon et al. (2007, 2010b) and Deike et al. (2015). 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. 2007; Deike et al. 2015). 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. (2007, 2010b), 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 1994; Jou 1997; St-Jean 2005; Yakhot 2006; Mininni & Pouquet 2009; Chibbaro & Josserand 2016). We begin with the conventional She–Leveque (SL) model (also known as the hierarchical structure model) (She & Leveque 1994). 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 2015) which has been adapted for the scaling relations for surface wave turbulence. We follow the framework of Meyrand et al. (2015) 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. 1984). The MF model was initially proposed by Benzi et al. (1984) 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. (2007) and Deike et al. (2015), 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.

4. Conclusions

We have examined the intermittency of turbulence in deep-water surface gravity waves. The scaling exponent of the surface elevation structure-function, a crucial parameter characterising intermittency in this context, was assessed through numerical solutions of the velocity potential of the surface. This assessment was conducted without additional closure assumptions or simplifications to kinetic equations. Our investigation has been focused on the effect of various wave characteristics, such as wave steepness and directional spreading of wave energy, on the intermittency of gravity wave turbulence. The observed increase in intermittency in directionally narrower sea states and higher steepness is consistent with previous research trends, indicating the influence of quasi-resonant wave–wave interactions and wave breaking.

Utilising our high-resolution numerical method for modelling nonlinear free surfaces, we have successfully determined scaling exponents for the surface elevation structure-function up to the 12th order. Our findings indicate that, at high orders, the conventional SL model demonstrated a better fit compared with the MF model. However, it is noteworthy that the MF model also performed well in capturing the intermittency behaviour in our numerical data. Comparative analyses were performed, aligning our findings with other analytical and numerical studies on wave turbulence, encompassing different natures and mechanisms of excitation, such as turbulence in bending waves on a plate and magnetohydrodynamic turbulence. These outcomes substantiate the notion of universality in intermittency phenomena within wave turbulence.