3.1. Simulations involving recurrence

For our initial test case, we draw inspiration from the previous simulation of Madsen et al. (Reference Madsen, Bingham and Liu2002). Following their set-up, we utilise an initial carrier wave steepness $k_0a_0=0.1$ , with dimensionless depth $k_0h=2\pi$ (i.e. deep water). The perturbations are prescribed with $p=0.2$ (corresponding to $\delta = k_0a_0$ ) and relative steepness $\epsilon = 0.05$ . The initial waves are determined and specified using the streamfunction solution of Fenton (Reference Fenton1988). The computational domain has length $10\lambda _0$ , where $\lambda _0=2\pi /k_0$ is the carrier wavelength, utilising periodic lateral boundary conditions. Simulations utilise $320$ horizontal grid points such that the carrier wave field is resolved by $\lambda _0/\Delta x = 32$ points per wavelength, coupled with eight points distributed in the vertical direction. The time step is set such that $T_0/\Delta t= 50$ , where $T_0=2\pi /\omega _0$ is the carrier wave period. The simulation is carried out over a duration of $500T_0$ .

The primary objective of these simulations is to verify the present fully nonlinear model against the theoretical predictions of McLean (Reference McLean1982b ). Specifically, we focus on assessing the initial growth rate of sidebands as a means of validation. To facilitate a spatial domain comparison, we adopt the method outlined by Benjamin & Feir (Reference Benjamin and Feir1967) and convert our simulated results from the time domain to the spatial domain using the following transformation:

(3.4) \begin{equation} \frac {x}{\lambda _0}=\frac {c_g t}{c_0 T_0}={\frac {1}{2}\left (1+\frac {2k_0h}{\sinh (2k_0h)}\right )\frac {t}{T_0}}, \end{equation}

where $t$ is the time, $c_g$ represents the group velocity of the primary wave train and $c_0=\omega _0/k_0$ , the celerity of the carrier wave.

During the initial growth stage, the theoretical sideband amplitudes will follow (McLean Reference McLean1982a )

(3.5) \begin{equation} \frac {a}{a_0}=\epsilon \exp (\textrm{Im}\{\beta \}\sqrt {gk_0}x/c_g), \end{equation}

where $\textrm{Im}\{\boldsymbol{\cdot }\}$ represents the imaginary part and $\beta$ is the unstable eigenvalue, representing a dimensionless complex frequency in a frame of reference moving with the primary wave, stemming from the analysis of McLean (Reference McLean1982a , Reference McLeanb ). For the present case McLean’s (Reference McLean1982a ) analysis, as implemented in Fuhrman, Madsen & Bingham (Reference Fuhrman, Madsen and Bingham2004), yields the dimensionless complex frequency $\beta =-0.0977+0.00364i$ , where $i$ is the imaginary unit. Figure 1 displays the spatial evolution of the carrier and sideband amplitudes from the present case. It is evident that the initial growth rate of the sidebands is in good agreement with that predicted by McLean’s theoretical analysis. This serves as validation, particularly for the initial stage of the simulation. Following the initial stage, the exponential growth slows, with the lower sideband amplitude slightly exceeding that of the upper sideband, as also found e.g. by Lo & Mei (Reference Lo and Mei1985). Following their respective peaks, the simulation results in a nearly symmetric recurrence cycle, where energy is passed from the sidebands back to the carrier wave, eventually resulting in a situation resembling the initial condition, and the phenomenon repeats. The length of each cycle, defined as the distance between two minima of the carrier wave, measures approximately $138\lambda _0$ . This finding agrees well with the results of Landrini et al. (Reference Landrini, Oshri, Waseda and Tulin1998) and Madsen et al. (Reference Madsen, Bingham and Liu2002), who reported $140\lambda _0$ and $141\lambda _0$ , respectively. With this simulation, we conclude that the present fully nonlinear model is capable of simulating the Benjamin–Feir instability with low initial carrier wave steepness.

Figure 1. Computed spatial evolution of the carrier wave and sideband amplitudes from the simulated Benjamin–Feir instability for $k_0h = 2\pi$ with low initial steepness $k_0 a_0=0.1$ . The theoretical sideband growth rate corresponds to the evolution predicted by McLean (Reference McLean1982b ).

Following the deep-water case, we next consider a simulation in finite water depth with $k_0h =2$ , again starting with an initial carrier wave steepness of $k_0a_0 =0.1$ . The perturbations are prescribed with $p=0.1$ and $\epsilon =0.1$ . The total simulation duration is extended to $900T_0$ , while all other computational parameters are identical to those used in the deep-water case. The dimensionless complex frequency associated with this case is $\beta = -0.0410 + 0.00218i$ . Figure 2 presents the computed spatial evolution of the carrier wave and sideband amplitudes for this finite-depth case, considering the transformation described by (3.4). As shown, the sideband amplitudes initially exhibit exponential growth, closely matching the theoretical prediction of McLean (Reference McLean1982a ), which can be taken as model validation. After reaching their peak values, a recurrence phenomenon develops, characterised by the transfer of energy from the sidebands back to the carrier wave, similar to the behaviour observed in the deep-water recurrence case. These results further demonstrate the capability of the present fully nonlinear model to accurately simulate the Benjamin–Feir instability under varying depth conditions.

Figure 2. Computed spatial evolution of the carrier wave and sideband amplitudes from the simulated Benjamin–Feir instability for $k_0h = 2$ with moderate initial steepness $k_0 a_0=0.1$ . The theoretical sideband growth rate corresponds to the evolution predicted by McLean (Reference McLean1982a ).

3.2. Simulations involving frequency downshift

