2.1. The ESBI wave model: configuration, verification and validation

The Cartesian coordinate system is adopted here with $\boldsymbol {x} = (x,y)$ being horizontal and $z$ being vertical coordinates. The still deep-water level is at $z=0$ . Unless otherwise specified, the variables are non-dimensionalised, i.e. the distances $\boldsymbol {x}$ , $z$ , and time $t$ are multiplied by the peak wavenumber $k_p$ and angular frequency $\omega _p = \sqrt {gk_p}$ , respectively. The parameter $g$ is the gravitational acceleration.

The potential flow theory assumes that the fluid is inviscid and irrotational, leading to velocities written as gradients of velocity potential $\phi$ , rescaled by $\sqrt {k_p^3/g}$ . The primary advantage of using the velocity potential is that it is a scalar quantity. Therefore, the number of unknowns is reduced compared with the Euler or Navier–Stokes equations, as the velocity vectors can be obtained directly by calculating the gradient of the velocity potential.

The free-surface boundary conditions for the potential flow wave theory consist of those on the water free surface $z=\eta (\boldsymbol {x},t)$ :

(2.1) \begin{equation} \frac {\partial \eta }{\partial {t} } + \nabla \phi \cdot \nabla \eta - \frac {\partial \phi }{\partial z} = 0, \end{equation}

(2.2) \begin{equation} \frac {\partial \phi }{\partial {t} } + \frac {1}{2}\left [ \left |\nabla \phi \right |^2 + \left (\frac {\partial \phi }{\partial z} \right ) ^2 \right ] + \eta = 0, \end{equation}

where $\nabla =(\partial _x,\partial _y)$ is the horizontal gradient operator.

Equations (2.1) and (2.2) are identical to the canonical pair derivable from the Hamiltonian water-wave system (Zakharov Reference Zakharov1968). They can be rewritten as a skew-symmetric form (Fructus et al. Reference Fructus, Clamond, Grue and Kristiansen2005a ):

(2.3) \begin{equation} \frac {\partial \varPsi }{\partial t} + \mathcal {A}\varPsi = \mathcal {N}, \end{equation}

where

(2.4) \begin{eqnarray} \varPsi = \begin{pmatrix} k \mathcal {F}\left \{ \eta \right \} \\ k\omega \mathcal {F}\left \{ \widetilde {\phi } \right \} \\ \end{pmatrix},\quad \mathcal {A} = \begin{bmatrix} 0 & - \omega \\ \omega & 0 \\ \end{bmatrix}, \quad \mathcal {N} = \begin{pmatrix} k\left ( \mathcal {F}\left \{ V \right \} - k\mathcal {F}\left \{ \widetilde {\phi } \right \} \right ) \\ \frac {k\omega }{2} \mathcal {F}\left \{ \frac {\left ( V + \nabla \eta \cdot \nabla \widetilde {\phi } \right )^{2}}{1 + \left | \nabla \eta \right |^{2}} - \left | \nabla \widetilde {\phi } \right |^{2} \right \} \\ \end{pmatrix}, \end{eqnarray}

and $\widetilde {\phi }=\phi (\boldsymbol {x},z=\eta ,t)$ denotes the velocity potential at the free surface, $V = \sqrt {1 + |\nabla \eta |^{2}}\partial _n\phi$ is the vertical velocity of the surface elevation, while $\mathcal {F}\{\psi \}$ is the Fourier transform of quantity $\psi$ in physical space defined as

(2.5) \begin{equation} \mathcal {F}\{\psi \} = \iint _{-\infty }^{\infty } \psi {\rm e}^{-\mathrm {i} \boldsymbol {k} \cdot \boldsymbol {x}}\, \mathrm {d} \boldsymbol {x}, \end{equation}

with $\mathcal {F}^{-1}\{\psi \}$ being its inverse transform, $\mathrm {i}=\sqrt {-1}$ and $k = |\boldsymbol {k}|$ and $\omega = \sqrt {k}$ are dimensionless wavenumber and frequency in Fourier space. The Fourier transform is implemented numerically by using the fast Fourier transform.

Equation (2.3) can be further reformulated as

