Lifetime characterisation of extreme wave localisations in crossing seas

Abstract

Rogue waves (RWs) can form on the ocean surface due to the well-known quasi-four-wave resonant interaction or superposition principle. The first is known as the nonlinear focusing mechanism and leads to an increased probability of RWs when unidirectionality and narrowband energy of the wave field are satisfied. This work delves into the dynamics of extreme wave focusing in crossing seas, revealing a distinct type of nonlinear RWs, characterised by a decisive longevity compared with those generated by the dispersive focusing (superposition) mechanism. In fact, through fully nonlinear hydrodynamic numerical simulations, we show that the interactions between two crossing unidirectional wave beams can trigger fully localised and robust development of RWs. These coherent structures, characterised by a typical spectral broadening then spreading in the form of dual bimodality and recurrent wave group focusing, not only defy the weakening expectation of quasi-four-wave resonant interaction in directionally spreading wave fields, but also differ from classical focusing mechanisms already mentioned. This has been determined following a rigorous lifespan-based statistical analysis of extreme wave events in our fully nonlinear simulations. Utilising the coupled nonlinear Schrödinger framework, we also show that such intrinsic focusing dynamics can be captured by weakly nonlinear wave evolution equations. This opens new research avenues for further explorations of these complex and intriguing wave phenomena in hydrodynamics as well as other nonlinear and dispersive multi-wave systems.

1. Introduction

Since the recording of the New Year or Draupner wave in 1995, fundamental research related to ocean rogue wave (RW) investigation has attracted much attention in recent decades due to its key relevance in coastal, ocean and arctic engineering applications (Kharif et al. 2008; Osborne 2010; Ducrozet et al. 2020; Mori et al. 2023; Klahn et al. 2024; Toffoli et al. 2024). Assuming the wave being unidirectional and ignoring the effect of bathymetry change (Li & Chabchoub 2023) and currents (Zheng et al. 2024), the formation of RWs can be explained as a result of wave superposition (Longuet-Higgins 1974; Fedele et al. 2016; McAllister et al. 2019; Häfner et al. 2021) or modulation instability (MI) (Benjamin & Feir 1967; Zakharov 1968; Tulin & Waseda 1999; Chabchoub et al. 2011; Bonnefoy et al. 2016). Both focusing mechanisms are equally important depending on the wave conditions at play (Dudley et al. 2019; Waseda 2020). The nonlinear mechanism in the form of MI, along with its manifestation in complex sea states (Tulin 1996; Waseda et al. 2009; Onorato et al. 2010; Gramstad et al. 2018; Toffoli et al. 2024), has been extensively studied as a key mechanism for wave group focusing. However, the predominance of quasi-four-wave resonant interactions for irregular ocean waves in crossing seas with strong directional spreading is considered to be less evident compared with unidirectional wave-field counterparts due to the violation of critical assumptions such as unidirectionality and narrowband spectral conditions (Janssen 2003; Mori et al. 2011; Fedele et al. 2016; Tang et al. 2021; Häfner et al. 2023). On the other hand, a recent experimental observation of nonlinear focusing dynamics in standing water waves (He et al. 2022b ) has shown that MI could still lead to notable amplifications in wave heights for such states. Here, standing waves can be indeed considered as a simplified specific type of crossing sea state with a crossing angle of $180^{\circ}$ . These findings further stress the involved nature of nonlinear wave interactions in complex configurations.