As demonstrated in the preceding case, a local frequency downshift occurs when the sideband amplitudes are near their respective peaks. Nevertheless, it is crucial to acknowledge that this downshift is temporary in nature, as the modulational process is cyclical. On the contrary, experiments conducted by Tulin & Waseda (Reference Tulin and Waseda1999), leading to wave breaking, have revealed a permanent downshift of the peak frequency. Simultaneously, the amplitudes of both the carrier wave and upper sideband decrease, while the lower sideband experiences a permanent increase in amplitude. In the following, we will examine one of their cases involving a brief wave breaking around the first modulational peak. The water depth is set as $h=1.2$ m, consistent with the value employed in previous studies such as Madsen et al. (Reference Madsen, Bingham and Liu2002) and Li & Fuhrman (Reference Li and Fuhrman2022), rather than the $2.1$ m utilised in the experimental set-up. A carrier wave characterised by $\lambda _0 =1.2$ m (corresponding to $\omega _0 = 1.14$ s $^{-1}$ and $k_0h=2\pi$ i.e. deep water) and $k_0a_0=0.133$ is used. Notably, the chosen steepness $k_0a_0$ exceeds the breaking threshold value of $0.1125$ according to Banner & Tian (Reference Banner and Tian1998) and Henderson et al. (Reference Henderson, Peregrine and Dold1999). Consequently, the occurrence of dissipative breaking would be anticipated from a physical perspective. The perturbations are generated using parameters $\delta = 0.785k_0a_0$ , which corresponds to $p = 0.2$ , from which the stability analysis of McLean (Reference McLean1982a ) yields the unstable eigenvalue $\beta =-0.0963+0.00639i$ , the perturbation strength is set to $\epsilon = 0.03$ . The total simulation time is $200T_0$ . The remaining computed parameters, such as the domain size and resolution, remain consistent with those utilised in § 3.1.

Figure 3. Computed spatial evolution of the carrier wave and sideband amplitudes from the simulated Benjamin–Feir instability with larger initial steepness $k_0a_0=0.133$ . The theoretical sideband growth rate corresponding to the evolution predicted by McLean (Reference McLean1982b ), and T $\&$ W 1999 represents the experimental data of Tulin & Waseda (Reference Tulin and Waseda1999).

The present fully nonlinear simulation is quantitatively validated in figure 3, again making use of (3.4). Here the spatial evolution of amplitudes (carrier wave plus both sidebands) is compared with the experimental data from Tulin & Waseda (Reference Tulin and Waseda1999). We likewise verify the initial exponential growth rate of the sidebands against the theoretical prediction of McLean (Reference McLean1982a ), based on the unstable eigenvalue reported above. Figure 3 demonstrates a good match of the initial sideband growth rate with the theoretical prediction. This match is maintained longer for the lower sideband than the upper. Furthermore, the initial and longer term evolutions of both the carrier wave and the sidebands exhibit reasonable agreement with the experimental data of Tulin & Waseda (Reference Tulin and Waseda1999). Around the vicinity of the first modulational peak, approximately at $x\approx 42\lambda _0$ , the primary wave undergoes a decrease, reaching a local minimum before gradually rising again. However, it does not fully regain its original strength attained prior to wave breaking. This phenomenon has also been confirmed through other experimental studies (Melville Reference Melville1982). It is also observed that the lower sideband consistently maintains the highest amplitude (around $0.75a_0$ ), while both the carrier wave and upper sideband experience fluctuations at much lower levels. Consequently, this confirms the occurrence of a permanent frequency downshift for larger initial carrier wave steepness. Physically, this phenomenon is related to dissipative wave breaking around the modulation peak. The present simulated results are likewise in close agreement with those reported previously, either utilising potential flow models (Madsen et al. Reference Madsen, Bingham and Liu2002) or more advanced computational fluid dynamics models (Li & Fuhrman Reference Li and Fuhrman2022) which directly simulated the wave breaking process.

Further validation of the present model’s ability to simulate strongly nonlinear wave evolution at finite depth is achieved through an additional numerical experiment with parameters $p=0.2, k_0a_0=0.2$ and $\epsilon =0.02$ , while maintaining all other computational parameters identical to the deep-water frequency downshifted case. The unstable eigenvalue $\beta = -0.0764 + 0.00902i$ , obtained from the finite-depth stability analysis of McLean (Reference McLean1982a ), is employed to verify the initial sideband growth rate. As shown in figure 4, the simulated sideband amplitudes initially exhibit exponential growth, closely matching the theoretical prediction up to approximately $x\approx 35\lambda _0$ for the lower sideband. Beyond the first modulational peak, a permanent frequency downshift is established, characterised by a sustained reduction in the carrier wave amplitude, associated with dissipative effects linked to wave breaking, and consistent with behaviour observed under deep-water conditions. These results further validate the present model’s accuracy in capturing steep wave evolution under finite-depth conditions.

Figure 4. Computed spatial evolution of the carrier wave and sideband amplitudes from the simulated Benjamin–Feir instability for $k_0h=2$ with initial steepness $k_0a_0=0.2$ . The theoretical sideband growth rate corresponding to the evolution predicted by McLean (Reference McLean1982a ).

3.3. Modelling energy dissipation

For the long-time simulations considered herein, it is essential to account for potential energy dissipation during wave evolution. To quantify the energy loss, the total mechanical energy is computed following the formulation of Klahn et al. (Reference Klahn, Madsen and Fuhrman2021a , Reference Klahn, Madsen and Fuhrmanb ) as

(3.6) \begin{equation} E = \frac {1}{2}\int _{0}^{L_x}\left ({\varPhi}_s \frac {\partial \eta }{\partial t}+g\eta ^2\right ){\rm d}x, \end{equation}