(2.6) \begin{eqnarray} \varPsi \left ( t \right ) = {\rm e}^{- \mathcal {A}\left ( t - t_{0} \right )} \left [ \varPsi \left ( t_{0} \right ) + \int _{t_{0}}^{t}{{\rm e}^{\mathcal {A}\left ( t - t_{0} \right )}\mathcal {N}\, \mathrm {d} t} \right ], \end{eqnarray}

and can be used as the prognostic equation for updating unknowns $\eta$ and $\widetilde {\phi }$ in time with the integration term evaluated by using a six-stage fifth-order Runge-Kutta method with embedded fourth-order solution (Clamond et al. Reference Clamond, Fructus and Grue2007; Wang & Ma Reference Wang and Ma2015). In general, to keep the difference below $1\,$ %, the time step size is automatically adjusted to about 1/20 of the peak wave period. For large spatio-temporal simulations of strongly nonlinear waves, a tolerance of $0.1\,$ % is selected corresponding to a time step size of about 1/50 of the peak wave period, which applies to the crossing sea simulations in the present study.

In order to update the solutions $(\eta ,\widetilde {\phi })$ , the velocity $V$ needs to be diagnosed by solving the boundary integral equation of Green’s theorem. A successive approximation approach can be adopted, and the total vertical velocity is expressed as $V = \sum_{m} V_m$ , where $m$ represents the order of the nonlinearity $\mathcal {O}(\varepsilon ^m)$ . Here, the expansion parameter $\varepsilon$ denotes the wave steepness. For simplicity, the recurrence formula for estimating $V_m$ in the fully nonlinear ESBI wave model in deep-water starts with

(2.7) \begin{equation} \mathcal {F}\left \{ V_{m=1} \right \} = k \mathcal {F}\left \{ \widetilde {\phi } \right \}. \end{equation}

Then, the remaining velocities in the Fourier domain are calculated by (Wang Reference Wang2025)

(2.8) \begin{eqnarray} \mathcal {F}\left \{ V_{m} \right \} &=& \sum _{j=1}^{m-1} - \frac {k^j}{j!}\mathcal {F}\left \{ \eta ^jV_{m-j} \right \} - \frac {k^{m-2}}{(m-1)!} \mathrm {i}\boldsymbol {k} \cdot \mathcal {F}\left \{ \eta ^{m-1}\nabla \widetilde {\phi } \right \} , \end{eqnarray}

for $m \geqslant 2$ . In this pseudo-spectrum method, the $2/(m+1)$ -rule is used here for anti-aliasing treatment, which is equivalent to the zero-padding method (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1987). We emphasise that a smoothing technique is not required here, and the present model is very stable for cases without appearance of breaking waves. The model has been comprehensively verified and validated for a variety of simulated strongly nonlinear and steep wave phenomena, crossing seas and laboratory experiments (Wang et al. Reference Wang, Ma and Yan2018, Reference Wang, Ma, Yan and Liang2021; Wang Reference Wang2025).

Before starting our study, we recall the commonly used definition of a RW, namely having a crest height $\eta _c$ exceeding at least 1.25 times the significant wave height $H_s$ , i.e. $\eta _c\gt 1.25 H_s$ . In a Gaussian sea state, $H_s$ is approximated as being four times the standard deviation of the entire water surface elevation. This threshold condition is commonly used and describes a significant deviation from the average wave height (Kharif et al. Reference Kharif, Pelinovsky and Slunyaev2008; Gramstad et al. Reference Gramstad, Bitner-Gregersen, Trulsen, Borge and Carlos2018; Mori et al. Reference Mori, Waseda and Chabchoub2023). In this context, we define $t_{LS}$ as the duration or lifetime of a series of observed sequences of extreme waves belonging to the same RW event. This is clarified and discussed later in the article.