Furthermore and in contrast to the dispersive focusing mechanism in directional wave fields, numerical studies have postulated an increased likelihood of RW formation in coupled two-wave systems when considering weak nonlinearity in the modelling of crossing seas, with the limitations that both wave fields have the same peak frequency and are narrowband (Grönlund et al. 2009; Liu et al. 2022). Notably, Liu et al. (2022) made a successful attempt to study the crossing RW shape under varying crossing angles and spectral shapes. This approach goes beyond traditional RW investigations, which mainly focus on spectral evolution, exceedance probability distributions and kurtosis progression. The latter work also highlights that the shape of freak waves is more influenced by the crossing angle between wave components rather than the frequency or directional spectral bandwidth. To theoretically understand such processes, the coupled nonlinear Schrödinger equation (CNLS) (Okamura 1984; Onorato et al. 2006a ) and its higher-order forms (Gramstad & Trulsen 2011; Gramstad et al. 2018; Li 2021) are commonly used frameworks for the description of wave envelope interaction dynamics for crossing wave states (Cavaleri et al. 2012). Indeed, the directional (2D + 1) NLS and CNLS frameworks are essential in understanding the fundamental wave hydrodynamics together with the emergence of localised and directional wave patterns, such as RWs (Chabchoub et al. 2019; Steer et al. 2019b ; He et al. 2022b ). However, it still remains unknown whether the commonly used unidirectional and weakly nonlinear wave evolution equations can sufficiently predict the occurrence of all possible ocean RWs, which can occur also in crossing seas. For example, theoretical studies such as that of Guo et al. (2020) based on the Davey–Stewartson equation (Davey & Stewartson 1974) underline that directional perturbation can trigger strong localisations in time and directional space. Fully nonlinear numerical simulations also predict the spanwise wave instability with some success (Fructus et al. 2005b ). These extreme events cannot obviously be predicted by the classical unidirectional NLS framework with the famed breather solutions (Akhmediev et al. 1985; Peregrine 1983; Chabchoub et al. 2011; Tikan et al. 2022). On the other hand, the directional NLS can be modified to accommodate exact solutions of the unidirectional and integrable NLS (Saffman & Yuen 1978; Chabchoub et al. 2019; Waseda et al. 2021) and can describe the dynamics of modulationally unstable short-crested waves. This also implies that MI may play an important role in wave focusing processes beyond the unidirectionality assumption.