where $L_x$ denotes the horizontal extent of the computational domain. The spatial evolution of the total energy for the simulations, using (3.4) and discussed in §§ 3.1 and 3.2, is presented in figure 5. Notably, there is no sign of energy loss in the simulations involving recurrence. In contrast, simulations exhibiting frequency downshift preserve $E/E_0\approx 1$ during the initial stage, followed by a rapid decline of approximately $15\,\%$ and $31\,\%$ over a spatial interval of roughly $10\lambda _0$ , for $k_0h=2\pi$ and $k_0h=2$ , respectively. In deeper water (figure 5 $a$ ), stronger dispersion delays sideband saturation and thus limits energy loss, whereas in finite depth (figure 5 $b$ ), reduced dispersion accelerates sideband growth and precipitates larger, earlier energy dissipation. Comparison with figures 3 and 4 confirms that irreversible energy dissipation contributes to the establishment of a permanent frequency downshift.

Figure 5. Spatial evolution of the total mechanical energy for wave fields exhibiting recurrence and frequency downshift, shown for $k_0h = 2\pi$ $(a)$ and $k_0h = 2$ $(b)$ .

4.1. Initial conditions and computational parameters

In the forthcoming simulations with the present model, the initial conditions for $(\eta , {\varPhi}_s)$ are constructed by linearly combining sinusoidal wave components with random phases uniformly distributed in the range of $[0,2\pi ]$ . Each component is characterised by its angular frequency $\omega$ and amplitude chosen based on the Texel, Marsen and Arsloe (TMA) spectrum proposed by Holthuijsen (Reference Holthuijsen2010). The TMA spectrum derives its name from the TEXEL, MARSEN and ARSLOE datasets that were used by Bouws et al. (Reference Bouws, Günther, Rosenthal and Vincent1985). It is characterised by the product of a JONSWAP (Joint North Sea Wave Observation Project) spectrum $S_J(\omega )$ and a depth factor $H(\omega ,h)$

(4.1) \begin{equation} S\left (\omega \right )=S_J\left (\omega \right )\times H\left (\omega ,h\right ). \end{equation}

The JONSWAP spectrum is defined as

(4.2) \begin{equation} S_{J}\left (\omega \right )=S_0\left (\frac {\omega }{\omega _p}\right )^{-5}\exp \left (-\frac {5}{4}\left (\frac {\omega }{\omega _p}\right )^{-4}\right )\gamma ^{\exp \left (-\left (\omega /\omega _p-1\right )^2/\left (2\sigma _s^2\right )\right )}, \end{equation}

where $\omega _p$ is the peak angular frequency of the wave field, $\gamma$ is the peak enhancement factor, $\sigma _s = 0.07$ if $\omega \lt \omega _p$ and $0.09$ otherwise and the constant $S_0$ is defined implicitly through the relation

(4.3) \begin{equation} \int _{0}^{\infty }S_{J}(\omega ){\rm d}\omega =\sigma ^2=\langle \eta ^2\rangle , \end{equation}

where angle brackets denote spatial averaging, $\langle \eta \rangle$ = 0 and $\sigma$ is the standard deviation of the surface elevation. Therefore, the value of $S_{0}$ is determined by the characteristic steepness of the wave field, which is expressed as $\varepsilon = 2k_p\sigma$ . To transform the deep-water JONSWAP spectrum into a spectrum applicable to arbitrary water depths, the depth factor introduced by Holthuijsen (Reference Holthuijsen2010) in (4.1) is defined as

(4.4) \begin{equation} H\left (\omega ,h\right )=\frac {\tanh ^{2}(kh)c}{2c_g},\def\lumina\hspace {0.2cm}c_g=\frac {c}{2}\left (1+\frac {2kh}{\sinh {2kh}}\right ),\def\lumina\hspace {0.2cm}{c=\frac {\omega }{k}}. \end{equation}

To compare with the wave statistics of Liu et al. (Reference Liu, Zhang, Yang and Yao2022), we employ the same spectral parameters and numerical resolution across all cases in the present model. The spectral significant wave height is taken as $H_{m0} = 4\sigma = 0.06$ m, with the peak wave period $T_p=1$ s (corresponding to the peak wavelength $\lambda _p = 1.56$ m, $\omega _p = 2\pi /T_p = 2\pi$ s $^{-1}$ and peak frequency $f_p=1/T_p=1$ Hz), and $\gamma = 6$ . In the computational set-up, we define the dimensions of the computational domain as $L_x = 128\lambda _p$ . The vertical extent of the fluid domain is reduced by using the artificial boundary condition with $b=1.5H_{m0}$ . To discretise the domain, a total of 4096 (horizontal) grid points are used corresponding to 32 points per peak wavelength, coupled with 11 points distributed in the vertical direction. Note that the resolution utilised implies that the spectrum is effectively cut off at $\omega _c\approx 4.86\omega _p$ , which is enough for simulating third-order effects, see (3.24) of Klahn et al. (Reference Klahn, Madsen and Fuhrman2021b ). Following the same procedure described in Klahn et al. (Reference Klahn, Madsen and Fuhrman2021b , Reference Klahn, Madsen and Fuhrmanc ), starting with linearised initial conditions, the nonlinear terms are ramped to fully on over a duration of $10T_p$ . Simulations then continue for a total duration of $200T_p$ . The characteristic wave steepness for the linearised initial conditions are set according to $\varepsilon$ . To characterise the waves in terms of bound wave nonlinearity, the nonlinear parameter defined by Toffoli et al. (Reference Toffoli, Onorato, Babanin, Bitner-Gregersen, Osborne and Monbaliu2007)

(4.5) \begin{equation} K^+=\frac {H_{m0}k_p}{4}\left (\frac {\left (4\tanh (k_ph)+\tanh {(2k_ph)}\right )\left (1-\tanh {(k_ph)^2}\right )}{\tanh {(k_ph)}\left (2\tanh {(k_ph)}-\tanh {(2k_ph)}\right )}+2\tanh {(k_ph)}\right ) \end{equation}