The numerical set-up is described next. The computational domain of the simulations covers 40 $\times$ 40 peak wavelengths, and is resolved into 1024 $\times$ 512 collocation points in $x$ (along-wave) and $y$ (cross-wave) directions, respectively. The selected domain size and resolution in space ensure that the Fourier modes up to seven and three times peak wavenumber in the $x$ and $y$ directions, respectively, are aliasing-free. The reference sea state and wave surface condition is based on the JONSWAP spectrum (Hasselmann et al. Reference Hasselmann1973) with different peakedness factors $\gamma = 1,\ 2,\ 3,\ 4,\ 6$ and $9$ and crossing wave field of steepness $k_pH_s = 0.28,\ 0.15$ and $0.06$ . On the selection of wave parameters, given that limited works have studied the effect of spectrum bandwidth of crossing sea RWs (Wang et al. Reference Wang, Ma, Yan and Liang2021), the selected peakedness values ensure that the spectrum widths in the current work are covered from broadband to narrowband conditions (Goda Reference Goda2010). The selected steepness also covers from large to small steepness cases, so that a wide range of Benjamin–Feir indices and kurtosis values can be expected (Mori & Janssen Reference Mori and Janssen2006; Janssen & Bidlot Reference Janssen and Bidlot2009) which also affects the RW formation. Nevertheless, to maximise the effect of nonlinearity and thus the occurrence of RWs, we mainly focus our analysis on the highest steepness $k_pH_s = 0.28$ case, which is consistent with, but not limited to, the wave condition previously adopted by Wang et al. (Reference Wang, Ma, Yan and Liang2021). To avoid the effects of potential wave breaking and high steepness values which can lower the accuracy of our model (Wang Reference Wang2025), we also tested lower steepness values for validation purposes. The two crossing wave systems with the same peak frequency $f_p=1$ Hz, peakedness parameter $\gamma$ , significant wave height $H_s$ and random phases cross-interact at an aperture angle of $28^\circ$ , $40^\circ$ and $53^\circ$ , in which the $40^\circ$ cases are mainly investigated, since these correspond to the most hazardous angle leading to the highest probability distribution tail (Toffoli et al. Reference Toffoli, Bitner-Gregersen, Osborne, Serio, Monbaliu and Onorato2011; Cavaleri et al. Reference Cavaleri, Bertotti, Torrisi, Bitner-Gregersen, Serio and Onorato2012; Bitner-Gregersen & Toffoli Reference Bitner-Gregersen and Toffoli2014). In table 1 we show the parameters of all 11 cases simulated by the fully nonlinear ESBI and weakly nonlinear CNLS frameworks (as introduced in § 2.2) in the data analysis results.

Table 1.

Table of all 11 cases tested numerically in § 3 with varying crossing sea angle $\beta$ , wave steepness $k_pH_s$ , JONSWAP peakedness factor $\gamma$ , simulation method(s) adopted and corresponding figures for analysis and comparison purposes. Note that in all cases, the simulation frequency is fixed at $f_p=1$ Hz.

To generate desired waves at $x=0$ as the boundary condition, a pneumatic wavemaker (Clamond et al. Reference Clamond, Fructus, Grue and Kristiansen2005) is adopted in the current study, which is a localised oscillating pressure on the boundary regime of the water surface to mimic the laboratory wavemaker. We refer the reader to the work by Wang et al. (Reference Wang, Ma, Yan and Liang2021) for further details on numerical wave generation. Mathematically, the pressure is the summation of all wave components after a bell-shaped transformation. According to our previous work (Wang Reference Wang2025), when given $J$ wave components in total, while the $j{\textrm {th}}$ component has amplitude $a_j$ , angular frequency $\omega _j$ , wavenumber $k_j$ , direction $\theta _j$ and phase $\varphi _j$ , the transformation from wave components to the pressure is written as

(2.9) \begin{equation} p_g=\sum _{j=1}^Ja_j \cos (\theta _j) \sqrt {\frac {{\rm e}}{2\pi }}{\rm e}^{-(k_jx)^2/2}\sin [-k_j\sin (\theta _j)y+\omega _jt+\varphi _j]. \end{equation}

The wave generation zone is deployed along $x=0$ and absorbed at a distance near the other end. Each simulation lasts for 1000 peak periods $T_p$ (equivalent to a typical three-hour sea state in a realistic ocean with around 10-second peak period), and four realisations are performed; thus, each case carried out by the ESBI simulations has a total length of 4000 $T_p$ .

A snapshot of the simulated cross free surface as described above is shown in figure 1. Several large-amplitude wave groups can be observed in directional space at a particular instant of time.

Figure 1.

An exemplified snapshot of a simulated crossing sea surface elevation for case 3 in table 1.