In fact, Onorato et al. (2006a )` discussed theoretically the role of MI in the wave focusing process in directional seas within the framework of the CNLS. Such a mechanism is active when the crossing angle between the two wave trains is below $70.53^\circ$ . For clarity, we define the wave crossing angle $\beta$ as the angle between the two crossing wave trains with wave direction vectors $(\kappa ,\iota )$ and $(\kappa ,-\iota )$ (Toffoli et al. 2011), i.e. $\beta =2\arctan (\iota /\kappa )$ . Numerical simulations solving the Euler equations using the high-order spectral method, an integration method originally proposed by West et al. (1987) and Dommermuth & Yue (1987), and physical wave tank experiments on crossing sea states (Toffoli et al. 2011) further confirmed that the maximum value of kurtosis is reached when the crossing angle is between $40^\circ$ and $60^\circ$ , which is in agreement with the weakly nonlinear theory (Onorato et al. 2006a ). The latter has also been applied by Cavaleri et al. (2012) to investigate the Louis Majesty cruise ship accident, which was likely caused by two almost identical wave states crossing at an angle of $50^\circ$ , as confirmed in the respective CNLS simulations. Furthermore, the occurrence of RWs within the crossing sea has been analysed by Bitner-Gregersen & Toffoli (2014) from hindcast data and validated by numerical high-order spectral method simulations, which shows that the maximum value of kurtosis is reached when the crossing angle $\beta =40^\circ$ , intriguingly, regardless of the wave directional spreading level (Bitner-Gregersen & Toffoli 2014). Independently of the valuable efforts mentioned above, which point to the experimental challenges in the comprehensive investigation of directional RWs, there have been limited attempts to classify directional RWs as emerged from nonlinear simulations. We are confident that our current study adds to the literature on bridging real ocean extreme waves with most recent RW theories, and future experimental studies with advanced observation techniques such as stereo imaging for ocean waves (Benetazzo et al. 2017; Guimarães et al. 2020) can further strengthen the impact of our findings.

In particular, to address the remaining key questions as discussed above, our numerical study, which is based on the fully nonlinear numerical framework developed by Wang et al. (2021), reveals the existence of a novel type of nonlinear, fully localised and directional RWs, i.e. extreme localised waves in directional space and time, which are distinct in their lifetime from the cases generated as a result of wave overlap. The procedure is initiated by accounting for the slow variation of the directional wave elevation field due to the effects of wave nonlinearity in the crossing sea state during wave data analysis, while the crossing NewWave theory (Taylor & Williams 2004) is adopted as the linear interference model dynamics for two generated and interfering JONSWAP sea states. In this context, we introduce a lifetime parameter, referred to as lifespan $t_{LS}$ of a RW event encompassing several consecutive RWs, which was adopted for similar purpose in previous works (Chabchoub et al. 2012; Kokorina & Slunyaev 2019; Slunyaev 2024) and as will be further elaborated in detail in the article. Based on the above, we reveal that fully localised RW elevations and their characteristic directional pattern are strongly correlated with their lifespans $t_{ LS}$ . The longer the $t_{ LS}$ of an extreme wave event, the more it differs from a large-amplitude wave created by the linear interference principle (Taylor & Williams 2004; Mathis et al. 2015; Birkholz et al. 2016; Fedele et al. 2016). At the same time, the wave energy exhibits a dual bimodal frequency evolution trend (Toffoli et al. 2010; Osborne & Ponce de León 2017) in the frequency domain. Such a key observation is also validated by considering a number of representative crossing angles (Liu et al. 2022), significant wave heights and JONSWAP peakedness factors. Finally, we identify a new type of fully localised RW envelope structure characterised by its unique long lifespan $t_{LS}$ and full-spatial localisation, i.e. in the $x$ direction, $y$ direction, time and considering the emergence of wave focusing recurrence. Despite having similar features to breathers, such dynamics cannot be predicted by the classical NLS or MI formalism that are applicable only for unidirectional waves and moderate carrier wave steepness.

The remainder of this paper is organised as follows: In § 2, the numerical methods adopted in the current investigation are introduced, including the fully nonlinear enhanced spectral boundary integral (ESBI) model and the weakly nonlinear CNLS framework. Section 3 provides a unique lifetime analysis of the fully and weakly nonlinear simulation data of crossing seas, including the spatio-temporal evolution of long-living RWs in such states and associated statistical analysis. Note that in § 3.1 and 3.2, only fully nonlinear numerical data are used for analysis, while in § 3.3 we compare the fully and weakly nonlinear numerical results to validate the key role of MI (Onorato et al. 2006b ) in a more subtle approach. The emerging directional structures of the long-living directional RW envelopes at the respective maxima are also revealed in this section using the CNLS framework. The results are summarised and discussed in § 4.

2. Numerical methods

Our study embraces fully and weakly nonlinear numerical schemes, which are described in detail below, together with details of the cross wave-field initialisation.

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 1968). They can be rewritten as a skew-symmetric form (Fructus et al. 2005a ):

(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. 2007; Wang & Ma 2015). 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 2025)

(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. 1987). 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. 2018, 2021; Wang 2025).

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. 2008; Gramstad et al. 2018; Mori et al. 2023). 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. 1973) 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. 2021), the selected peakedness values ensure that the spectrum widths in the current work are covered from broadband to narrowband conditions (Goda 2010). 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 2006; Janssen & Bidlot 2009) 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. (2021). To avoid the effects of potential wave breaking and high steepness values which can lower the accuracy of our model (Wang 2025), 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. 2011; Cavaleri et al. 2012; Bitner-Gregersen & Toffoli 2014). 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. 2005) 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. (2021) 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 2025), 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. 2022a ,c). 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. 2021). 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. 2013). 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 2025) 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. 2012), 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 1984; Onorato et al. 2006a ) 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 1968). 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 2010) 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. (2022b ) 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. Lifespan analysis of emerged RWs and categorisation

We begin by reasonably assuming that crossing RWs develop along the positive $x$ direction, i.e. at $\beta =0$ . During the prescribed evolution, when observing the wave elevation field, several localised peaks can be identified and tracked during their development before and after the steepest wave group focusing, and as figure 2 also indicates, a phenomenon also known as ‘three sisters’ (Kharif et al. 2008; Magnusson et al. 2019). This process can be referred to as a RW ‘event’ (Chabchoub et al. 2012; Sergeeva & Slunyaev 2013; Slunyaev et al. 2016; Kokorina & Slunyaev 2019; Slunyaev 2024). Considering the arguments above, we assume that crossing RWs within a same ‘event’ occur consecutively in both space (along the positive $x$ direction) and time (along the positive time axis). Thus, we trivially define that two RWs are considered part of one ‘event’ if their adjacent distance does not exceed two wavelengths in the positive $x$ direction and 0.8 wavelengths in the perpendicular direction. In addition to that, the latter RW must occur within a time interval not exceeding five wave periods. Note that sufficient redundancy has been reserved for the above spatial and temporal intervals, so that the lifespan may be slightly longer. That said, this hardly affects the upcoming ensemble analysis results, which are merely based on individual peak RWs of each event. Thus, the lifespan of an independent event $t_{LS}$ can be defined as:

(3.1) \begin{equation} t_{LS}=t_{RW,end}-t_{RW,ini}, \end{equation}

where $t_{RW,ini}$ and $t_{RW,end}$ are the occurrence time of the first and last identified RW in an identified extreme localisation event. In particular, we refer to the events with $t_{LS}/T_p=1$ as short-lived ones with a lifespan equal to or less than one peak wave period.

Figure 3. Comparison of two exemplary crossing RW events. Both events are extracted from case 3 in table 1 with different lifespans $t_{LS}$ , evolving from the top to bottom in each column. The averaged direction wave-field evolution is plotted every two peak wave periods $T_p$ . When $t/T_p=0$ , the current RWs reach their peak. The left-hand column shows a short-lifespan RW event with $t_{LS}=1T_p$ only. The right-hand column shows a long-lifespan RW event with $t_{LS}=40T_p$ . The orthogonal white dashed lines indicate the location of the compared RW events relative to the reference centre point.

According to the above-mentioned technique to group a series of RWs into RW events sharing the same lifespan, we identify and compare two representative RW events with different $t_{LS}$ in figure 3, which are extracted from the same random wave field, such as that computed and represented in figure 1. This suggests that RWs appearing at different time-scale lengths along the positive $x$ direction can be observed within the same crossing wave field. To simplify the further analysis, we treat such cases as independent events with the first class of these events being conjectured as a result of wave interference, due to the very short focusing time, as elaborated in § 3.2.

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. 2009), 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. 2010; Osborne & Ponce de León 2017), 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. (2010), which numerically observed the development of a similar bimodal pattern in the spectra. Meanwhile, the analysis by Osborne & Ponce de León (2017) 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 1982; Arthur et al. 2007; Cremonini et al. 2021), 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. (2007), which is more effective than the standard K-means algorithm (Lloyd 1982). 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 2004) 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. (2022), 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. 2021).

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. 2006a ) as well as with previous numerical and experimental evidence (Toffoli et al. 2010; Cavaleri et al. 2012).

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

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.

4. Discussion and conclusion

This paper investigates key physical and statistical properties of RW events in crossing JONSWAP-type sea states using a fully nonlinear ESBI framework (Wang et al. 2021) and a lifespan-based analysis approach (Chabchoub et al. 2012; Kokorina & Slunyaev 2019). Among the three selected crossing angles, we mainly focus on the most hazardous $\beta =40^\circ$ case and analyse RWs under different JONSWAP spectrum bandwidths and significant wave heights (Toffoli et al. 2011; Cavaleri et al. 2012; Bitner-Gregersen & Toffoli 2014).

Our fully nonlinear numerical ESBI results reveal RW events developing along the positive $x$ direction with lifespans ranging from $1T_p$ , satisfying the RW threshold criterion, up to $40 T_p$ and beyond. For these long-lasting nonlinear RW events, we observe a clear focusing and decay process of the averaged directional wave elevation field, along with a distinct recurrence of the averaged directional wave spectrum. This suggests that wave superposition and the classical (unidirectional) MI (Benjamin & Feir 1967; Mori et al. 2011) are insufficient for the prediction of all extreme wave events in such realistic directional seas.

On the other hand, we further analyse the statistical characteristics of such crossing RW events by varying the JONSWAP peakedness factor $\gamma$ and notice that the probability of long-lifespan events increases with increasing $\gamma$ , highlighting the potential role of quasi-four-wave resonant interactions in one of the two colliding wave beams, which might be not substantially influenced by the other wave-field component. Moreover, previous studies point to the significance of spreading effects (Toffoli et al. 2010) as well as directional weakly nonlinear effects (Osborne & Ponce de León 2017) in RW generation. The emergence of the localised severe and typical spectral broadening followed by a dual bimodality further supports the role of nonlinear wave interactions in the formation of RWs in cross seas.

Furthermore, classifying RW events into three clusters based on their corresponding lifespans reveals a gradual deviation of the three-dimensional elevation and energy spectra from the NewWave theory (superposition principle) with increasing longevity of the extreme wave events. Similar deviation and dual bimodality of the long-living RW events are also observed under a lower significant wave height and different crossing angles to exclude wave-breaking events, further supporting our observations above.

In order to explore the role of weakly nonlinear effects in the lifetime of RWs, we compare the fully nonlinear ESBI numerical results with the weakly nonlinear CNLS simulations. Both frameworks exhibit very good qualitative agreement, with consistent changes in wave elevation patterns and directional spectral broadening, particularly observed for longer extreme event lifespans. Interestingly, cluster 3, which regroups the long-lifespan extreme wave envelopes, reveals unique discrepancies in the RW coherence when modelled by ESBI and CNLS. This is likely due to the limitations of the CNLS approach (Liu et al. 2022) and of the Hilbert transformation for such wave systems.

Note that our study also differs from the experimental results reported by Steer et al. (2019a ). In the latter work, the authors superimpose two oblique and non-dispersive Gaussian wave envelopes to a unidirectional wave train, while our current work investigates two crossing random JONSWAP wave fields having a different wave vector. Moreover, the directional extreme events considered in our work do not remain stationary. Although Steer et al. (2019a ) and the current work both report ‘ ${\rm X}$ -like’ wave shape, the formation mechanisms and amplification factors are distinct. Furthermore, our study is in agreement with previous studies investigating random crossing seas (Toffoli et al. 2010; Cavaleri et al. 2012;Wang et al. 2021), suggesting that maximum kurtosis can be achieved for crossing angles ranging between $40^\circ$ and $60^\circ$ , and further revealed the distinguished crossing RW structures under both fully nonlinear ESBI and weakly nonlinear CNLS frameworks.

While predicting such extreme localisations is to date not analytically fully understood, the spatio-temporal localised and directional RW solutions (Qiu et al. 2016; Guo et al. 2020) offer promising avenues for future investigation.

In conclusion, this work unveils a characteristic type of directional nonlinear and coherent RW structure in crossing seas, highlighting the role of nonlinear wave interaction in the formation of extreme events in colliding two-wave systems. The detected RWs are characterised by short and long lifespans developing along the positive $x$ direction while the corresponding directional spectrum broadening is followed by a distinct dual bimodal pattern. We also confirm that CNLS is sufficient to qualitatively describe such long-lived freak waves for moderate steepness values. Further theoretical and experimental studies are necessary to fully comprehend and predict these nonlinear waves, beyond the limitations of the wave set-up and the CNLS framework as adopted in this work.