is used. This parameter transitions to the characteristic steepness as $k_ph\rightarrow \infty$ and to the Ursell number as $k_ph\rightarrow 0$ . The values of $k_ph$ , $\varepsilon$ and $K^+$ for each of the cases to be considered are listed in table 1. The time step is set such that $\Delta t=T_p/50$ , with artificial damping applied at each time step to ensure stability. For each case, we have carried out the time integration with 500 independent initial conditions, with the phase of each wave component determined randomly. Each individual simulation considered, running on a single processor, requires approximately three hours to complete in terms of wall-clock time. We note that the simulations have thus cumulatively taken approximately one year of computation time in total. Due to the presence of artificial damping, the total energy $E$ , computed according to (3.6), exhibits a gradual decrease of approximately $4\,\%$ – $5\,\%$ over a duration of $200T_p$ , as shown in figure 6. This level of energy reduction is comparable in magnitude to that reported by Xiao (Reference Xiao2013) and Klahn et al. (Reference Klahn, Madsen and Fuhrman2021a , Reference Klahn, Madsen and Fuhrmanb ).

Table 1. The values of the relative water depth, the characteristic wave steepness and the nonlinear parameter given by (4.5) adopted in the present simulations.

Figure 6. The total mechanical energy as a function of time for wave fields for different values of $k_ph$ .

Figure 7. Evolution of the computed spectra for the cases involving six dimensionless water depths. (a) $k_{p}h = 10$ , (b) $k_{p}h = 4$ , (c) $k_{p}h = 3$ , (d) $k_{p}h = 2$ , (e) $k_{p}h = 1.5$ and (f) $k_{p}h = 1$ .

4.2. Wave spectra and significant wave height

Here, we consider the evolution of the wavenumber spectrum at four different non-dimensional times: $t/T_p = 0$ , 20, 50 and 80, as shown in figure 7. Although the total simulation duration is $200T_p$ , we present results only up to $80T_p$ , as the spectral characteristics exhibit negligible variation beyond approximately $50T_p$ , indicating that the spectrum has reached a statistically steady state within the considered wavenumber range. As expected, the initial spectra exhibit minor discrepancies at various water depths due to the change in dispersion relation induced by variation in water depth. As described in Holthuijsen (Reference Holthuijsen2010), the derivation of the TMA spectrum depends critically on the assumption that the high-frequency tail of the JONSWAP spectrum in deep water is proportional to $f^{-5}$ . Consequently, the spectrum tail in deep water should be close to the power law $k^{-3}$ , as clearly depicted in figure 7. Similar to previous simulations in finite depth (Xiao Reference Xiao2013), our observations reveal that in deeper water ( $k_ph\geqslant 2$ , see figure 7 a–d), the tails of ${\textit{S}}(k)$ are close to $k^{-3}$ and remain time invariant. In shallower water ( $k_ph\leqslant 1.5$ , see figure 7 e, f), the spectral tails at $t/T_p = 0$ are already steeper than $k^{-3}$ . This can be attributed to finite-depth effects, which suppress high-wavenumber (short-wave) components due to enhanced bottom influence and reduced dispersive capacity. Furthermore, the tails at $t/T_p = 20, 50, 80$ exhibit a steeper slope relative to the initial state at $t/T_p = 0$ . A downward shift in the spectral peak is observed in most cases relative to the initial spectrum; however, for $k_ph=1.5$ , the peak frequency remains nearly unchanged. This is because, as $k_ph$ approaches the critical value of 1.363, MI becomes weak and can no longer effectively transfer energy to lower wavenumbers. This finding is in accordance with the results from both experiments conducted by Onorato et al. (Reference Onorato2009) and numerical simulations presented in Dysthe et al. (Reference Dysthe, Trulsen, Krogstad and Socquet-Juglard2003). In these simulation studies, a spectrum shift is already observed on the scale of the Benjamin–Feir instability, and this pattern is also noted in the study of Onorato et al. (Reference Onorato, Osborne, Serio, Resio, Pushkarev, Zakharov and Brandini2002). Furthermore, in the shallowest case ( $k_ph=1$ ), the peak frequency exhibits a rapid decrease during the initial stage, accompanied by the emergence of a secondary wave system at lower frequencies, as illustrated in figure 8. This is attributed to enhanced nonlinear interactions, which become increasingly effective as dispersion weakens in shallow water.

Figure 8. Spectral evolution over time for the case with $k_ph = 1$ .

4.3. Skewness, kurtosis and the PDF of the surface elevation

We now investigate some statistical properties of the surface elevation. In figures 9 and 10, we respectively present the skewness ${\mathcal{S}}=\langle \eta ^3\rangle /\sigma ^3$ and kurtosis ${\mathcal{K}}=\langle \eta ^4\rangle /\sigma ^4$ as a function of non-dimensional time in various water depths. As reference the results from Liu et al. (Reference Liu, Zhang, Yang and Yao2022) are additionally shown in these two figures. As shown in figure 9, ${\mathcal{S}}$ , serving as a descriptor for the vertical asymmetry of the wave profile, increases significantly from the Gaussian value of zero and settles into a relatively steady state of around ${\mathcal{S}}$ = 0.2 over the period of $10T_p$ , due to the nonlinear ramping procedure employed in the present simulations. Beyond this time frame, there is a relatively good agreement between the fully nonlinear model and the HOS method adopted in Liu et al. (Reference Liu, Zhang, Yang and Yao2022), except for the results of $k_ph = 1$ , where the fully nonlinear model yields results approximately $25\,\%$ higher than those of the HOS method (differences for $t/T_p\lt 10$ are again due to the ramping procedure employed, and are thus not substantial). Moreover, the stable values are almost constant, even for various $k_ph$ . This can be attributed to the similar values of $\varepsilon$ for these five cases, see table 1. This is in line with the well-known fact that skewness is proportional to steepness (Fedele & Tayfun Reference Fedele and Tayfun2009).