From equation (2.9), it is clear that linear waves will be generated from the domain boundary, whereas the second and higher harmonics appear almost instantaneously when the waves propagate away from the generation zone; see for instance experimental studies (He et al. Reference He, Ducrozet, Hoffmann, Dudley and Chabchoub2022a ,Reference He, Witt, Trillo, Chabchoub and Hoffmannc). Hereafter, the velocity potential and free-surface elevation satisfy the fully nonlinear kinematic and dynamic boundary conditions, attributed to the fully nonlinear formalism and numerical approach for solving the boundary integral equation.

Figure 2.

Spatio-temporal distribution of crossing sea RW locations across the ESBI numerical domain from case 3 in table 1. The sampling duration is taken as 250 wave periods for clarity, and both sides (each 10 wavelengths wide) of the numerical domain are neglected in the analysis due to the absence of crossing wave field. Each circle represents one RW, and its colour indicates the appearance time of the current RW.

The generated wave field under the fully nonlinear ESBI framework is shown to rapidly evolve to a fully developed state after approximately 20 wavelengths of wave propagation for crossing wave trains according to the previous kurtosis analysis (Wang et al. Reference Wang, Ma, Yan and Liang2021). In figure 2, the long-living and localised RW structures can be clearly observed at a distance far from the wave generation zone, located at $x=0$ , and appear throughout the rest and majority of the domain.

We further emphasise that wave breaking is inevitable under wave condition $k_pH_s = 0.28$ . To stabilise the simulations, a low-pass filter is employed to suppress the breaking (Xiao et al. Reference Xiao, Liu, Wu and Yue2013). This filter is shown to well represent the energy dissipation quantitatively over a broad range of wave steepness, breaker types and directional spreading. We refer to a previous work (Wang Reference Wang2025) for a direct comparison of the model simulations with laboratory experiments, involving the evolution of kurtosis and wave crest exceedance probability trends. In order to fully exclude the possible effect of breaking RWs for very steep states, cases 7–11 in table 1 with reduced wave steepness have been adopted.

2.2. The hydrodynamic CNLS

The purpose of introducing the CNLS is twofold. It is adopted to compare the strong RW localisations, as obtained from the fully nonlinear ESBI simulations, with a weakly nonlinear framework which applies for cross sea modelling (Cavaleri et al. Reference Cavaleri, Bertotti, Torrisi, Bitner-Gregersen, Serio and Onorato2012), and to test if a weakly nonlinear wave framework is sufficient to characterise all measured RWs observed in the fully nonlinear simulations. We carried out benchmarking numerical simulations by means of the CNLS, as previously reported and parametrised for water waves (Okamura Reference Okamura1984; Onorato et al. Reference Onorato, Osborne and Serio2006a ) and mentioned in § 1. The two-wavetrain deep-water coupled framework is written as

(2.10) \begin{align} &\frac {\partial u_1}{\partial t}+c_x\frac {\partial u_1}{\partial x}+c_y\frac {\partial u_1}{\partial y}-\mathrm {i}\alpha _x \frac {\partial ^2 u_1}{\partial x^2}-\mathrm {i}\alpha _y \frac {\partial ^2 u_1}{\partial y^2}+\mathrm {i}\alpha _{xy} \frac {\partial ^2 u_1}{\partial x \partial y}+\mathrm {i}(\xi |u_1|^2 + 2\zeta |u_2|^2 ) u_1=0,\nonumber\\ &\frac {\partial u_2}{\partial t}+c_x\frac {\partial u_2}{\partial x}-c_y\frac {\partial u_2}{\partial y}-\mathrm {i}\alpha _x \frac {\partial ^2 u_2}{\partial x^2}-\mathrm {i}\alpha _y \frac {\partial ^2 u_2}{\partial y^2}-\mathrm {i}\alpha _{xy} \frac {\partial ^2 u_2}{\partial x \partial y}+\mathrm {i}(\xi |u_2|^2 + 2\zeta |u_1|^2) u_2=0, \end{align}

where the coefficients are defined as

