1. Introduction
As a result of cheaper computational storage and improved sensors, the number of surface waves included in databases of field measurements has soared over recent decades, going from fifty thousand at the end of the 1970s to hundreds of millions in 2020 (Forristall Reference Forristall1978; Karmpadakis, Swan & Christou Reference Karmpadakis, Swan and Christou2020). They allow for systematic correlation studies with hindcast data, evidencing, for instance, that the probability of occurrence of rogue waves (RWs) is independent of the instantaneous wind speed and direction (Christou & Ewans Reference Christou and Ewans2014). These approaches are undoubtedly valuable as they single out the environmental conditions that favour the occurrence of RW but remain far from being exhaustive. For instance, the overwhelming majority of deep-water waves discussed in this context in Christou & Ewans (Reference Christou and Ewans2014) share the same directions of swell and current, precluding the possible evidence of generation of RW by wave–current interactions, a phenomenon yet recognized as a promising outlook (Adcock & Taylor Reference Adcock and Taylor2014; Toffoli et al. Reference Toffoli, Waseda, Houtani, Calaveri, Greaves and Onorato2015; Ducrozet et al. Reference Ducrozet, Abdolahpour, Nelli and Toffoli2021). More fundamentally, drawing a comprehensive theory of RWs based on these results is complicated by the lack of statistically stationary states: in practice, wave elevation time series from different storms are spliced into 20 min samples then recombined with others sharing similar proxies (e.g. wave mean frequency, mean direction of propagation, etc.), which unavoidably introduces a bias and explains why the distribution of rare events such as RWs is still discussed.
To assess theoretical models, laboratory experiments nicely complement field experiments since they provide long-time statistics under controlled conditions. Most of them take place in long flumes in which unidirectional waves, also referred to as ‘long-crested waves’, are randomly generated by a wave maker and propagate over more than a hundred metres before being damped by a beach. Such experiments typically report a transient overshoot of the kurtosis, of the spectral width and of the RW probability associated with the emergence of high-amplitude structures locally akin to the so-called Peregrine soliton (PS) (Onorato et al. Reference Onorato, Osborne, Serio, Cavaleri, Brandini and Stansberg2004, Reference Onorato, Osborne, Serio and Cavaleri2005, Reference Onorato, Osborne, Serio, Cavaleri, Brandini and Stansberg2006; Shemer & Sergeeva Reference Shemer and Sergeeva2009; Shemer, Sergeeva & Slunyaev Reference Shemer, Sergeeva and Slunyaev2010b; Shemer, Sergeeva & Liberzon Reference Shemer, Sergeeva and Liberzon2010a; Cazaubiel et al. Reference Cazaubiel, Michel, Lepot, Semin, Aumaître, Berhanu, Bonnefoy and Falcon2018; Dematteis et al. Reference Dematteis, Grafke, Onorato and Vanden-Eijnden2019; Michel et al. Reference Michel, Bonnefoy, Ducrozet, Prabhudesai, Cazaubiel, Copie, Tikan, Suret, Randoux and Falcon2020). This dynamics can be modelled by the nonlinear Schrödinger equation (NLSE), an exact solution of the latter, localized in both space and time, being the PS. The instability of a continuous wave train, called the ‘modulation instability’ and generating RWs, can also be studied in long one-dimensional flumes and described by the NLSE, see, e.g. Lighthill (Reference Lighthill1965), Benjamin & Feir (Reference Benjamin and Feir1967), Lake et al. (Reference Lake, Yuen, Rungaldier and Ferguson1977), Melville (Reference Melville1982), Chabchoub et al. (Reference Chabchoub, Genty, Dudley, Kibler and Waseda2017) and references therein. All these results strongly depend on the directionality of the wave field, as shown both theoretically through the existence of transverse instabilities (Badulin & Ivonin Reference Badulin and Ivonin2012; Ablowitz & Cole Reference Ablowitz and Cole2021), numerically (Onorato, Osborne & Serio Reference Onorato, Osborne and Serio2002; Soquet-Juglard et al. Reference Soquet-Juglard, Dysthe, Trulsen, Krogstad and Liu2005; Gramstad & Trulsen Reference Gramstad and Trulsen2007; Toffoli et al. Reference Toffoli, Bitner-Gregersen, Onorato and Babanin2008) and experimentally (Waseda Reference Waseda2006; Onorato et al. Reference Onorato2009), questioning their relevance in accounting for in situ RWs.
On the other hand, another set of experiments investigates the theory of weak wave turbulence (WWT), which predicts how energy spreads among random waves in nonlinear interaction (Falcon & Mordant Reference Falcon and Mordant2022). They take place in basins with reflecting walls and deal with isotropic or at least strongly multidirectional waves (‘short-crested waves’). Until recently, they essentially consisted of generating waves with a wavelength a fraction of the length of the basin and measuring the energy cascade toward small scales (Denissenko, Lukaschuk & Nazarenko Reference Denissenko, Lukaschuk and Nazarenko2007; Lukaschuk et al. Reference Lukaschuk, Nazarenko, McLelland and Denissenko2009; Nazarenko et al. Reference Nazarenko, Lukaschuk, McLelland and Denissenko2010; Deike et al. Reference Deike, Miquel, Gutiérrez, Jamin, Semin, Berhanu, Falcon and Bonnefoy2015; Aubourg et al. Reference Aubourg, Campagne, Peureux, Ardhuin, Sommeria, Viboud and Mordant2017; Campagne et al. Reference Campagne, Hassaini, Redor, Sommeria, Valran, Viboud and Mordant2018). A breakthrough occurred in 2020, when it was evidenced that forcing multidirectional random waves of short wavelengths in such basins not only generates even shorter wavelengths but also larger ones, corresponding to the inverse cascade of WWT (Falcon et al. Reference Falcon, Michel, Prabhudesai, Cazaubiel, Berhanu, Mordant, Aumaître and Bonnefoy2020). Such wave fields are valuable for the study of RWs since the waves involved in their dynamics are spontaneously generated by nonlinear interactions rather than directly forced by the wave maker. Moreover, they verify isotropy, homogeneity and stationarity, and therefore offer a unique framework to confront theoretical predictions on RWs to a simplified though strongly nonlinear model of the sea state. The present study reports the statistics of thousands of RWs measured in such a state and investigates the effect of high-order nonlinearities.
2. Experimental set-up
Experiments are carried out in the large-scale basin (40 m long 
$\times$ 30 m wide 
$\times$ 5 m deep) of Ecole Centrale de Nantes, France. At one end of the basin, 48 flaps of width 
$\ell = 0.62\ \mathrm {m}$ are driven independently by different realizations of white noise filtered in the 
$[f_0-\Delta f, f_0+\Delta f]$ frequency range, with 
$f_0=1.8 \mathrm {Hz}$ the central frequency and 
${\Delta f = 0.2\ \mathrm{Hz}}$ the bandwidth. Therefore, each flap generates independent waves of frequency around 
$f_0$ (wavelength 
$\lambda _0 = 0.48\ \mathrm {m}$, group velocity 
$v_g=0.43\ \mathrm {m}\ \mathrm {s}^{-1}$) with a directional spread that can be estimated as 
$\theta = 2\times (\lambda _0/\ell ) = 88^{\circ }$. Three forcing amplitudes are considered, hereafter referred to, in increasing order, as Runs 1 to 3. At the other end, a solid vertical wall is built ahead of the usual beach. This set-up is sketched in figure 1(a).
(a) Experimental set-up showing the 48 flap wave generator, the end wall and the 23 probes. Probes 1 and 2 are used to verify the wave maker behaviour and are not included in the data analysis. (b) Photograph of a typical wave field (Run 3, the horizontal field of view is approximately one metre).
As reported in Falcon et al. (Reference Falcon, Michel, Prabhudesai, Cazaubiel, Berhanu, Mordant, Aumaître and Bonnefoy2020), a statistically stationary, homogeneous and isotropic nonlinear steady state is reached after a transient of up to twenty minutes. The general picture is as follows: during this transient, the waves generated at 
$f_0$ by the flaps travel nearly 70 times the length of the basin (
$20\ \mathrm {min}/v_g = 2.8\ \mathrm {km}$). As they propagate, nonlinear effects such as four-wave resonant interactions and very steep structures spread energy in all directions. Some of these strongly nonlinear effects visible from the shore are found to occur homogeneously in the basin, e.g. capillary waves generated by large gravity waves. Note that white capping is not observed, see figure 1(b).
The surface elevations 
$\lbrace \eta _i(t)\rbrace _{i=1 \dots 23}$ are recorded by 23 resistive probes of vertical resolution 0.1 mm and frequency resolution 100 Hz during 27 to 30 h depending on the run. These measurements can be used to verify the claims of stationarity, homogeneity and isotropy. Stationarity is confirmed through the time evolution of statistical measurements of the wave field, e.g. the standard deviation of surface elevations computed over one minute samples, and is achieved after up to twenty minutes, see figures in Falcon et al. (Reference Falcon, Michel, Prabhudesai, Cazaubiel, Berhanu, Mordant, Aumaître and Bonnefoy2020). The transients are not investigated in this study and only measurements performed in the steady-state regimes are hereafter discussed. All probes are found to measure a similar standard deviation of surface elevation up to a maximum relative difference of 
$10\,\%$: homogeneity is closely achieved, and to remove the small remaining bias each signal is normalized by the standard deviation of the corresponding probe. Isotropy is the most challenging assumption to test since it cannot be investigated from a single elevation signal. The cross-correlation between pairs of elevation signals is therefore introduced. For each run, it is computed as
\begin{equation} R_{i,j}(\tau) = \frac{\langle \eta_i(t) \eta_j(t+\tau) \rangle}{\sqrt{\langle \eta_i(t)^{2} \rangle \langle \eta_j(t)^{2} \rangle}}, \end{equation}
where 
$\langle {\cdot } \rangle$ denotes a temporal averaging. Over all runs, all lags 
$\tau$ and all probes 
$i \neq j$, 
$\vert R_{i,j} \vert$ remains less than 0.2 and the probes are therefore largely uncorrelated, as expected from their large spatial separation. Nevertheless, the remaining correlations evidence that 
$R_{i,j}(\tau )$ is almost symmetric, i.e. that the wave field is essentially isotropic (
$R_{i,j}(\vert \tau \vert )$ and 
$R_{i,j}(-\vert \tau \vert )$, respectively, account for signals propagating from 
$i$ to 
$j$ and from 
$j$ to 
$i$). Quantitatively, with 
$(i \neq j) \in [13,14,15]$ standing for the three close central probes and 
$\tau _{max}$ such that 
$R_{i,j}(\tau _{max})$ is maximum,
\begin{equation} \left \vert \frac{R_{i,j}(\tau_{max}) - R_{i,j}(-\tau_{max})}{R_{i,j}(\tau_{max})} \right \vert < 0.16, \end{equation}
a strong indication toward isotropy. Finally, note that the power spectrum density 
$S_\eta (f)$ (PSD), reported in figure 2, reveals that most of the energy is located at frequencies smaller than the forcing range, corresponding to waves forced by nonlinear interactions. These PSDs present some features theoretically predicted for the inverse cascade of WWT and derived under the assumption of stationarity, homogeneity and isotropy (Falcon et al. Reference Falcon, Michel, Prabhudesai, Cazaubiel, Berhanu, Mordant, Aumaître and Bonnefoy2020).
PSD of surface elevation for the three different steady states considered. Thin lines correspond to experimental results and thick ones to their numerical model. The forcing bandwidth is also displayed.
It is instructive to detail the energy budget of this wave field. Energy is injected in the wave system at a rate 
$\mathcal {P}_{inj}$ that can be measured through decay experiments and is of several watts (see Falcon et al. (Reference Falcon, Michel, Prabhudesai, Cazaubiel, Berhanu, Mordant, Aumaître and Bonnefoy2020), note that this power is much smaller than the one supplied to the wave maker). Conversely, the power dissipated by viscosity at high frequency (>2 Hz) at the surface boundary layer can be estimated from the experimental PSD and reads (Miles Reference Miles1967)
\begin{equation} \mathcal{P}_{diss} = 2 S \rho g \int_{2\ \mathrm{Hz}}^{\infty} S_\eta(f) \alpha(f) \,\mathrm{d}f, \end{equation}
with 
$S= 30\times 40\ \mathrm {m}^{2}$ the surface of the basin, 
$\rho =10^{3}\ \mathrm {kg}\ \mathrm {m}^{-3}$ the density, 
$g=9.81\ \mathrm {m}\ \mathrm {s}^{-2}$ the acceleration due to gravity, 
$\alpha (f) = 2\nu k^{2} = 2 \nu (2{\rm \pi} f)^{4}/g^{2}$ the damping coefficient for clean water and 
$\nu = 10^{-6}\ \mathrm {m}^{2}\ \mathrm {s}^{-1}$ the kinematic viscosity. We find typically 
$\mathcal {P}_{diss} \sim \mathcal {P}_{inj}/10$, meaning that most of the energy is dissipated by another mechanism than viscous dissipation of high-frequency waves in the bulk. We believe that this mechanism is linked with the nonlinear dynamics at large scales, which involves very steep structures acting as localized sources of dissipation (e.g. cusps of very steep slope).
3. Numerical model
To identify high-order nonlinear effects in the experimental data, these wave fields are reproduced numerically up to second-order nonlinearities. The elevation at a given location is computed as 
$\eta (t) = \eta ^{(1)}+ \eta ^{(2)}$, where the linear contribution 
$\eta ^{(1)}$ is the sum of 
$N_\omega \times N_\theta =512$ independent progressive waves (
$N_\omega =16$ angular frequencies, each of them associated with 
$N_\theta =32$ directions), and 
$\eta ^{(2)}$ is the nonlinear correction. More precisely, 
$\eta ^{(1)}$ reads
\begin{equation} \eta^{(1)}(t) = \sum_{n_\omega=1}^{N_\omega} \sum_{n_\theta=1}^{N_\theta} a_{n_\omega,n_\theta} \cos \left( - \omega_{n_\omega} t + \phi_{n_\omega,n_\theta}\right),\end{equation}
where 
$a_{n_\omega,n_\theta }$ are random numbers drawn from normal distributions of zero mean and standard deviations 
$A_{n_\omega }$. The phase constants 
$\phi _{n_\omega,n_\theta }$ are uniformly distributed in the range 
$[0, 2{\rm \pi} ]$. The leading-order nonlinear correction 
$\eta ^{(2)}$ stems from Longuet-Higgins (Reference Longuet-Higgins1977) (up to a correction factor of one half, see Srokosz Reference Srokosz1986). In particular, it involves the wave vectors of the linear waves, set to model an isotropic wave field as
\begin{equation} \boldsymbol{k}_{n_\omega,n_\theta} = \frac{\omega_{n_\omega}^{2}}{g} \left[\cos \left(\frac{2{\rm \pi} n_\theta}{N_{\theta}}\right) \boldsymbol{e}_x+ \sin\left(\frac{2{\rm \pi} n_\theta}{N_{\theta}}\right) \boldsymbol{e}_y \right]. \end{equation}
The angular frequencies 
$\lbrace \omega _{n_\omega } \rbrace$ are linearly distributed in a given range with 
$\Delta \omega = 2{\rm \pi} \times 0.1\ \mathrm {rad}\ \mathrm {s}^{-1}$. Both this range and the constants 
$\lbrace A_{n_\omega } \rbrace$ are adjusted to reproduce the experimental spectra at large scale, see figure 2. For each run, 
$5\times 10^{7}$ values of 
$\eta (t=0)$ and millions of waves from time series of 
$\eta (t)$ with a time step of 
$0.01\ \mathrm {s}$ are computed from independent drawings of 
$\lbrace a_{n_\omega,n_\theta },\phi _{n_\omega,n_\theta } \rbrace$. The former are used to obtain the data reported in table 1 and figures 2–3 whereas waves are documented in figures 4–8.
The p.d.f.s of the normalized surface elevation 
$\eta /\sigma$ from experiments (symbols) and numerical models (thick coloured lines), compared with (a) a normal law and a Tayfun law and (b) the second- and third-order Gram–Charlier series computed with 
$\mathcal {S} = K-3 = 0.2$.
The p.d.f.s of the normalized wave troughs and crests compared with numerical models (thick coloured lines), a Rayleigh distribution (black solid line) and the first nonlinear corrections for unidirectional and narrowband waves (dotted lines).
The p.d.f.s of the normalized wave height 
$H/H_S$ from experiments (symbols) and numerical models (thick coloured lines), compared with a Rayleigh distribution.
Standard deviation 
$\sigma$, skewness 
$\mathcal {S}$, kurtosis 
$K$, dimensionless spectral bandwidth 
$\nu$, mean frequency 
$f_0$, peak frequency 
$f_{p}$, Tayfun frequency 
$f_T$ and steepness 
$\varepsilon _T$ based on 
$f_T$. Here, ‘exp’ denotes experimental measurements, ‘num’ numerical models and ‘th’ theoretical estimates given by (4.1) and (4.3) and computed based on the experimental PSD. The number of waves 
$N_{tot}$ and RWs 
$N_{RW}$, defined as 
$H>2 H_S$ with 
$H_S = 4 \sigma$, are also reported.
4. Moments
The first moments of 
$\eta (t)$ from experiments and numerical models are reported in table 1. The standard deviation 
$\sigma = \langle \eta ^{2} \rangle ^{1/2}$ is found to increase with the forcing amplitude, while the skewness 
$\mathcal {S} = \langle \eta ^{3} \rangle / \sigma ^{3}$ and the kurtosis 
$K = \langle \eta ^{4} \rangle / \sigma ^{4}$ remain roughly constant. Other characteristics of sea states are computed, namely the dimensionless spectral bandwidth 
$\nu = (m_0m_2/m_1^{2}-1)^{1/2}$, with 
$m_n = \int S_\eta (f) f^{n} \,\mathrm {d}f$ the spectral moments, the mean frequency 
$f_0 = m_1/m_0$, the peak frequency 
$f_{p}$, the Tayfun frequency 
$f_T = f_0/[1+\nu ^{2}(1+\nu ^{2})^{-3/2}]$ discussed later in the manuscript (Tayfun Reference Tayfun1993; Tayfun & Fedele Reference Tayfun and Fedele2007) and the steepness 
$\varepsilon _T = (2 {\rm \pi}f_T)^{2} \sigma /g$ based on 
$f_T$, with 
$g$ the acceleration due to gravity. The dimensionless parameters measured experimentally (
$\mathcal {S}$, 
$K$, 
$\nu$ and 
$\varepsilon _T$) correspond to typical values observed in the ocean, although field measurements yield 
$f_{0,{p},T}= O(0.1)\ \mathrm {Hz}$ and 
$\sigma = O(1)\ \mathrm {m}$ (Christou & Ewans Reference Christou and Ewans2014). This confirms that the wave field under study shares the complex dynamics at work in the ocean while allowing the recording of ten times more waves over the same acquisition time.
The skewness 
$\mathcal {S}$ can be compared with theoretical predictions. The linear model reduces surface elevation to a sum of independent progressive waves of various frequencies and amplitudes (
$\eta ^{(1)}(t$) in (3.1)), for which 
$\mathcal {S}$ vanishes. In the 1960s, Longuet-Higgins computed the second-order nonlinear correction 
$\eta ^{(2)}(t)$ and showed that it only involves non-resonant interactions, mathematically of the form of progressive waves that do not verify the linear dispersion relation, the so-called ‘bound waves’ (Longuet-Higgins Reference Longuet-Higgins1977). The skewness then becomes non-zero and can be inferred from 
$S_\eta (f)$: simplified under the assumption of an isotropic wave field, it reads
\begin{equation} \mathcal{S} = \int \frac{3k_1}{2{\rm \pi} \sigma^{3}} S_\eta(k_1) S_\eta(k_2) I\left( \frac{k_2}{k_1} \right) \mathrm{d}k_{1,2}, \end{equation}
where 
$I$ is an explicit function, see Appendix A. Further, assuming a narrowband frequency spectrum (
$\nu \ll 1$, i.e. 
$f_0=f_T=f_{p}$) numerically yields 
$\mathcal {S} = 2.07 \varepsilon _T$, in contrast to 
$\mathcal {S} = 3 \varepsilon _T$ for unidirectional waves of a narrowband frequency spectrum. The theoretical prediction of 
$\mathcal {S}$ computed from (4.1) together with the experimental PSD 
$S_\eta$ is reported in table 1: it accounts for both numerical models and experimental results.
Several decades later, Janssen built on the canonical transformation introduced in Zakharov (Reference Zakharov1968) to derive the surface elevation up to the next order and to consistently compute the deviation of the kurtosis from three (Janssen Reference Janssen2009). Disentangling resonant and non-resonant interactions, he obtained
\begin{equation} K=3(1+C_4^{dyn}+C_4^{can}), \end{equation}
where 
$C_4^{dyn}$ results from four-wave resonant interactions and only allows analytic expressions for spectra that are narrow in frequency and direction (Fedele Reference Fedele2015; Janssen & Fedele Reference Janssen and Fedele2019). In contrast, 
$C_4^{can}$ is associated with bound waves and can be inferred directly from 
$S_\eta (f)$: for an isotropic wave field,
\begin{equation} C_4^{can}=\int \frac{k_1^{2}}{{\rm \pi}^{2}\sigma^{4}} S_\eta(k_1) S_\eta(k_2) S_\eta(k_3) \psi \left( \frac{k_2}{k_1}, \frac{k_3}{k_1}\right) \mathrm{d}k_{1,2,3},\end{equation}
where 
$\psi$ is another explicit function, see Appendix B. Furthermore, if the spectrum is narrowband in frequency, it reduces to 
$C_4^{can} = 2.75\varepsilon _T^{2}$. The theoretical values of 
$3C_4^{can}$ computed from (4.3) and the experimental PSD 
$S_\eta$ are reported in table 1. They match our numerical models, in which no resonant interaction takes place, but strongly differs from experimental measurements. This demonstrates that four-wave interactions not only generate the low-frequency waves under study but also crucially affect their statistics. Note that a similar conclusion has been reached in a regime of capillary wave turbulence dominated by four-wave interactions (Shats, Punzmann & Xia Reference Shats, Punzmann and Xia2010; Xia, Shats & Punzmann Reference Xia, Shats and Punzmann2010).
5. Probability density functions (p.d.f.s)
The p.d.f.s of experimental and numerical normalized surface elevations 
$f(u=\eta /\sigma )$ are reported in figure 3, along with a normal law of zero mean and unit variance, a Tayfun law and two Gram–Charlier series. The normal distribution describes linear waves and accounts neither for the finite skewness nor for a kurtosis other than three. The Tayfun law corresponds to unidirectional and narrowband waves with second-order nonlinearities (Tayfun Reference Tayfun1980), see Appendix C for its analytic expression. It only depends on the steepness 
$\varepsilon _T$ and has been shown empirically to provide a fair estimate of the tails of 
$f(u)$ for isotropic and broadbanded waves as well, provided that 
$\varepsilon _T$ is artificially tuned to generate the observed skewness (0.24 in the case of figure 3(a)) (Aubourg et al. Reference Aubourg, Campagne, Peureux, Ardhuin, Sommeria, Viboud and Mordant2017; Falcon et al. Reference Falcon, Michel, Prabhudesai, Cazaubiel, Berhanu, Mordant, Aumaître and Bonnefoy2020). It is found here to fit the tails of the numerical p.d.f.s and to underestimate the experimental ones. This difference in the probability of extreme surface elevations translates into the difference in kurtosis discussed before. The p.d.f.s are also compared with the low-order Gram–Charlier series commonly used in theoretical work on surface waves, see Appendix D for definitions. They are reported in figure 3(b) based on the typical experimental values 
$\mathcal {S} = K-3 = 0.2$ from table 1. As observed in Klahn, Madsen & Fuhrman (Reference Klahn, Madsen and Fuhrman2021), they both underestimate large positive values and fail to capture large negative ones (for which the p.d.f. is either undefined, as for the second-order Gram–Charlier approximation, or largely above the experimental data, as for the third-order approximation).
Time series are then analysed in terms of zero down-crossing waves, i.e. events separated by zero crossings (
$\eta = 0$) in which 
$\eta$ assumes negative then positive values (IAHR 1989). By definition, the wave height 
$H$ is the sum of the wave trough 
$\eta _{T}$ (taken positive) and wave crest 
$\eta _{C}$, the duration of the wave being the period 
$T$. In this manuscript, RWs are defined as waves for which 
$H > 2 H_S$, with 
$H_S = 4\sigma$ the significant wave height, whereas large crests are defined by 
$\eta _{C} > 1.25 H_S$. The threshold 
$1.25 H_S=5\sigma$ corresponds to an alternative definition of RWs in the literature (Fedele et al. Reference Fedele, Brennan, de León, Dudley and Dias2016). The numbers of recorded waves and RWs are reported in table 1.
Consider first the p.d.f. of 
$\eta _{C}$ and 
$\eta _{T}$. For unidirectional and narrowband wave fields, they have been explicitly computed by Tayfun up to second-order nonlinearities (Tayfun Reference Tayfun1980), see Appendix E for their analytic expressions. Similar to surface elevation, the p.d.f. of 
$\eta _{C}$ has been empirically found to fit the tails of multidirectional wave fields as well (Soquet-Juglard et al. Reference Soquet-Juglard, Dysthe, Trulsen, Krogstad and Liu2005; Denissenko et al. Reference Denissenko, Lukaschuk and Nazarenko2007; Klahn et al. Reference Klahn, Madsen and Fuhrman2021). Both experimental and numerical p.d.f.s are reported in figure 4 along with the theoretical Rayleigh distribution (
$f_{R}(\xi ) = \xi \exp (-\xi ^{2}/2)$, capturing linear waves) and the Tayfun distributions with the steepness parameter tuned to describe a skewness of 
$0.24$. Our numerical models with bound waves only indicate that the fortuitous agreement between Tayfun's predictions for unidirectional waves and data from isotropic wave fields is restricted to crests. Moreover, one of the main outcomes of this work is that large crests are much more likely to be found experimentally than numerically or expected from the Tayfun law.
The wave height 
$H = \eta _{C}+\eta _{T}$ is then investigated. As routinely observed, the distribution of 
$H/H_S$ as a function of the wave period 
$T$ peaks close to the inverse Tayfun frequency 
$f_T^{-1}$ (Tayfun Reference Tayfun1993; Tayfun & Fedele Reference Tayfun and Fedele2007), see the additional figures in Appendix F. The experimental and numerical p.d.f.s of 
$u=H/H_S$, reported in figure 5, are compared with the Rayleigh distribution 
$f_R(u) = 4 u \exp (-2u^{2})$, which describes narrowband waves with no assumption on directionality and is valid even when the second-order nonlinearities are included (Longuet-Higgins Reference Longuet-Higgins1952; Tayfun Reference Tayfun1980). These data are all found to be similar. The wave height 
$H = \eta _{C}+\eta _{T}$ is therefore not only independent of second-order nonlinearities, as can be shown theoretically, but also seems to be largely independent of higher-order corrections. This is in sharp contrast with the statistics of 
$\eta _{C}$ and 
$\eta _{T}$ detailed above.
6. Shape of large crests
The mean surface elevation at a given position right before/after a large crest occurs (identified as 
$\eta (0)$ with time origin shifted such as the crest manifests at 
$t=0$) is approximated at second order in the joint limit of small amplitude and frequency bandwidth as
\begin{equation} \eta(t) = \eta(0)\left[ \frac{\varPsi(t) + \dfrac{\eta_C \mathcal{F}(t)}{H_s}}{1+ \dfrac{\eta_C \mathcal{F}(0)}{H_s}}\right] ,\end{equation}
where 
$\varPsi (t) = \langle \eta (0) \eta (t) \rangle / \sigma ^{2}$ is the autocorrelation function, 
$\mathcal {F}$ is a function of 
$S_\eta$ detailed in Appendix G and 
$\eta _C$ is the linear component of 
$\eta (0)$ (Fedele & Tayfun Reference Fedele and Tayfun2009). Previous studies have only tested this result in the linear limit in which 
$\eta _C = 0$ (Soquet-Juglard et al. Reference Soquet-Juglard, Dysthe, Trulsen, Krogstad and Liu2005; Klahn et al. Reference Klahn, Madsen and Fuhrman2021). The normalized elevation 
$\eta (t)/\eta (0)$ computed from (6.1) with both 
$\eta _C = 0$ and 
$\eta _C = 1.25 H_s$ is reported in figure 6, along with experimental and numerical values for crests such that 
$\eta _C>1.25 H_S$. Our data confirm that the linear approximation overestimates the depths of the troughs preceding and following the crest, a discrepancy fixed with second-order corrections. However, both theoretical models are symmetric in time reversal (since 
$\varPsi (t) = \varPsi (-t)$ and 
$\mathcal {F}(t) = \mathcal {F}(-t)$) whereas experimental measurements before and after the crest occurs persistently differ. This asymmetry also manifests in steeper slopes before the crests 
$(t<0)$ than after 
$(t>0)$. The numerical simulations of Fujimoto, Waseda & Webb (Reference Fujimoto, Waseda and Webb2019) have shown that, at a fixed time and for directional wave fields, high crests are not symmetric in space as a result of the four-wave resonant interactions not captured by the second-order model reported in (6.1).
Shape of the large crests (
$\eta (0) > 1.25 H_S$) for Run 2. The coloured area corresponds to experiments (mean value 
$\pm$ standard deviation), the black dashed line to numerical models and solid lines to first- and second-order theories. Similar figures for Run 1 and Run 3 are reported in Appendix F.
7. Conclusion
Laboratory experiments with simplified directional spectra provide useful hints about the various processes taking place in the ocean without the usual bias of, e.g. wave breaking regularization in numerical simulations or varying environmental conditions in field measurements. In this study, more than two thousand RWs were observed in statistically homogeneous, isotropic and steady wave fields, allowing the predictions of commonly used theoretical models to be confronted with data in which strongly nonlinear events take place. To highlight the consequences of these high-order nonlinearities, numerical simulations associated with similar PSDs and valid up to second order were carried out. Therefore, they include the leading-order bound wave correction but not the resonant interactions.
The third and fourth normalized moments of surface elevation are compared with theoretical results in which the leading-order bound wave correction is accounted for. These analytic expressions are found to accurately describe the skewness of both experimental and numerical data. However, they significantly underestimate the experimental kurtosis while being in agreement with the numerical ones, evidencing a first consequence of resonant interactions on the statistics. This discrepancy is also manifest in the tails of the normalized surface elevation p.d.f.s.
The surface elevation time series are then split into individual waves whose heights, crests and troughs are analysed. The wave height is found to be robust to high-order effects, the experimental p.d.f.s being similar to the numerical ones and to the Rayleigh distribution. A similar conclusion cannot be drawn regarding the wave crests and troughs, for which large values are much more likely experimentally than numerically, indicating that four-wave resonant interactions strongly affect their statistics. The impact of high-order nonlinearities on large crests is further evidenced through the comparison of their mean shape with first- and second-order theoretical predictions, none of them being able to capture the asymmetry under time reversal. Therefore, the phenomenology of rogue waves crucially depend on how they are defined: high-order nonlinear effects do not seem to play a significant role if the wave height criterion 
$H> 8 \sigma$ is used, whereas for RW depicted as 
$\eta _{C}> 5 \sigma$ (referred to as ‘large crests’ in this paper) they significantly enhance their probability of occurrence. This finding demonstrates the current need for higher-order theoretical models that disentangle troughs and crests.
Note that, as reported in previous studies (Soquet-Juglard et al. Reference Soquet-Juglard, Dysthe, Trulsen, Krogstad and Liu2005; Denissenko et al. Reference Denissenko, Lukaschuk and Nazarenko2007; Aubourg et al. Reference Aubourg, Campagne, Peureux, Ardhuin, Sommeria, Viboud and Mordant2017; Falcon et al. Reference Falcon, Michel, Prabhudesai, Cazaubiel, Berhanu, Mordant, Aumaître and Bonnefoy2020; Klahn et al. Reference Klahn, Madsen and Fuhrman2021), some features of our second-order numerical model of isotropic waves are surprisingly well fitted by theoretical models derived for unidirectional and narrowband wave fields, provided that the single parameter they depend on, the steepness 
$\varepsilon _T$, is tuned to generate the observed skewness. This applies to the tails of the p.d.f.s of both the normalized surface elevation and wave crests, but not to the wave troughs.
Many geophysical processes that are both challenging to model theoretically and to disentangle from other effects in field experiments could benefit from similar investigations with these isotropic nonlinear steady states. This includes, but is not limited to, the impact of waves on mixing and air–sea fluxes, the effect of rain in calming the sea and the effective parameters of random nonlinear waves (diffusion of a pollutant, damping and scattering of a wave train, etc.).
Funding
We thank the technical team at the ECN facilities for their help and support on the experimental set-up. Part of this work was supported by the French National Research Agency (ANR DYSTURB Project No. ANR-17-CE30-0004), and by a grant from the Simons Foundation MPS No. 651463-Wave Turbulence.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Detail on (4.1)
Following Janssen (Reference Janssen2009) and its notations, the third moment of the surface elevation 
$\mu _3$ is related to the standard deviation 
$\sqrt {\mu _2}$ and to the skewness parameter 
$C_3$ through its (51) and (52), that are
\begin{equation} C_3 = \frac{\mu_3}{\mu_2^{3/2}} = \frac{3}{m_0^{3/2}} \int \mathrm{d}\boldsymbol{k}_{1,2} E_1 E_2 \left( \mathcal{A}_{1,2} + \mathcal{B}_{1,2} \right),\end{equation}
where 
$m_0 = \int \mathrm {d}\boldsymbol {k}_1 E_1$ and 
$E(\boldsymbol {k})$ is the first-order spectrum. After lengthy but straightforward computations using various equations of Janssen (Reference Janssen2009), we obtain the transfer coefficients 
$\mathcal {A}_{1,2} (\boldsymbol {k}_1,\boldsymbol {k}_2)$ and 
$\mathcal {B}_{1,2} (\boldsymbol {k}_1,\boldsymbol {k}_2)$ for deep-water gravity waves
$$\begin{gather}
\mathcal{A}_{1,2}= \frac{1}{\sqrt{k_1k_2}} \left[\!
\frac{\left(\sqrt{k_1}+\sqrt{k}_2\right)^{2}(
\boldsymbol{k}_{1} \boldsymbol{{\cdot}} \boldsymbol{k}_2
-k_1k_2)}{(\sqrt{k_1} +\sqrt{k_2})^{2}-\vert
\boldsymbol{k}_1+ \boldsymbol{k}_2 \vert}-
\left(\!\frac{\boldsymbol{k}_{1} \boldsymbol{{\cdot}}
\boldsymbol{k}_2 - k_1k_2 -
\sqrt{k_1k_2}(k_1+k_2)}{2}\!\right)\! \right],
\end{gather}$$
$$\begin{gather}\mathcal{B}_{1,2}= \frac{1}{\sqrt{k_1k_2}}\left[ \frac{(\sqrt{k_1}- \sqrt{k_2})^{2} (\boldsymbol{k}_1\boldsymbol{{\cdot}} \boldsymbol{k}_2 +k_1k_2)}{ (\sqrt{k_1}- \sqrt{k_2})^{2}- \vert \boldsymbol{k}_1 - \boldsymbol{k}_2 \vert} - \left( \frac{ \boldsymbol{k}_{1} \boldsymbol{{\cdot}} \boldsymbol{k}_2 + k_{1}k_2- \sqrt{k_1k_2}(k_1+k_2) }{2} \right) \right]. \end{gather}$$
It can be readily confirmed that this expression of the skewness corresponds to the one initially derived by Longuet-Higgins ((3.11) of Longuet-Higgins (Reference Longuet-Higgins1977) corrected by a misprint of one half). Given that the wave field is assumed isotropic, 
$E(\boldsymbol {k}) \,\mathrm {d}\boldsymbol {k}= S_\eta (k)/(2 {\rm \pi}) \,\mathrm {d}k \,\mathrm {d}\theta$, with 
$S_\eta (k)$ the surface elevation PSD. Moreover, since the transfer coefficients are invariant by a simultaneous rotation of 
$\boldsymbol {k}_1$ and 
$\boldsymbol {k}_2$, (A1) reduces to
\begin{align} C_3 &= \frac{3}{m_0^{3/2}}\int_0^{2{\rm \pi}} \mathrm{d} \theta \iint \mathrm{d}k_1\,\mathrm{d}k_2 \frac{S_\eta(k_1) S_\eta(k_2)}{2 {\rm \pi}\sqrt{k_1k_2}}\left[ \frac{\left(\sqrt{k_1}+\sqrt{k}_2\right)^{2}(\boldsymbol{k}_{1} \boldsymbol{{\cdot}} \boldsymbol{k}_2 -k_1k_2)}{(\sqrt{k_1} +\sqrt{k_2})^{2}-\vert \boldsymbol{k}_1+ \boldsymbol{k}_2 \vert}\right.\nonumber\\ &\quad \left.\vphantom{\frac{\left(\sqrt{k_1}+\sqrt{k}_2\right)^{2}(\boldsymbol{k}_{1} \boldsymbol{{\cdot}} \boldsymbol{k}_2 -k_1k_2)}{(\sqrt{k_1} +\sqrt{k_2})^{2}-\vert \boldsymbol{k}_1+ \boldsymbol{k}_2 \vert}}+ \frac{(\sqrt{k_1}- \sqrt{k_2})^{2} (\boldsymbol{k}_1\boldsymbol{{\cdot}} \boldsymbol{k}_2 +k_1k_2)}{ (\sqrt{k_1}- \sqrt{k_2})^{2}- \vert \boldsymbol{k}_1 - \boldsymbol{k}_2 \vert} - \boldsymbol{k}_{1} \boldsymbol{{\cdot}} \boldsymbol{k}_2 + \sqrt{k_1k_2}(k_1+k_2) \right], \end{align}
with 
$\boldsymbol {k}_1 = k_1 \boldsymbol {e}_x$ and 
$\boldsymbol {k}_2 = k_2 ( \cos \theta \boldsymbol {e}_x + \sin \theta \boldsymbol {e}_y )$. Define a function 
$I$ such that
\begin{align} I(\alpha) &= \int_0^{2{\rm \pi}} \mathrm{d} \theta \left[ \frac{\sqrt{\alpha}\left(1+\sqrt{\alpha}\right)^{2}(\cos \theta - 1)}{(1 +\sqrt{\alpha})^{2}-\sqrt{1+\alpha^{2} + 2 \alpha \cos \theta}}\right.\nonumber\\ &\quad \left.+ \frac{\sqrt{\alpha}(1- \sqrt{\alpha})^{2} (\cos \theta + 1)}{ (1- \sqrt{\alpha})^{2}- \sqrt{1 + \alpha^{2} - 2 \alpha \cos \theta}} + (1+\alpha) \right], \end{align}and (A4) then reads
\begin{equation} C_3 = \frac{3}{m_0^{3/2}}\iint \frac{S_\eta(k_1) S_\eta(k_2) k_1}{2 {\rm \pi}} I\left( \frac{k_2}{k_1} \right)\mathrm{d}k_1\,\mathrm{d}k_2, \end{equation}
which corresponds, with 
$C_3 \rightarrow \mathcal {S}$ and 
$m_0 \rightarrow \sigma ^{2}$ (our notations), to (4.1).
Appendix B. Detail on (4.3)
A similar procedure can be applied to compute the canonical contribution to the kurtosis from (59) of Janssen (Reference Janssen2009),
\begin{equation} C_4^{can} = \frac{4}{m_0^{2}} \int E(\boldsymbol{k}_1)E(\boldsymbol{k}_2)E(\boldsymbol{k}_3)\varPsi(\boldsymbol{k}_1,\boldsymbol{k}_2,\boldsymbol{k}_3) \,\mathrm{d}\boldsymbol{k}_1\,\mathrm{d}\boldsymbol{k}_2\,\mathrm{d}\boldsymbol{k}_3, \end{equation}
where 
$\varPsi$ is an explicit interaction coefficient not detailed here. With 
$E(\boldsymbol {k}_i) = S_\eta (k_i)/(2{\rm \pi} ) \,\mathrm {d}k_i \,\mathrm {d}\theta _i$ and 
$\sigma ^{2} = m_0$,
\begin{equation} C_4^{can} = \frac{4}{(2{\rm \pi})^{3}\sigma^{4}} \int S_\eta(k_1)S_\eta(k_2)S_\eta(k_3)\varPsi(\boldsymbol{k}_1,\boldsymbol{k}_2,\boldsymbol{k}_3) \,\mathrm{d}k_1\,\mathrm{d}k_2\,\mathrm{d}k_3\,\mathrm{d}\theta_1\,\mathrm{d}\theta_2\,\mathrm{d}\theta_3. \end{equation}
Since 
$\varPsi$ is invariant under a simultaneous rotation of 
$\boldsymbol {k}_1$, 
$\boldsymbol {k}_2$ and 
$\boldsymbol {k}_3$, a first integration can be performed
\begin{equation} C_4^{can} = \frac{4}{(2{\rm \pi})^{2}\sigma^{4}} \int S_\eta(k_1)S_\eta(k_2)S_\eta(k_3)\varPsi(k_1 \boldsymbol{e}_x,\boldsymbol{k}_2,\boldsymbol{k}_3) \,\mathrm{d}k_1\,\mathrm{d}k_2\,\mathrm{d}k_3\,\mathrm{d}\theta_2\,\mathrm{d}\theta_3, \end{equation}
with 
$\boldsymbol {k}_{2,3}= k_{2,3} ( \cos \theta _{2,3} \boldsymbol {e}_x + \sin \theta _{2,3} \boldsymbol {e}_y )$. Finally, note that the function 
$\varPsi$ is such that
\begin{equation} \varPsi(k_1 \boldsymbol{e}_x,\boldsymbol{k}_2,\boldsymbol{k}_3) = k_1^{2} \varPsi\left(\boldsymbol{e}_x,\frac{\boldsymbol{k}_2}{k_1},\frac{\boldsymbol{k}_3}{k_1}\right), \end{equation}
and define a function 
$\psi$ by
\begin{equation} \psi\left(\alpha, \beta \right) = \int \varPsi\left( \boldsymbol{e}_x, \alpha \left[ \cos \theta_2 \boldsymbol{e}_x + \sin \theta_2 \boldsymbol{e}_y \right], \beta \left[ \cos \theta_3 \boldsymbol{e}_x + \sin \theta_3 \boldsymbol{e}_y \right] \right) \mathrm{d}\theta_2 \,\mathrm{d}\theta_3. \end{equation}
The coefficient 
$C_4^{can}$ then reads
\begin{equation} C_4^{can} = \frac{4}{(2{\rm \pi})^{2}\sigma^{4}} \int S_\eta(k_1)S_\eta(k_2)S_\eta(k_3)k_1^{2} \psi\left(\frac{k_2}{k_1},\frac{k_3}{k_1}\right) \mathrm{d}k_1\,\mathrm{d}k_2\,\mathrm{d}k_3, \end{equation}which corresponds to (4.3).
Appendix C. Tayfun p.d.f. of surface elevation
The p.d.f. of surface elevation can be explicitly computed in the case of a unidirectional and narrowband wave field in which only the first nonlinear correction is computed. However, several misprints make the expression of this p.d.f. difficult to obtain from the literature. In particular, the original derivation of Tayfun (Reference Tayfun1980) must be corrected as follows: his (24) should read
\begin{equation} F_\xi (u) = \left( 2{\rm \pi}\right)^{{-1/2}} \int_{\alpha(u)}^{\infty} {\rm e}^{-\tau^{2}/2}\left\lbrace \mathrm{erf}\left[A(\tau,u) + \beta \right]+\mathrm{erf}\left[A(\tau,u) - \beta \right] \right\rbrace {\rm d}\tau, \end{equation}and his corrected (27) is
\begin{equation} A(\tau, u) = \beta \sqrt{1+\frac{\sqrt{2}\gamma u}{\beta}+ \frac{\tau^{2}}{{2}\beta^{2}}}. \end{equation}
Note also that only approximate expressions of this p.d.f. are reported in Soquet-Juglard et al. (Reference Soquet-Juglard, Dysthe, Trulsen, Krogstad and Liu2005): indeed, their (7) becomes undefined for large negative values of the surface elevation (if their 
$1+2\sigma z <0$, their 
$C(0)$ required in the integral is no longer real valued).
For completeness, the full set of equations required to compute the p.d.f. 
$f(u)$ of the normalized surface elevation 
$u = \eta /\sigma$ (
$\sigma = \langle \eta ^{2} \rangle ^{1/2}$) is reported below
\begin{equation} f(u) = \frac{\mathrm{d} F}{\mathrm{d}u}, \quad F(u) = \frac{1}{\sqrt{2{\rm \pi}}} \int_{\alpha(u)}^{\infty} \exp\left({- \frac{\tau^{2}}{2}}\right) \left[ \mathrm{erf}\left(A(\tau,u)+ \beta \right) + \mathrm{erf}\left(A(\tau,u)- \beta \right) \right] \mathrm{d}\tau, \end{equation}with
\begin{align}
A(\tau,u) &= \beta \sqrt{1+
\frac{\sqrt{2}\gamma u}{\beta}+ \frac{\tau^{2}}{2\beta^{2}}},\nonumber\\
\beta &= \frac{1}{\sqrt{-1 + \sqrt{1+ 4 \sigma^{2} k^{2}}}}, \nonumber\\
\gamma &= \sqrt{\frac{1+ \sqrt{1+4 \sigma^{2} k^{2}}}{2} },
\end{align}and
\begin{equation} \alpha \left(u \geqslant - \frac{\beta}{\sqrt{2}\gamma} \right) = 0, \quad \alpha \left(u <{-} \frac{\beta}{\sqrt{2}\gamma} \right) = \beta\sqrt{-2\left(1+ \frac{\sqrt{2}\gamma u}{\beta}\right)}. \end{equation}Appendix D. Second- and third-order Gram–Charlier series
A theoretical approach to the p.d.f. of surface elevation consists in using low-order Gram–Charlier series. Following Klahn et al. (Reference Klahn, Madsen and Fuhrman2021), we define in this manuscript the second-order approximation as
\begin{equation} f_{GC2}\left(u = \frac{\eta}{\sigma}\right) = \frac{1}{\sqrt{2{\rm \pi}}}{\rm e}^{- u^{2}/2} \left[ 1 + \frac{\mathcal{S}}{6}H_3(u) \right], \quad H_3(u) = u^{3} - 3u, \end{equation}and the third-order one
\begin{equation} f_{GC3}\left(u = \frac{\eta}{\sigma}\right) = \frac{1}{\sqrt{2{\rm \pi}}}{\rm e}^{- u^{2}/2} \left[ 1 + \frac{\mathcal{S}}{6}H_3(u) + \frac{1}{24}(K-3) H_4(u) + \frac{1}{72}\mathcal{S}^{2} H_6(u) \right], \end{equation}with
\begin{equation} H_4(u) = u^{4} - 6u^{2}+3, \quad H_6(u) = u^{6} - 15u^{4} + 45u^{2} - 15. \end{equation}Appendix E. Tayfun p.d.f. of the crests and troughs
For unidirectional and narrowband wave fields, the p.d.f. of crests accounting for second-order nonlinearities reads (Tayfun Reference Tayfun1980)
\begin{equation} f_C(\xi_{C}) = \frac{2 \varepsilon}{-1 + \sqrt{1 + 4 \varepsilon^{2}}} \left(1 - \frac{1}{\sqrt{1+ 2 \varepsilon \xi_{C}}} \right)\exp \left[- \frac{\left({-}1+\sqrt{1+2 \varepsilon \xi_{C}} \right)^{2}}{-1+\sqrt{1+4 \varepsilon^{2}}} \right], \end{equation}
with 
$\xi _{C} = \eta _{C} / \sigma$ and 
$\sigma = \langle \eta ^{2} \rangle ^{1/2}$ (note that the p.d.f. reported in Tayfun (Reference Tayfun1980) considers instead the wave crest normalized by the standard deviation of the linear component). The steepness parameter 
$\varepsilon = \sigma k$, with 
$k$ the central wavenumber of the narrowband wave fields, is in that case related to the skewness 
$\mathcal {S}= 3 \varepsilon + O(\varepsilon ^{3})$. Similarly, for the troughs,
\begin{equation} f_T(\xi_{T}) = \frac{-2 \varepsilon}{-1 + \sqrt{1 + 4 \varepsilon^{2}}} \left(1 - \frac{1}{\sqrt{1- 2 \varepsilon \xi_{T}}} \right)\exp \left[- \frac{\left({-}1+\sqrt{1+2 \varepsilon \xi_{T}} \right)^{2}}{-1+\sqrt{1+4 \varepsilon^{2}}} \right], \end{equation}
with 
$\xi _{T} = \eta _{T} / \sigma$.
Appendix F. Additional wave features
The raw data of the normalized wave height 
$H/H_S$ plotted vs the wave period 
$T$ are reported in figure 7, while the shape of the large crests for Runs 1, 2 and 3 are shown in figure 8.
Experimental normalized wave height as a function of the wave period. The mean value is plotted in thick black and peaks close to the Tayfun period 
$1/f_T$. Vertical lines indicate 
$f_0^{-1}$, 
$f_T^{-1}$ and 
$f_{p}^{-1}$ (
$f_0^{-1}< f_T^{-1}< f_{p}^{-1}$).
Shape of large crests 
$(\eta (0)> 1.25 H_S)$ for Runs 1, 2 and 3. The coloured area corresponds to experiments (mean value 
$\pm$ standard deviation), the black dashed line to numerical models and solid lines to first- and second-order theories. The central figure corresponds to figure 6.
Appendix G. Expected shape of large waves
From (5.7) of Fedele & Tayfun (Reference Fedele and Tayfun2009), define in the deep-water and isotropic limit the function
\begin{align}
\mathcal{F} (t)&= \frac{2}{(2{\rm \pi})^{2}\sigma^{3}} \int S_\eta(k_1) S_\eta(k_2)\nonumber\\
&\quad \times \left[(\mathcal{A}_{1,2} + \mathcal{B}_{1,2}) \cos(\omega_1t)
\cos(\omega_2t)\right.\nonumber\\
&\quad \left.- (\mathcal{A}_{1,2} - \mathcal{B}_{1,2})
\sin(\omega_1t) \sin(\omega_2t)\right] {\rm d}k_1 \,{\rm
d}k_2 \,{\rm d}\theta_1 \,{\rm d}\theta_2,
\end{align}
with 
$\omega _{1,2} = \sqrt {g k_{1,2}}$. A first angular integration can be performed to obtain
\begin{align} \mathcal{F} (t)&= \frac{1}{{\rm \pi}\sigma^{3}} \int S_\eta(k_1) S_\eta(k_2)\nonumber\\ &\quad \times\left[\!(\mathcal{A}_{1,2} + \mathcal{B}_{1,2}) \cos(\omega_1t) \cos(\omega_2t) - (\mathcal{A}_{1,2} - \mathcal{B}_{1,2}) \sin(\omega_1t) \sin(\omega_2t)\!\right] {\rm d}k_1 \,{\rm d}k_2 \,{\rm d}\theta, \end{align}
with 
$\boldsymbol {k}_1 = k_1 \boldsymbol {e}_x$ and 
$\boldsymbol {k}_2 = k_2 (\cos \theta \boldsymbol {e}_x + \sin \theta \boldsymbol {e}_y )$. Further assume
\begin{align} J(\alpha) &= \sqrt{\alpha} \int_0^{2{\rm \pi}} \mathrm{d} \theta \left[ \frac{\left(1+\sqrt{\alpha}\right)^{2}(\cos \theta - 1)}{(1 +\sqrt{\alpha})^{2}-\sqrt{1+\alpha^{2} + 2 \alpha \cos \theta}}\right.\nonumber\\ &\quad \left. - \frac{(1- \sqrt{\alpha})^{2} (\cos \theta + 1)}{ (1- \sqrt{\alpha})^{2}- \sqrt{1 + \alpha^{2} - 2 \alpha \cos \theta}} + 2 \right], \end{align}to obtain
\begin{align} \mathcal{F}(t) &= \frac{1}{{\rm \pi} \sigma^{3}} \int S_\eta(k_1) S_\eta(k_2) k_1\nonumber\\ &\quad \times \left[ I\left( \frac{k_2}{k_1}\right) \cos(\omega_1t) \cos(\omega_2t) - J\left( \frac{k_2}{k_1}\right) \sin(\omega_1 t) \sin(\omega_2t) \right] \mathrm{d}k_1 \,\mathrm{d}k_2, \end{align}
which allows simple numerical integration. The elevation profile 
$\eta (t)$ close to a crest of linear elevation 
$\xi _c$ then follows from (5.10) of Fedele & Tayfun (Reference Fedele and Tayfun2009) and reads at leading order
\begin{equation} \eta(t) = \xi_c \varPsi(t) + \frac{\xi_c^{2} \mathcal{F}(t)}{4\sigma}, \end{equation}
with 
$\varPsi (t) = \langle \eta (0) \eta (t) \rangle / \sigma ^{2}$ the autocorrelation function.
Open access