Figure 9. Computed variation of the skewness ${\mathcal{S}}$ over time with various dimensionless water depths. Asterisks represent HOS results from Liu et al. (Reference Liu, Zhang, Yang and Yao2022), full lines denote results from the fully nonlinear model and vertical dashed lines indicate the end of nonlinear ramping at $t/T_p=10$ . (a) $k_{p}h = 10$ , (b) $k_{p}h = 4$ , (c) $k_{p}h = 3$ , (d) $k_{p}h = 2$ , (e) $k_{p}h = 1.5$ and (f) $k_{p}h = 1$ .

Kurtosis is predominantly influenced by the nonlinear dynamics of free waves (see e.g. Mori & Janssen Reference Mori and Janssen2006; Onorato et al. Reference Onorato2009), which plays a crucial role in the generation of extreme events. The kurtosis ${\mathcal{K}}$ is regarded to be a good indicator of the importance of MI (see e.g. Xiao et al. Reference Xiao, Liu, Wu and Yue2013; Liu et al. Reference Liu, Zhang, Yang and Yao2022). The resulting temporal evolution of ${\mathcal{K}}$ from the fully nonlinear simulations for all five cases are presented as lines in figure 10. In the case of larger water depths (figure 10 a–c), it is evident that ${\mathcal{K}}$ approaches a quasi-steady level after first reaching a local maximum. This indicates that the nonlinear interactions among free waves have been significant in these three scenarios, which has likewise been reported in many studies in deep water (see e.g. Onorato et al. Reference Onorato2009; Toffoli et al. Reference Toffoli, Gramstad, Trulsen, Monbaliu, Bitner-Gregersen and Onorato2010; Xiao et al. Reference Xiao, Liu, Wu and Yue2013). At the instant of maximum ${\mathcal{K}}$ , the wave fields may therefore be considered to be at their most nonlinear state during the simulations. On the other hand, for shallower water depths (figure 10 d, e), it is observed that ${\mathcal{K}}$ converges towards a steady value without necessarily first reaching a clear local maximum. The time at which the roughly steady state is reached depends on water depth. This observation implies that the wave fields considered in these two cases are driven away from their initial Gaussian state by bound wave nonlinearities. It is worth highlighting that, as water depth decreases, the four-wave nonlinear interactions among the free waves become weaker. Assessing all cases considered in this work, it becomes apparent that the quasi-steady value is strongly related to the water depth. More specifically, it tends to be higher for larger water depths. Compared with the results from Liu et al. (Reference Liu, Zhang, Yang and Yao2022), the maximum ${\mathcal{K}}$ from the present fully nonlinear simulations are around $6\,\%$ larger, whereas it is approximately $6\,\%$ smaller for $k_ph=1$ . This difference can likely be attributed to the presence of full nonlinearity.

Figure 10. Computed variation of the kurtosis ${\mathcal{K}}$ over time for various dimensionless water depths. Asterisks represent HOS results from Liu et al. (Reference Liu, Zhang, Yang and Yao2022), full lines denote results from the fully nonlinear model and vertical dashed lines indicate the end of nonlinear ramping at $t/T_p=10$ . (a) $k_{p}h = 10$ , (b) $k_{p}h = 4$ , (c) $k_{p}h = 3$ , (d) $k_{p}h = 2$ , (e) $k_{p}h = 1.5$ and (f) $k_{p}h = 1$ .

Figure 11. Comparison of PDFs computed from data generated using the present fully nonlinear model (Klahn et al. Reference Klahn, Madsen and Fuhrman2021d , circles, with error bars) with simulated results from the third-order HOS method (Liu et al. Reference Liu, Zhang, Yang and Yao2022, asterisks), second-order theory (Fuhrman et al. Reference Fuhrman, Klahn and Zhai2023, black lines, referred to as FKZ23; Tayfun & Alkhalidi Reference Tayfun and Alkhalidi2020, red lines), as well as third- through sixth-order theoretical solutions (dashed lines, following the methodology of Klahn et al. Reference Klahn, Zhai and Fuhrman2024, referred to as KZF24 above). (a) $k_{p}h = 10$ , (b) $k_{p}h = 4$ , (c) $k_{p}h = 3$ , (d) $k_{p}h = 2$ , (e) $k_{p}h = 1.5$ and (f) $k_{p}h = 1$ .

We now consider PDFs of the surface elevation at $t/T_p = 50$ for various water depths. The specific time, $t/T_p = 50$ , was selected to correspond directly with the results presented by Liu et al. (Reference Liu, Zhang, Yang and Yao2022). Although PDFs of the surface elevation inherently evolve with time, and thus the choice of the instant analysed could influence quantitative outcomes, adopting $t/T_p = 50$ ensures consistency and facilitates meaningful comparisons. Other instants would likely yield quantitatively distinct distributions, but the chosen reference time provides a robust basis for evaluating and interpreting our results within an established framework. For convenience, we scale the surface elevation $\eta$ by its standard deviation $\sigma$ , i.e. as $\zeta = \eta /\sigma$ . For each case the numerous space series segments have been collectively analysed, with a $\zeta$ bin size corresponding to 0.2, resulting in the six PDFs (circles) shown in figure 11. Error bars are also included, estimated as $\pm p(\zeta )/\sqrt {N_b}$ , where $N_b$ is the number of samples in each bin, following Onorato et al. (Reference Onorato2009). Also shown for comparison are the PDFs from (i) the third-order HOS results of Liu et al. (Reference Liu, Zhang, Yang and Yao2022) (asterisks), (ii) the (first-order) Gaussian distribution (dotted lines)

(4.6) \begin{equation} p(\zeta ) = \frac {1}{\sqrt {2\pi }}\exp \left (-\frac {\zeta ^2}{2} \right ), \end{equation}

and (iii) the exact second-order theoretical PDF derived recently by Fuhrman et al. (Reference Fuhrman, Klahn and Zhai2023) (black lines)