(2.11) \begin{align} &c_x=\frac {\hat {\omega }}{2\hat {k}^2}\kappa ,\ c_y=\frac {\hat {\omega }}{2\hat {k}^2}\iota ,\nonumber \\ \alpha _x=\frac {\hat {\omega }}{8\hat {k}^4}(2\iota ^2-\kappa ^2&),\ \alpha _y=\frac {\hat {\omega }}{8\hat {k}^4}(2\kappa ^2-\iota ^2),\alpha _{xy}=-\frac {3\hat {\omega }}{4\hat {k}^4}\iota \kappa , \nonumber\\ \xi =\frac {1}{2}\hat {\omega }\hat {k}^2,\ \zeta =\frac {\hat {\omega }}{2\hat {k}}&\frac {\kappa ^5-\kappa ^3\iota ^2-3\kappa \iota ^4-2\kappa ^4\hat {k}+2\kappa ^2\iota ^2\hat {k}+2\iota ^4\hat {k}}{(\kappa -2\hat {k})\hat {k}},\nonumber\\ &\kappa =\hat {k} \cos {\frac {\beta }{2}},\ \iota =\hat {k} \sin {\frac {\beta }{2}}. \end{align}

Here, $u_1(x,y,t)$ and $u_2(x,y,t)$ are two crossing complex wave envelopes with wavenumbers of $k_1=(\kappa ,\iota )$ and $k_2=(\kappa ,-\iota )$ , respectively, as mentioned in § 1; $\beta$ is the crossing angle; and the dimensional angular frequency $\hat {\omega }$ and the modulus of the wavenumbers $\hat {k}$ obey the deep-water dispersion relation $\hat {k}=\hat {\omega } ^2/g$ . Note that when $\beta =0$ and $u_2$ are inactive, equation (2.10) is uncoupled and naturally reduces to the classical NLS (Zakharov Reference Zakharov1968). The first-order approximation of a two-wave-field crossing elevation $\eta (x,y,t)$ is given by

(2.12) \begin{equation} \eta (x,y,t)=\frac {1}{2}\left (u_1(x,y,t){\rm e}^{\mathrm {i}(\kappa x+\iota y-\hat {\omega } t)} +u_2(x,y,t){\rm e}^{\mathrm {i}(\kappa x-\iota y-\hat {\omega } t)}+\text{c.c.}\right ), \end{equation}

where $\text{c.c.}$ denotes the complex conjugate. To ensure highest numerical accuracy, fourth-order Runge-Kutta and pseudo-spectral methods (Yang Reference Yang2010) are adopted to advance equation (2.10) in time. The two equations are coupled in a staggered manner, as already adopted by He et al. (Reference He, Slunyaev, Mori and Chabchoub2022b ) to validate laboratory observations. The corresponding numerical results are reported and compared with the fully nonlinear ESBI results in § 3.3. Note that the same numerical domain size and time have been adopted in the CNLS simulations as in the ESBI simulations.

3.1. Spatio-temporal evolution of long-lived directional coherent RW structures

To further examine whether the observed long-lived localised wave structures experience both significant growth and decay in wave height, and thus obey the definition of a RW, i.e. being localised in both space and time, and intuitively speaking ‘appearing from nowhere and disappearing without a trace’ (Akhmediev et al. Reference Akhmediev, Ankiewicz and Taki2009), we depict in figure 4 fully nonlinear numerical results of the averaged long-lifespan RW elevation field and its corresponding directional spectra $20T_p$ before the peak, during the maximal focusing and $20T_p$ after the peak.

Figure 4.

Spatio-temporal evolution of the averaged directional RW events generated from case 5 in table 1 with long lifespan ( $t_{LS}\gt 35 T_p$ ) and at three different stages of extreme focusing evolution: $20T_p$ before reaching the peak (left-hand panels), at the focusing peak ( $t=0$ ) (centre panels) and $20T_p$ after the peak (right-hand panels). Upper row: the averaged elevation field. Lower row: corresponding average wave energy spectra.

In figure 4, we can observe a strong wave-amplitude focusing followed by a decay in both $x$ and $y$ directions, as well as a severe broadening and narrowing of the directional spectrum. The current observations also show that a clearly detectable long-lived breathing-type process can occur in directional seas, featured by a directional spatial localisation at $t=0$ as well as the temporal localisation throughout the evolution process $20T_p$ before and after the peak along the positive $x$ direction. Thus, these observed directional and particular RW structures are likely to be a result of nonlinear interactions, which indeed yield a full localisation in time and directional space. To the best of our knowledge, such a pulsating wave phenomenon with longevity features has so far not been reported and not discussed in detail by means of fully nonlinear numerical simulations.