(4.7) \begin{equation} p(\zeta )= \left (\frac {2}{{\mathcal{S}}}\right )^{1/3} \exp \left [\frac {1}{3{\mathcal{S}}^2} + \frac {\zeta }{{\mathcal{S}}}\right ]{Zi}\left (\chi \right ), \end{equation}

where

(4.8) \begin{align} \chi &= \left (\frac {2}{{\mathcal{S}}}\right )^{1/3}\left (\frac {1}{2{\mathcal{S}}}+\zeta \right ), \end{align}

(4.9) \begin{align} {Zi}(\chi ) &\equiv \begin{cases} {Ai}(\chi ) & \text{for} \, \, \, \chi \geqslant -1.17371, \\ \mbox{Ci}(\chi ) & \text{for} \, \, \, \chi \lt -1.17371, \end{cases} \end{align}

where

(4.10) \begin{equation} \mbox{Ci}(\chi )\equiv \big[{Ai}(\chi )^2+{Bi}(\chi )^2\big]^{1/2}. \end{equation}

Note that the switch to the $\mbox{Ci}(\chi )$ function in (4.9) below the indicated threshold for $\chi$ was suggested by Fuhrman et al. (Reference Fuhrman, Klahn and Zhai2023), provided that ${\mathcal{S}} \leqslant 0.2$ , to eliminate non-physical oscillatory behaviour for negative $\zeta$ in the otherwise theoretical solution. Above, ${Ai}(\chi )$ and ${Bi}(\chi )$ are respectively the Airy functions of the first and second kinds. Figure 11 also presents (iv) the second-order Tayfun & Alkhalidi (Reference Tayfun and Alkhalidi2020) distribution (red lines)

(4.11) \begin{equation} p(\zeta )= \begin{cases} \alpha p_G(\zeta _s)/(1+\epsilon _1\zeta _s) & \text{for} \, \, \, \zeta _s\gt 0 ,\\ \alpha p_G(\zeta _s)/\exp {(\epsilon _1\zeta _s)} & \text{for} \, \, \, -2/\epsilon _1\lt \zeta _s\leqslant 0 ,\\ 0 & \text{for} \, \, \, \zeta _s\leqslant -2/\epsilon _1 ,\end{cases} \end{equation}

where

(4.12) \begin{equation} \zeta = \begin{cases} (\zeta _s+\epsilon _1\zeta _s^2/2-\eta _m)/\sigma _s & \text{for} \, \, \, \zeta _s\gt 0 ,\\ (\zeta _se^{\epsilon _1\zeta _s/2}-\eta _m)/\sigma _s & \text{for} \, \, \, -2/\epsilon _1\lt \zeta _s\leqslant 0, \end{cases} \end{equation}

where $\alpha = \sigma _s/P_G(2/\epsilon _1)$ , $p_G$ is the standard normal PDF, see (4.6), $P_G$ is the standard normal cumulative distribution function, $\epsilon _1 = 0.3377{\mathcal{S}}+0.0174{\mathcal{S}}^2 +0.0259{\mathcal{S}}^3$ , $\eta _m = 0.1687{\mathcal{S}}-0.0012{\mathcal{S}}^2+0.0101{\mathcal{S}}^3$ and $\sigma _s = 1+0.0025{\mathcal{S}}+0.0396{\mathcal{S}}^2+0.0104{\mathcal{S}}^3$ . Additionally shown in figure 11 are (v) numerically determined PDFs at third- through sixth-order, based on the theoretically based solutions of Klahn, Zhai & Fuhrman (Reference Klahn, Zhai and Fuhrman2024), who showed that the theoretical PDFs could be computed to essentially any desired order, which was newly posed as solutions to an ordinary differential equation. The sixth-order results based on their methodology are newly presented in this study, with the required coefficients used in the solution detailed in Appendix A. The statistical moments required for the calculation of the second- through sixth-order theoretical PDFs are as indicated in table 2, calculated from all model results at $t/T_p=50$ . Space series of the surface elevation containing the largest computed rogue waves are likewise provided in figure 12, as examples, where it is seen that these rogue waves often appear to be surrounded by an otherwise rather ordinary wave field for the given conditions. From these PDFs, as well as the depicted space series, it is seen that all fully nonlinear simulations result in isolated rogue waves, with crest elevations in the most extreme cases reaching nearly $\zeta = 8$ i.e. eight standard deviations above the mean.

Table 2. Summary statistical moments of cases considered in the present work ( ${\mathcal{S}}_h$ is the hyperskewness, ${\mathcal{K}}_h$ is the hyperkurtosis and $m_7$ is the seventh statistical moment).

From figure 11, it is observed that the third-order HOS results of Liu et al. (Reference Liu, Zhang, Yang and Yao2022) exhibit noticeable deviations from the present fully nonlinear results, which generally predict significantly increased probability density of both extreme wave crests and troughs, as is clear in both positive and negative tail regions. It is well known that, as $k_ph$ approaches 1.363, the four-wave nonlinear interactions (in the weakly nonlinear regime, at least within the narrow-band approximation) become weaker. Therefore, the two primary factors contributing to this observed discrepancy are seemingly limitations associated with higher-order nonlinear dynamics of free waves and stronger effects arising from bound waves. Additionally, we further observe that none of the current theoretical distributions fully capture the characteristics of the simulated results. To illustrate this, we proceed by comparing our findings with a range of analytical distributions, detailing how each aligns with or diverges from the simulated outcomes. Notably, the heavy positive tailed regions computed from the present fully nonlinear results are reasonably predicted by e.g. the fifth and sixth-order solutions based on Klahn et al. (Reference Klahn, Zhai and Fuhrman2024), although even these under-predict the probability density in these regions, perhaps with the exception of results with $k_ph=1.5$ . This stands in contrast to the findings of Toffoli et al. (Reference Toffoli, Onorato, Babanin, Bitner-Gregersen, Osborne and Monbaliu2007), who reported no significant deviations from second-order theory. The discrepancy is likely attributable to the fact that their data are characterised by a directional spreading, whereas our computations are based on long-crested, unidirectional waves, which represent an idealised scenario providing an upper bound for wave statistics. On the other hand, none of the depicted theoretical PDFs do an adequate job of predicting the probability density of the negative tail region (where the resulting theoretical PDFs turn oscillatory at third and higher orders), which is considered as an open problem. Furthermore, at $k_ph=1$ , ${\mathcal{K}}$ is 2.6244 (see table 2), which is less than 3. Therefore, the theoretical PDFs above second order proposed by Klahn et al. (Reference Klahn, Zhai and Fuhrman2024) are not applicable (and hence not shown), as cumulants used in their theory should be positive.

Figure 12. Example snapshots of the computed surface elevation surrounding the largest crests generated by the fully nonlinear wave model of Klahn et al. (Reference Klahn, Madsen and Fuhrman2021c ). Insets depict the region immediately surrounding the largest crest. Variable $x_p$ denotes the $x$ position of the highest crest. (a) $k_{p}h = 10$ , (b) $k_{p}h = 4$ , (c) $k_{p}h = 3$ , (d) $k_{p}h = 2$ , (e) $k_{p}h = 1.5$ and (f) $k_{p}h = 1$ .

4.4. Wave crest distribution and rogue wave occurrence

Here, we investigate the exceedance probability of the wave crest, defined as the highest elevation of each individual wave with respect to the mean water level using zero-up crossing analysis. In the context of the crest amplitude, linear theory predicts a Rayleigh distribution

(4.13) \begin{equation} P\left (\frac {\eta _c}{H_{m0}}\right )= \exp \left (-8\frac {\eta _c^2}{H_{m0}^2}\right ), \end{equation}

whereas second-order interactions should play a role in causing deviations from this result. Assuming deep water and narrow-bandedness of the spectrum, Tayfun (Reference Tayfun1980) formulated a second-order wave crest distribution, which is also reasonably suitable in intermediate water depth according to Fedele et al. (Reference Fedele, Herterich, Tayfun and Dias2019). The exceedance probability is expressed as

(4.14) \begin{align} P \!\left(\frac {\eta _c}{H_{m0}}\right ) =\exp \left [-\frac {8}{H_{m0}^2k_p^2}\left (\sqrt {1+2k_p\eta _c}-1\right )^2\right ] = \exp \left [-\frac {2}{\varepsilon ^2}\left (\sqrt {1+4\varepsilon \frac {\eta _c}{H_{m0}}}-1\right )^2\right ]\!,\nonumber\\ \end{align}

where $\eta _c$ is the crest height. Extending this framework, Tayfun & Fedele (Reference Tayfun and Fedele2007) later proposed a third-order wave crest distribution model, further improving the description of nonlinearity

(4.15) \begin{align} P\left (\frac {\eta _c}{H_{m0}}\right ) = \exp \left [-\frac {2}{\varepsilon ^2}\left (\sqrt {1+4\varepsilon \frac {\eta _c}{H_{m0}}}-1\right )^2\right ]\left [1+\varLambda \left (\frac {\eta _c}{H_{m0}}\right )^2\left (\left (\frac {2\eta _c}{H_{m0}}\right )^2-1\right )\right ], \nonumber\\ \end{align}

where $\varLambda = \lambda _{40}+2\lambda _{22}+\lambda _{04}$ is a relative measure of third-order nonlinearities. The surface cumulants $\lambda _{mn}$ are defined by

(4.16) \begin{equation} \lambda _{mn} = \frac {\langle \eta ^m\hat {\eta }^n\rangle }{\sigma ^{m+n}}+(-1)^{m/2}(m-1)(n-1),\,\, \text{for}\,\,\, m + n = 4, \end{equation}

where $\hat {\eta }$ is the Hilbert transform of $\eta$ . In figure 13, the theoretical (linear) Rayleigh, second-order Tayfun and third-order Tayfun–Fedele distributions, along with the simulated results of Liu et al. (Reference Liu, Zhang, Yang and Yao2022), are utilised as references for comparing with the numerical results obtained from the present fully nonlinear model at $t/T_p = 50$ .

Figure 13. Exceedance probability of wave crests for various dimensionless depths $k_p h$ . Also shown are the Rayleigh distribution (dotted lines), the second-order Tayfun distribution (full lines), results from the third-order HOS model from Liu et al. (Reference Liu, Zhang, Yang and Yao2022) (asterisks) and those from the present fully nonlinear simulations (circles). (a) $k_{p}h = 10$ , (b) $k_{p}h = 4$ , (c) $k_{p}h = 3$ , (d) $k_{p}h = 2$ , (e) $k_{p}h = 1.5$ and (f) $k_{p}h = 1$ .