To further characterise and distinguish these RW events in such irregular cross sea states, it is important to extract key statistical characteristics, which are analysed and discussed as next.

3.2. Statistical analysis of RW events in crossing seas

In the following, we extract all independent RW events from the fully nonlinear ESBI numerical data and calculate the corresponding probability density function (PDF) with a categorisation with respect to both $t_{LS}$ and maximum wave crest height $\eta _{{max}}$ of the extreme events. Figure 5 shows the PDFs of all independent events, which are generated with respect to the maximum amplification factor and the normalised lifespan for all six JONSWAP- $\gamma$ cases with the identical $H_s$ values, as mentioned earlier.

Figure 5.

Statistical PDF results of the correlation between the amplification factors, defined as $\eta _{{max}}/H_s$ , their lifespans $t_{LS}$ as RWs along the positive $x$ -directions and the PDF while considering different JONSWAP spectral peakedness parameter $\gamma$ values, i.e. cases 1–6 in table 1.

The results indicate that most independent events have a short lifespan, suggesting that the wave superposition principle is a dominant focusing mechanism for our modelled cross sea states. That being said, around one in $10^5$ to $10^3$ events exhibits significantly longer lifespans, depending on the JONSWAP peakedness parameter $\gamma$ value considered, with a trend of increasing amplification factor by gradually narrowing the initial energy spectrum. Such observations indicate the existence of a nonlinear focusing mechanism responsible for such wave amplifications, despite having a low probability of emergence. However, the fact that such RWs are recurrent, i.e. experience recurrent focusing, these can indeed still pose a significant threat in the ocean.

To interpret these nonlinear RW events, we analyse in figure 6 the events according to their lifespans and the spectral directionality measured by the full area at half maximum (FAHM), which is simply a two-dimensional generalisation of the well-known concept of full width at half maximum measuring the bandwidth of a power spectrum. Similar concepts or approaches have been involved in previous works when studying the nonlinear evolution of directional wave fields (Toffoli et al. Reference Toffoli, Onorato, Bitner-Gregersen and Monbaliu2010; Osborne & Ponce de León Reference Osborne and Ponce de León2017), yet may have not considered a variety of cases or quantified the physical features, which are different from the current dual bimodal system studied. Notably, we did not introduce any initial directional spreading and the two wave trains are perfectly unidirectional initially. Since the distributions of the RWs are highly non-uniform along the $t_{LS}/T_p$ axis, a good linear fit and reasonable clustering can be achieved upon averaging all original RW elevation data along the $t_{LS}$ axis within every certain interval. It is chosen to be $t_{LS}/T_p=0.6$ in our case.

Figure 6.

The correlation between the RWs’ lifespans $t_{LS}/T_p$ and their directional FAHM factor. The left $y$ axis represents the FAHM values of all RWs visualised by dots, and the right $y$ axis shows the FAHM values of the averaged RW elevation fields (instead of the averaged FAHM value) depicted by larger circles. Here, an averaging interval of $t_{LS}=0.6$ is used. Considering the FAHM of the averaged RW elevation, the data are categorised into three K-means clusters in different colours. The linear fit and $95\,\%$ confidence interval are also given for each case, i.e. cases 1–6 in table 1.

Surprisingly, those nonlinear RW events with longer lifespans show greater FAHM levels, suggesting the occurrence of a long-term broadening of the directional wave spectrum. Here, we highlight the work by Toffoli et al. (Reference Toffoli, Onorato, Bitner-Gregersen and Monbaliu2010), which numerically observed the development of a similar bimodal pattern in the spectra. Meanwhile, the analysis by Osborne & Ponce de León (Reference Osborne and Ponce de León2017) showed that an initial and single JONSWAP spectrum could also become directionally unstable along the modulation channel. However, these strongly relevant studies were not directly related to the formation of full-spatial localised and directional RW structures. For further ease of the analysis, the RW events under each $\gamma$ case considered in figure 6 are categorised into three K-means clusters (Lloyd Reference Lloyd1982; Arthur et al. Reference Arthur2007; Cremonini et al. Reference Cremonini, De Leo, Stocchino and Besio2021), by using different colours, each representing a group of independent RW events with different FAHM and lifespan ranges. Note that the clustering is mainly used to group the many RW data observed from the simulations, and is achieved by adopting the K-means++ algorithm as developed by Arthur et al. (Reference Arthur2007), which is more effective than the standard K-means algorithm (Lloyd Reference Lloyd1982). As a result, the current K-means clustering in figure 6 derives three groups of RW events, each with distinct ranges of FAHM (two-dimensional bandwidth) and $t_{LS}$ . These can be approximately understood as linear, quasi-linear and nonlinear RW events. The specific $t_{LS}$ intervals of the clusters under different $\gamma$ values are summarised in table 2.

Table 2.

The $t_{LS}$ intervals of three clusters under different JONSWAP $\gamma$ values corresponding to figure 6.

Figure 7.

Averaged RW local elevation field, truncated according to the three K-means clusters from figure 6 (top three rows), and compared with the averaged RW elevation field calculated from all $t_{LS}=1T_p$ cases corresponding to one-RW events (fourth row), with all RW events as reference (fifth row), and the corresponding second-order NewWave theory spectrum (bottom row). Cases 4, 5 and 6 in table 1 are considered in the figure.

Upon averaging the RW events within each cluster and tracking the directional wave elevation field evolution, we can notice in figure 7 deviations between the long-living RW elevation fields and the NewWave theory (Taylor & Williams Reference Taylor and Williams2004) with lifetime $t_{LS}/T_p=1$ (bottom row), which increases with the increase of $t_{LS}$ or the cluster number. This is particularly the case for cluster 3, in which the elevation fields are nearly independent of the initial JONSWAP peakedness $\gamma$ values and form four dark wave trough holes around the centre peak. In fact, by averaging all RW events, the ‘all cluster’ results are almost consistent with the corresponding linear NewWave theory due to the dominance of short-living RW events. In such case, the role of long-living nonlinear RWs is completely overlooked.

In order to further understand and provide a more comprehensive picture of the observations above, figure 8 shows the three-dimensional shapes and characteristics of RW surface elevation fields as derived from the NewWave theory and numerical data clusters.

Figure 8.

Three-dimensional surface elevation corresponding to case 5 in table 1. Top left: the corresponding NewWave theory prediction; top right: the average of all RWs; bottom left: the shortest-living RWs with $t_{LS}/T_p=1$ only; bottom right: the longest-living RWs (i.e. cluster 3).

Whether such a coherent structure implies a possible deterministic description of so-called higher-dimensional nonlinear RWs, as is the case for NLS breathers for unidirectional wave states, needs future attention.

Figure 9.

Averaged RW local spectra corresponding to figure 7.

Notably, the overall three-dimensional shapes of these extreme waves is qualitatively similar to those already reported by Liu et al. (Reference Liu, Waseda, Yao and Zhang2022), yet the main focus of the latter work was not on the RW event longevity under variable JONSWAP peakedness parameters while varying crossing angles and wave steepness values. Considering the above, we next show the corresponding spectra of all cases discussed and shown in figure 9, confirming once again the significant and distinctive broadening of the directional spectrum from short- to long-lived $t_{LS}$ RW events. Note that with longer $t_{LS}$ , the directional spectrum gradually differs from the wave interference case, which is rather characterised by remaining narrowband in the two wave directions after the extreme wave focusing. Despite this, from both figures 7 and 9, one can notice that averaging all RW events can easily lead to the neglect of clusters 2 and 3, which are most affected by nonlinearity and directional spreading. Further, these clusters are characterised by the development of a dual bimodal structure in the spectrum. It is therefore recommended to consider and adopt an appropriate aggregated approach to analyse the real-world RW data (Häfner et al. Reference Häfner, Gemmrich and Jochum2021).

Figure 10.

Statistical PDF results of the correlation between the amplification factors, defined as $\eta _{{max}}/H_s$ , their lifespans $t_{LS}$ as RWs along the positive $x$ directions and the PDF while considering different crossing angles, i.e. cases 7, 8 and 9 from table 1.

Figure 11.