The wave crest distributions from the present fully nonlinear simulations exhibit a substantial departure from the Tayfun distribution (4.14) at the tails for cases having larger water depths (see figure 13 a–d), although this solution still notably provides much better accuracy than the Rayleigh distribution. These deviations are primarily influenced by the nonlinear dynamics of free waves, which dominates the statistical characteristics of the wave crests. The aforementioned phenomenon of deviation from the Tayfun distribution has already been substantiated in the context of long–crested waves through experimental and numerical models, as demonstrated in Onorato et al. (Reference Onorato, Osborne, Serio, Cavaleri, Brandini and Stansberg2006) and Onorato et al. (Reference Onorato2009). As we move towards shallower water depth (see figure 13 e), the crest amplitude attenuates, leading to reduced (but still apparent) deviations from the second-order theory. Similar trends with respect to water depth are also observed in the simulations conducted by Liu et al. (Reference Liu, Zhang, Yang and Yao2022), confirming that the primary contributing factor is likely the MI enhanced by third-order nonlinearity. For each case considered in the present work, comparing with the third-order HOS simulation results of Liu et al. (Reference Liu, Zhang, Yang and Yao2022), the results of the fully nonlinear simulations demonstrate much larger probability of large crests, especially for $k_ph\geqslant 1.5$ . This discrepancy can be attributed to the impact of full nonlinearity in the interactions among free waves, which will result in and enhance the MI. It should be noted that, at $k_ph = 1.5$ , Liu et al. (Reference Liu, Zhang, Yang and Yao2022) showed that the difference between their results and the second-order Tayfun distribution almost disappears, which indicates that their simulated wave fields are dominated by bound wave effects. Conversely, the deviation of the present results from both the Tayfun and Tayfun–Fedele distributions is still significant. This seemingly indicates that the nonlinear dynamics of free waves and higher-order bound wave effects remain important in the statistical properties of the fully nonlinear wave fields. The third-order Tayfun–Fedele model yields improved predictions in the tail region compared with the second-order Tayfun model. However, it still underestimates the probability of extreme crests when compared with fully nonlinear simulations. This highlights the limitations of weakly nonlinear approximations in capturing the strongly nonlinear wave dynamics, consistent with the findings of Karmpadakis, Swan & Christou (Reference Karmpadakis, Swan and Christou2019). Furthermore, at $k_ph=1$ , $\lambda _{40}$ takes a negative value (−0.3756), which violates the underlying assumption of the Tayfun–Fedele model that $\lambda _{40}$ must be positive. As such, the model is not applicable in this case and is therefore omitted.

Figure 14. Comparison of the probability of rogue wave occurrence from the fully nonlinear model (circles) with third-order HOS results reported in Liu et al. (Reference Liu, Zhang, Yang and Yao2022) (asterisks). (a) $k_{p}h = 10$ , (b) $k_{p}h = 4$ , (c) $k_{p}h = 3$ , (d) $k_{p}h = 2$ , (e) $k_{p}h = 1.5$ and (f) $k_{p}h = 1$ .

The simplest rogue wave definition is one having $\eta _c \gt 1.25H_{m0}$ . Figure 14 depicts the progression of the probability of rogue wave occurrence $P(\eta _c\gt 1.25H_{m0})$ for various dimensionless water depths. As observed, both the present model and the HOS method truncated at third order reveal an initial increase in $P(\eta _c\gt 1.25H_{m0})$ , attributed to the robust MI of the wave fields. Subsequently, a decline ensues, and eventually a relatively stable stage will occur if the simulation is sufficiently long. Here, we focus on comparing the present results with those obtained using the third-order HOS method (Liu et al. Reference Liu, Zhang, Yang and Yao2022). Notably, the present model exhibits a more rapid increase during the ascending phase. ( $20 \leqslant t/T_p \leqslant 50$ ) and produces higher $P(\eta _c\gt 1.25H_{m0})$ values compared with the third-order HOS results, albeit of the same order of magnitude. For instance, in the deepest case with $k_ph = 10$ , the maximum $P(\eta _c\gt 1.25H_{m0})$ reaches up to 2.93 $\times 10^{-3}$ , significantly exceeding the third-order HOS (Liu et al. Reference Liu, Zhang, Yang and Yao2022) peak of $1.68\times 10^{-3}$ by nearly a factor of two. This highlights the discernible impact of full nonlinearity on MI, further contributing to the likelihood of rogue wave formation. As the depth $k_ph$ decreases from 10 to 1.5, both the initial growth rate and $P( \eta _c\gt 1.25H_{m0})$ decrease, yet the maximum values still significantly exceed the third-order HOS predictions. This further illustrates that the MI induced by the higher-order nonlinearity weakens as water depth decreases.

The time evolution of $P(\eta _c\gt 1.25H_{m0})$ and ${\mathcal{K}}$ is plotted simultaneously in figure 14 to enable comparison and analysis of their relationship. This comparison reveals that the evolution of $P(\eta _c\gt 1.25H_{m0})$ is consistent with that of ${\mathcal{K}}$ . This connection emerges as the kurtosis is closely tied to MI, serving as a primary indicator for the occurrence of rogue waves. Kurtosis as such an indicator has been widely used in deep water (e.g. Mori & Janssen Reference Mori and Janssen2006; Xiao et al. Reference Xiao, Liu, Wu and Yue2013). To better elucidate and quantify the relationship in finite depth, we plot $P(\eta _c\gt 1.25H_{m0})$ for various water depths at distinct times (ranging from $t=10T_p$ to $100T_p$ in $5T_p$ increments) as a function of the excess kurtosis ${\mathcal{K}}-3$ in figure 15. Based on these scatter distributions, we compute a linear fit

(4.17) \begin{equation} P( \eta _c\gt 1.25H_{m0})=0.0026({\mathcal{K}}-3), \end{equation}

with a coefficient of determination ( $R^2$ ) of approximately 0.93 from the fully nonlinear model, surpassing the value obtained from the third-order HOS results of Liu et al. (Reference Liu, Zhang, Yang and Yao2022), who found the fit

(4.18) \begin{equation} P( \eta _c\gt 1.25H_{m0})=0.0010({\mathcal{K}}-3), \end{equation}

with $R^2 = 0.56$ based on their results. This indicates that the present model gives a stronger correlation than the third-order HOS method. Comparison of (4.17) and (4.18) indicates that full nonlinearity increases the probability of rogue waves for a given excess kurtosis by a factor $0.0026/0.0010=2.60$ for the present conditions.

Figure 15. Correlation between the probability of rogue wave occurrence $P(\eta _c\gt 1.25H_{m0})$ and the excess kurtosis ${\mathcal{K}}-3$ . Depicted are results from third-order HOS modelling of Liu et al. (Reference Liu, Zhang, Yang and Yao2022) (asterisks) and the present fully nonlinear simulations (circles), with best-fit lines shown for both.