Averaged RW local elevation field, truncated according to the three K-means clusters in the same approach as for figure 7, and compared with the averaged RW elevation field calculated from all $t_{LS}=1T_p$ cases corresponding to one-RW events (fourth row), with all RW events as reference (fifth row) and the corresponding second-order NewWave theory spectrum (bottom row). Cases 7, 8 and 9 from table 1 are analysed in the figure.

In order to further validate our observation on such novel RW structure in crossing seas and avoid any potential wave breaking, we next investigate the effect of varying crossing angle on the formation of such RW structure. Figure 10 shows statistical PDF results of all RW events after lowering the wave steepness from the previous $k_pH_s=0.28$ to $k_pH_s=0.15$ and consider three different crossing angles $\beta$ = $28^\circ$ , $40^\circ$ and $53^\circ$ , i.e. cases 7, 8 and 9 in table 1. One can notice from the current statistical results that RW events have a significantly smaller amplification factor and $t_{LS}$ at $\beta$ = $28^\circ$ compared to the other two crossing angles, suggesting a much weakened nonlinear interaction between the two wave trains at a small crossing angle. This can be further confirmed in figure 11 showing the averaged directional RW elevation field. From figure 11 revealing the corresponding averaged RW elevation field, we can indeed notice that the change of the RW’s directional surface elevation with $t_{LS}$ is not significant when $\beta$ = $28^\circ$ compared with the other two angles that are within the most unstable region of $40^\circ {-} 60^\circ$ . The above observation is in qualitative agreement with the CNLS predictions (Onorato et al. Reference Onorato, Osborne and Serio2006a ) as well as with previous numerical and experimental evidence (Toffoli et al. Reference Toffoli, Onorato, Bitner-Gregersen and Monbaliu2010; Cavaleri et al. Reference Cavaleri, Bertotti, Torrisi, Bitner-Gregersen, Serio and Onorato2012).

Figure 12.

Averaged RW local spectra corresponding to figure 11.

Furthermore, the spectral evolution in figure 12, as well as that with fixed angle in figure 9, exhibits again a clear dual bimodal frequency evolution trend with a crossing angle $\beta$ within $40^\circ{-}60^\circ$ (Toffoli et al. Reference Toffoli, Onorato, Bitner-Gregersen and Monbaliu2010).

3.3. Comparison with the CNLS framework

To gain deeper insight into the role of nonlinearity, specifically the degree and order, in the manifestation of coherent, directional, large-amplitude waves in a crossing sea set-up, we compare the fully nonlinear results with simulations based on the CNLS framework.

Figure 13.

Comparison of directional and averaged RW local elevation field simulated by means of the fully nonlinear (ESBI) framework (left-hand panels) and the corresponding weakly nonlinear (CNLS) framework (right-hand panels), i.e. cases 10 and 11 from table 1.

Figure 14.

Comparison of the averaged RW local spectra corresponding to figure 13.

Figure 15.

Comparison of the magnified and directional-averaged localised maximal envelope peak shapes corresponding to figure 13.

After being classified into three K-means clusters, the directional and localised RW elevation results in figure 13 show a trend similar to those in figure 7 for a low significant wave height $k_pH_s=0.06$ within the JONSWAP wave field generated in each direction with $\gamma =9$ . Such wave height scales have been chosen to motivate future laboratory experiments.

Differences in averaged wave envelope shapes are noticeable in figure 13, particularly for the long-lived extreme waves, which belong to clusters 2 and 3. The directional spectra in figure 14 further confirm a typical and expected spectral broadening in both fully and weakly nonlinear frameworks. In fact, the CNLS predicts a broader spectrum and reduced dual bimodality trend, compared with the fully nonlinear ESBI simulations.

Intriguingly, one can also observe the differences in the extreme wave envelope peaks, as magnified in figure 15. While the bottom cluster 1 wave envelopes within the figure appear to be the result of wave overlap, showing a quasi-similar pattern, the long-lived extreme events in clusters 2 and 3 exhibit a completely different and distinct shape when simulated by the fully and weakly nonlinear frameworks. Such directional localised wave patterns at these different angles cannot be analytically modelled so far. Consequently, a complete and quantitative characterisation of these coherent extreme waves requires comprehensive future numerical and experimental explorations.