1 INTRODUCTION

Rogue, freak, or extreme waves are highly destructive ocean waves that represent serious threats to various marine activities, such as sea voyages, ocean fishing, and oil exploitation. Several physical mechanisms based both on linear and on nonlinear theories have been proposed to explain the formation of such waves (Kharif and Pelinovsky, 2003; Kharif et al., 2009). Comparing to linear theories, it has been suggested that nonlinear mechanisms related to second- and third-order nonlinear wave-wave interactions (Janssen, 2003; Fedele, 2008; Fedele and Tayfun, 2009) can evidently prove the prediction of rogue wave in both aspects of occurrence probability and wave heigh magnitude.

Nonlinear wave-energy focusing caused by such nonlinearities can lead to the formation of rogue waves more frequently and wave surface elevations that generally deviate from the Gaussian (normal) distribution, resulting in non-Gaussian sea states (Longuet-Higgins, 1963). Commonly used measures of the non-Gaussianity of sea state are skewness μ₃ and (excess) kurtosis μ₄:

where η denotes the wave surface elevation, and the terms in angled brackets denote statistical averages. With the consideration of nonlinearities up to third-order, it is clear that, skewness is contributed entirely by second-order nonlinear interactions involving bound waves (Tayfun, 1980; Tayfun and Fedele, 2007; Fedele and Tayfun, 2009; Janssen, 2009), and that kurtosis includes a dynamic component (μ₄^free) attributable to third-order quasi-resonant interactions (Janssen, 2003; Mori and Janssen, 2006) between free waves and another bound component (μ₄^bound) induced by both second- and third-order bound-wave nonlinearities (Tayfun, 1980; Tayfun and Lo, 1990; Fedele, 2008; Fedele and Tayfun, 2009; Janssen, 2009); these variables can be expressed as follows:

The non-Gaussianity of sea state is also sensitive to the geometries of the corresponding wave spectrum, and this relation has been well established using theoretical models (e.g., Janssen, 2003, 2009; Fedele and Tayfun, 2009) and confirmed by laboratory and numerical experiments (e.g., Onorato et al., 2009a, b; Toffoli et al., 2009; Waseda et al., 2009; Fedele, 2015). Generally, at least three geometric factors—wave steepness (hereafter, SP), bandwidth (BW), and directional spreading (DS)—can influence skewness, kurtosis, or both. For steeper SP and narrower BW and DS, it can be concluded that the deviation from Gaussianity will be greater, resulting in higher probability of rogue wave occurrence in such a sea state.

Following the above investigations, it has been established that two types of approach can be used to estimate skewness and kurtosis using the three spectral geometric factors. One is based on formulae/equations derived from theoretical wave models such as the Zakharov equation (Zakharov, 1968; Janssen, 2003) and the Tayfun-Lo model (Tayfun and Lo, 1990); and the other is based on phase-resolving models such as the high-order spectral method (HOSM) (Dommermuth and Yue, 1987; West et al., 1987) and the modified nonlinear Schrödinger equations (MNLS) (Dysthe, 1979; Trulsen and Dysthe, 1996).

Estimation of non-Gaussianity based on theoretical formulae and equations is highly convenient, and such estimates can certainly be applied to assess the non-Gaussianity associated with the evolution of sea state. Some methods have even been adopted in operational rogue wave forecasting systems, such as the ECMWF-IFS WAM (Janssen and Bidlot, 2009; ECMWF, 2016; Janssen, 2018) and the Space-Time Extremes (STE) forecasting system included in version 5.16 of WWⅢ (Barbariol et al., 2015; The WAVEWATCH Ⅲ Development Group, 2016; Barbariol et al., 2017). However, the indicators of skewness and kurtosis obtained using theoretical models, must first be obtained from theoretical models derived under the narrowband assumption in an environment with unidirectional wave propagation. When applying these indicators in scenarios with real sea conditions, specific parameters representing the spectral geometry are selected and undetermined parameters are introduced. Moreover, for the indicators adopted in rogue wave forecasting systems, calibrations are performed according to the final forecasting results; therefore, estimation of the indicators must balance the consideration of the spectral geometric factors and the other factors. Details regarding the process of calibrating operational non-Gaussianity indicators can be found in Barbariol et al.(2015, 2017) for the WWⅢ-STE system and in ECMWF Technical Memoranda (Janssen and Bidlot, 2009; Janssen, 2018) for the ECMWF-IFS system; the latter can also be found in the Appendix of this paper.

In comparison with theoretical models, phase-resolving models can be used to assess arbitrary two-dimensional (2D) spectra without any spectral shape limitations. These models adopt a given 2D spectrum as the initial wave field and simulate the evolution of surface elevations within a designed duration. Then, according to Eq.1, the corresponding skewness and kurtosis can be obtained based on the statistics of the surface elevations in the simulated wave field. The accuracy and feasibility of such models have been well verified by laboratory experiments (e.g., Toffoli et al., 2010). However, phase-resolving simulations are cumbersome, and performing such simulations based on wave spectrum time series obtained during sea state evolution is very time consuming. Therefore, the "phase-resolved" non-Gaussianity with the evolution of the sea state of interested can hardly be obtained and analyzed.

This paper presents a new approach that makes estimation of the phase-resolved non-Gaussianity more convenient. Using the estimation approach, phase-resolved non-Gaussianity can now be illustrated throughout the evolution of sea states of interest, not just at a few specific times. The remainder of this paper is organized as follows. The establishment of the estimation approach is introduced in Section 2. The application and validation of the approach based on certain sea states with rogue waves are described in Section 3. Finally, a discussion and the derived conclusions are presented in Sections 4 and 5, respectively.

2 ESTABLISHMENT OF THE ESTIMATION APPROACH 2.1 Experimental environment based on the HOSM

The HOSM numerical algorithm (Dommermuth and Yue, 1987; West et al., 1987) is dedicated to solve the Laplace equation for the velocity potential ϕ:

in a fluid, which is assuming to be inviscid and incompressible and the flow is irrotational. And surface elevations η can be obtained by solving Eq.3 at surface level z=η using the following equations:

where , and the ▽ operators in the equations above are 2D. Using the HOSM, ϕ can be expanded to a prescribed order M, and the complicated Dirichlet problem for ϕ at level z=η is transformed into M simpler Dirichlet problems for ϕ^(m) at z=0 (Ducrozet et al., 2016). The HOSM can obtain solutions for both ϕ and η with high accuracy and acceptable efficiency. Skewness and kurtosis can then be obtained using the solved η.

In this study, the HOSM experimental environment was established based on the open-source software package HOS-ocean ver. 1.5 (Ducrozet et al., 2016), which was developed at the Laboratoire de recherche en Hydrodynamique, Énergétique et Environnement Atmosphérique of the École Centrale de Nantes (France). HOS-ocean has been extensively validated in terms of nonlinear regular wave propagation (Bonnefoy et al., 2010), and it has also been widely adopted in previous research for modeling rogue wave sea state (e.g., Ducrozet et al., 2007; Ducrozet and Gouin, 2017; Jiang et al., 2019). HOS-Ocean can ensure the stability and convergence of the calculations. For example, all aliasing errors generated in the nonlinear terms can be removed. Time integration is performed with an efficient fourth-order Runge-Kutta Cash-Karp scheme, and the time step can be automatically selected to achieve a desired level of accuracy (or "tolerance", with typical values in the range of 10^-7–10^-5; in this study, the accuracy achieved was 10^-7). Moreover, a relaxation period of 10T_p (where T_p denotes the peak wave period) is often considered (as was the case in this study) to remove numerical instabilities attributable to fully nonlinear computations that start with linear initial conditions (Dommermuth, 2000). Additional details can be found in Ducrozet et al.(2016, 2017) and Ducrozet and Gouin (2017).

The HOSM experimental environment established in this study considers a physical space of size L_x=xlen×λ_p and L_y=ylen×λ_p, where xlen=ylen=25.5, λ_p is the wavelength at the peak frequency, and the discretization of the space is set as N_x×N_y=256×256. The HOSM is a pseudospectral method which means certain operators including the nonlinear terms are treated in the physical space instead of the spectral space, and the transforms between the physical and spectral spaces are handled efficiently using fast Fourier transforms. Thus, the spectral resolution of the pseudospectra can be determined as Δk_{x, y}=2π/L_{x, y}, and the spectral space extends from zero to k_max=(N_{x, y}–1)/2×Δk_{x, y}, where kmax is the cutoff frequency. Additionally, the HOSM cannot deal with wave breaking issues; therefore, a typical cutoff frequency of k_max=5k_p is generally adopted in simulations to obtain accurate solutions for the most energetic parts of the spectrum and to restrict the breaking of waves in the wave field to a very limited range (Ducrozet et al., 2017). For both the physical and the spectral spaces, the x-direction is taken as the principal direction of the pseudospectra, the y-direction is horizontal and perpendicular to the x-direction, and the z-direction in the physical space points vertically down with a consideration of infinite water depth.

The effects of spectral geometries on non-Gaussianity can be associated only with nonlinear wave-wave interactions, and therefore the HOSM simulations in this study considered only the wave-wave nonlinearities and ignored other factors that might influence the non-Gaussianity of the simulated wave fields, that is, wind force and energy dissipation due to breaking. According to current research, the non-Gaussianity of the sea state is mainly related to second- and third-order nonlinearities, nonlinearities up to the third order (i.e., M=3) were considered in the simulations.

2.2 Initial condition

To ensure that the initial conditions of the HOSM simulations were close to those of real sea states, the Joint North Sea Wave Project (JONSWAP) spectrum and specific directional spreading were introduced. The JONSWAP spectrum (Hasselmann et al., 1973; Goda, 1988, 1999) can be expressed as follows:

where H_s (m) is the significant wave height and f_p (Hz) is the peak frequency, which is also the reciprocal of the peak period T_p. Parameter σ=0.07 when f≤f_p and σ=0.09 when f > f_p. Parameter γ is the peak enhancement factor, which is related to the BW of Eq.6. Parameter β in Eq.6 can also be expressed in terms of γ as follows:

The adopted DS relations (Socquet-Juglard et al., 2005) were as follows:

where the spreading parameter Θ in degrees (°) or radians (rad) indicates that energy is distributed within the range of Θ/2 on both sides of the principal direction (0° or 0 rad). Finally, the 2D initial spectra in the HOS simulations are generated as follows:

The three spectral geometric factors in Eqs.6 & 8 are adjustable, and various initial conditions can be generated by considering different combinations of these factors. Based on Eq.6, the expression for SP is as follows:

where H_s is the significant wave height and λ_p is the wavelength at the peak frequency, which can be calculated from f_p according to the linear dispersion relation. Thus, SP is determined based on both H_s and f_p(T_p) in Eq.6. According to observations in the northern North Sea (deep water) from 1973 to 2001 (Haver, 2002), ε is generally < 0.06 and most frequently in the range of 0.01–0.02 (Fig. 1). To cover all observed sea conditions, the value of ε for the initial spectra was set to vary within the range of 0.01–0.06, and considering the number of calculations involved in the simulations, the interval was set to 0.01. Furthermore, without loss of generality, each value of ε in the range was associated with a unique combination of H_s and T_p. The settings of SP, together with the corresponding simulated physical and pseudospectral spaces, which are closely related to λ_p(T_p), are listed in Table 1.

The BW and DS values of the initial spectra can be determined from γ and Θ, respectively. For BW, parameter γ was set to vary within the range of 1.0–7.0 in accordance with the JONSWAP observations at an interval of 1.0. For DS, the range of the spreading parameter Θ was set to 8°–340°, and the interval of Θ in the range of 8°–176° (180°–340°) was set at 8° (20°).

Finally, the total number of initial spectra with adjustable ε, γ, and Θ values was 6×7×31=1 302, and the HOSM simulations were performed for each individual case. Regardless of how the combination of H_s and T_p might vary, the number of waves involved in each initial wave field will always be 25.5×25.5 (Section 2.1).

2.3 Skewness and kurtosis indicators

One complete HOSM simulation with a given initial spectrum is called one realization. Notably, previous experiments revealed that, with the nonlinearities that are up to third order, the most significant variations in skewness and kurtosis in the simulated wave fields all occurred in the first 180T_p of the simulation duration, and subsequent changes in both skewness and kurtosis were minimal because the simulated wave field tended to be Gaussian. Thus, the simulation duration of each realization was set to 180T_p. The skewness μ₃ and kurtosis μ₄ were calculated using Eq.1 based on the 256×256 surface elevations for each of the output simulated wave fields, and the output interval was 1T_p. Notably, the temporal evolutions of μ₃ and μ₄ obtained from one realization were random and unstable, and that obtaining averages of them from a number of realizations with the same initial spectrum would significantly reduce uncertainty, see Section 2.4.

We considered the indicators μ₃ and μ₄, which can represent the average characterization of nonGaussianity during each 180T_p simulation. Therefore, the temporal evolutions of μ₃ and μ₄ were first averaged for the final 170T_p (denoted and ), with the first 10T_p treated as the relaxation period, as noted in Section 2.1. Second, to compare the skewness and kurtosis obtained from different initial conditions, we introduced benchmarks denoted B_μ3 and B_μ4. The benchmarks B_μ3 and B_μ4 were the and results obtained from the initial condition of ε=0.06, γ=10, and Θ=4°, representing an extremely steep, near-narrowband, and near-unidirectional sea condition that would hardly ever be observed in a real ocean. Finally, we computed the ratios of the 'averages' to the 'benchmarks' as follows:

which could be used as non-Gaussianity indicators to identify the temporal evolution of skewness and kurtosis. Thus, the non-Gaussianity of the simulated wave fields corresponding to the initial spectra with certain combinations of ε, γ, and Θ was comparable.

2.4 Number of repetitions

Incorporation of a large number of repetitions in an averaging process generally produces results that are stable and convergent. Therefore, eight convergence tests were performed to determine the minimum number of repetitions to be averaged. Each test involved an individual initial condition, as listed in the right part of Fig. 2. For each initial condition, 100 realizations were conducted, of which n were then selected at random to produce a collection named C_n. For n=30, 40, 50, …, 100, there were eight collections to establish. The indicators R_μ3 and R_μ4 for each collection are shown in Fig. 2.

It can be seen from Fig. 2 that satisfactory convergence could be achieved for both R_μ3 and R_μ4 after 60–70 realizations. Accordingly, the number of repetitions involved in the averaging procedure was set to 80.

2.5 Spectral geometries and non-Gaussianity indicators

The obtained non-Gaussianity indicators constituted two reference sets for skewness and kurtosis, here denoted R_μ3 (ε, γ, Θ) and R_μ4 (ε, γ, Θ), respectively. Parts of the two references are illustrated in Fig. 3. As shown in Fig. 3a–d, the overall trends of R_μ3 and R_μ4 are affected by ε, γ, and Θ, confirming that as the value of ε increased, the value of γ increased, or as the value of Θ decreased, the value of R_μ3 and R_μ4 increased. Moreover, the values of R_μ3 and R_μ4 change continuously with the three geometric parameters. As illustrated in Fig. 3a–b, R_μ3 is affected most evidently by the SP parameter ε; moreover, Θ tends to cause slight increase in R_μ3 as it approaches 0°, as shown in the upper-left corner of Fig. 3a. Additionally, γ can influence R_μ3 when Θ is narrow, as shown in the upper-left corner of Fig. 3b. In Fig. 3c–d, within the range where Θ is extremely narrow (e.g., Θ≤20°), R_μ4 decreases markedly as Θ widens; when Θ is beyond the extremely narrow range, the value of R_μ4 decreases, but this decrease is much slower than that when Θ widens within the extremely narrow range. Parameters γ and ε can also have certain effects on R_μ4. For example, larger values of γ and ε result in larger values of R_μ4, and the corresponding increase becomes significant when Θ is extremely narrow. Same results can be found in the previous research (Gramstad and Trulsen, 2007), but in which the spectral geometries that can influence kurtosis were discussed in the form of spectral width in the direction of the main wave propagation (or alternatively the Benjamin-Feir Index) and in the direction perpendicular to the main wave propagation (or alternatively the crest length).

It should be noted that Rμ could be negative for some combinations of (ε, γ, Θ) owing to two possible factors. First, some of the simulated kurtosis values were rather small, and they oscillated slightly and randomly near 0, resulting in a negative average value. Second, μ₄^free might be influenced by defocusing of wave energy as the ratio of BW to DS increases (Fedele, 2015). Therefore, negative R_μ4 values represent sea states with inactive MI or the defocusing of energy, both factors might lead to decreased possibility of rogue waves forming in such sea states.

2.6 Development of the extraction approach

As R_μ3 and R_μ4 change continuously with ε, γ, and Θ, the non-Gaussianity indicators can be obtained via three-dimensional interpolation using an arbitrary combination of (ε_i, γ_i, Θ_i), where subscript i indicates an arbitrary position in each of the ranges of ε, γ, and Θ, respectively. The crucial step is to determine how the parameters ε_i, γ_i, and Θ_i that appear in the non-Gaussianity references might be related to the geometries of SP, BW, and DS in a given spectrum.

The SP for any spectrum can be calculated as follows:

where is the significant wave height based on the zeroth-order spectral moment m₀, and L_p is the wavelength at the peak frequency of the given spectrum. To avoid confusion, H_m0and L_p here denote wave characteristics obtained from a given spectra, while the H_s and λ_p in Eq.10 represent the same characteristics but calculated from the JONSWAP spectrum, i.e., Eq.6. Notably, S_p is always equal to ε_i for arbitrary spectral shapes.

For DS, the spreading parameter σ_θ (Kuik et al., 1988) can be adopted as follows:

By applying Eqs.8–13, we find that σ_θ increases monotonically with Θ within the range of Θ∈[8°, 340°], as indicated by the black solid line in Fig. 4a. This line can be fitted by a second-order polynomial, shown as the red dashed line in Fig. 4a, and expressed as follows:

After σ_θ has been calculated from the given spectrum, it is easy to solve for the root of Eq.14 in the range of 8°–340° to obtain the corresponding Θ_i.

For BW, the parameter Q_p (Goda, 1970) can be used:

Similarly, applying Eqs.6–15 reveals that Q_p increases monotonically with γ in the range of γ∈[1, 7], as indicated by the black solid line in Fig. 4b. This line can also be fitted by a second-order polynomial, shown as the red dashed line in Fig. 4b, and expressed as follows:

Thus, parameter γ_i can be obtained by solving for the root of Eq.16 in the range of (γ∈[1, 7]) with Q_p calculated from the given spectrum.

3 APPLICATION AND VALIDATION OF THE ESTIMATION APPROACH 3.1 Rogue wave events and wave modeling

The newly established approach is applied to some sea states with occurrence of rogue waves. The 10 rogue wave events discussed in this section were all observed using laser sensors installed on oil platforms in the North Sea. The locations of the platforms and UTC times at which the events occurred are listed in Table 2, together with the synchronously observed maximum wave heights (H_max) and significant wave heights (H_s), which were digitized from Magnusson and Donelan (2013) for Draupner and Andrea and from Guedes Soares et al. (2003) for the Alwyn events.

In this study, the wave fields containing the selected events were reproduced using the WWⅢ wave model ver.5.16 (The WAVEWATCH Ⅲ Development Group, 2016). The modeled spectral space was set with 36 directions at intervals of 10° and 35 frequencies spaced from 0.042 to 1.05 Hz as a geometric progression with a ratio of 1.1. The computational grids of the physical area modeled, as illustrated in Fig. 5, include an outer grid named NS4 (50°N–78°N, -18°E–15°E) with 0.25°×0.25° resolution and an inner grid named NS8 (52°N–68°N, -6°E–7°E) with 0.125°×0.125° resolution. Bathymetric data were obtained from ETOPO1 of the National Oceanic and Atmospheric Administration (NOAA) National Geophysical Data Center (NGDC) (National Centers for Environmental Information, 2017), and the wind force adopted in the simulations were derived from ECMWF-ERA5 reanalysis hourly data (Copernicus, 2018). As indicated by the Case IDs listed in Table 2, the modeling started 30 days before and ended 1 day after the occurrence of the Draupner and Andrea events, while for the Alwyn events, the modeling commenced 30 days before Alwyn_r1 and concluded 1 day after Alwyn_r8.

Some of the simulated bulk wave parameters are illustrated in Fig. 6a–d together with observed data digitized from published research for comparison. Comparison of the black solid lines denoting the digitized observations with the blue lines of the simulations reveals acceptable reproduction of the observed sea states; thus, these data provide a reliable foundation for the following analyses. The simulated significant wave height H_s, peak wave period T_p, peak wavenumber k_p, peak wavelength L_p, and parameter k_pd (k_p×d, where d is the water depth at the sites; Table 2) at the times of occurrence of all the rogue wave events are listed in Table 3.

3.2 Validation of the estimation approach

To validate the estimation approach, additional HOSM simulations were performed based on the same experimental environment as that established in Section 2.1. However, the initial conditions were replaced with modeled wave spectra for the 10 events, and simulations that considered the second-order nonlinearities (M=2) only and both the second- and third-order nonlinearities (M=3) were performed separately to elucidate the dominant mode in the wave fields. A similar approach was used in previous research (Fedele et al., 2016) in which the Draupner and Andrea events were studied.

The results of the additional simulations are illustrated in Fig. 7 and listed in Table 4. As shown in the upper-right corner of Fig. 7, each panel titled with a Case ID contains upper and lower parts exhibiting the temporal evolutions of skewness μ₃ and kurtosis μ₄ of the simulated wave field, respectively, and the x-axis in each panel represents the time duration with a dimensionless form of time (T_p). The blue and red lines shown in both parts of each panel depict the μ₃ and μ₄ simulated when considering M=2 and M=3, respectively. The corresponding non-Gaussian indicators (here, denoted R_μ3^rw and R_μ4^rw) obtained via the newly proposed estimation approach are presented by the black solid lines, and for comparison with the simulated blue and red lines, these two indicators have been multiplied by the benchmarks of B_μ3 and B_μ4, respectively (i.e., black lines in both parts of each panel represent the values of R_μ3^rw×B_μ3 and R_μ4^rw×B_μ4). For each event identified by a Case ID in Table 4, the non-Gaussianity indicators (R_μ3^rw×B_μ3, R_μ4^rw×B_μ4), mean values of μ₃ and μ₄ for M=2 and M=3 (, , , ), respectively, and the corresponding deviations (S_{μ3, M2}, S_{μ3, M3}, S_{μ4, M2}, S_{μ4, M3}) are all shown, and the deviations are defined as follows:

where N=170 denotes the number of μ₃ or μ₄ samples in the final 170T_p of each of the HOSM simulations (see Section 2.3).

As shown in the upper part of each panel in Fig. 7, the R_μ3^rw (×B_μ3) indicator agrees well with the skewness obtained from the simulated wave fields for both M=2 and M=3 nonlinearities, and in the lower part, as expected, R_μ3^rw(×B_μ3) is best matched to the kurtosis obtained from the simulated wave fields for M=3. Table 4 also shows that the values of R_μ3^rw×B_μ3, , and are similar for each event, whereas is generally larger than and R_μ4^rw×B_μ4 best matches , as illustrated by the lower values of S_{μ4, M3} in comparison with those of S_{μ4, M2}. Under the current experimental conditions, the acceptable goodness of fit of the cyan lines to the red lines in Fig. 7, together with the corresponding parameters illustrated in Table 4, can be considered to verify the feasibility of the newly developed approach in a real wave environment.

3.3 Evolution of non-Gaussianity in rogue wave sea states

Throughout the evolution of the simulated sea states, the indicators R_μ3^rw and R_μ4^rw estimated via the newly proposed approach are shown in Fig. 8, and the indicators for the Draupner (Fig. 8a and d), Andrea (Fig. 8b & e), and Alwyn (Fig. 8c & f) events are depicted by the magenta, blue, and red solid lines, respectively. In each panel of Fig. 8, the vertical dashed lines denote the times of rogue wave occurrence, and the horizontal lines represent (from top to bottom) the 90%, 75%, 50%, 25%, and 10% quantiles of each indicator (for R_μ3^rw the values are 0.485 5, 0.439 3, 0.362 5, 0.218 1, and 0.148 4, and for R_μ4^rw the values are 0.037 5, 0, 030 7, 0.021 9, 0.007 6, and 0.003 8) obtained within all the modeled durations at the three locations (about three months). The quantiles can help determine the level of the parameter and indicator in the evolution of sea state, and it is noted that the quantiles of the phase-resolved non-Gaussianity indicators now can be obtained because of the application of the newly proposed approach. The values of each indicator at the times of event occurrence are listed in Table 5.

As shown in Fig. 8a–c, the R_μ3^rw values at the times of event occurrence are reasonably similar and greater than the 75% quantiles (even over 90% at times during the Alwyn events), and such a situation can persist for several to dozens of hours both before and after event occurrence. For R_μ4^rw, Fig. 8d–f reveals that the kurtosis of the sea state does not exhibit any abnormality in comparison with that at the other times simulated, and the values of R_μ4^rw for the Draupner and Andrea events are even lower than the 50% quantile.

Thus, it can be concluded that the selected rogue wave events all occurred in non-Gaussian sea states with relatively large skewness, but normal kurtosis. Finally, it should be noted that high values of both R_μ3^rw and R_μ4^rw according to the quantiles could be found at times before the occurrence of the selected events, suggesting that the probability of rogue wave generation in such sea states was high. However, large values of R_μ3^rw and R_μ4^rw just indicate higher probability, rather than evidence of rogue wave occurrence; and on the other hand, rogue waves might have been generated in those sea states, but have not been observed.

4 DISCUSSION

To assess the differences in non-Gaussianity indicators estimated using theoretical formulas and the newly proposed approach, several operational indicators adopted in the extreme wave forecasting system of the Integrated Forecasting System operated by European Centre for Medium-Range Weather Forecasts (ECMWF-IFS) are adopted here as an example. The meaning of each indicator and the methods used for calculating them are given in the Appendix. The evolution trends of the operational indicators are illustrated in Fig. 9, where the colored, horizontal, and vertical lines depict the same events, quantiles, and times just the same as Fig. 8.

Figure 9a–c shows that indicator C₃, representing skewness, exhibits higher values (greater than the 75% quantile) at the time of event occurrence than at other times; moreover, such a situation also persist for several to dozens of hours around the occurrence of the events, just confirming the R_μ3^rw results in Fig. 8. However, different from the normal values of kurtosis indicator R_μ4^rw found at the times of rogue wave occurrence, the kurtosis indicator C₄ (Fig. 9d–f) displays abnormal behavior (even greater that the 90% quantile) in comparison with that of the kurtosis obtained at other moments. Furthermore, as can be observed in the lower panels of Fig. 9, indicators C₄^bound (Fig. 9g–i) and C₄^free (Fig. 9j–l), which are associated with the bound and free parts of kurtosis, respectively, both exhibit extremely high values at the times of event occurrence; nevertheless, the contribution of C₄^free to parameter C₄ was much smaller than the contribution of C₄^bound, and the abnormality of C₄ at the times of rogue event occurrence was mainly influenced by C₄^bound.

As expressed in Eqs.2 & A2 (from Appendix), kurtosis comprises both a dynamic (free) part and a bound part. The dynamic contribution is relatively small; for example, it can be < 10% of the bound component in a normal sea state with broad BW and DS (Annenkov and Shrira, 2014); however, it can reach high levels in specific wave environments, suggesting that MI are active (Janssen, 2003; Fedele, 2015). The C₄^free at the occurrences of the events are found being > 20% of the C₄^bound, which is slightly larger than a "normal" value expected. However, the lack of notable discrepancies between the black (M=2) and red (M=3) lines in the lower part of each panel of Fig. 7 indicates that the kurtosis of the selected rogue wave sea states is dominated by second-order nonlinearities, thus, the third-order nonlinearities and the associated MI are inactive; similar conclusions can also be found in Fedele et al. (2016). As mentioned previously, the operational indicators derived from the theoretical formulae might not be suitable for real sea states, thus, the prominent C₄ values found at the times of rogue wave occurrences might be caused by an improper calibration. It should be noted that the estimation results obtained with the phase-resolving model, including the R_μ3 and R_μ4 values obtained in this study, might also be incomplete. It is difficult to discuss which kind of estimation is correct or not with the current results, and that is not the main purpose of this study.

The non-Gaussianity indicators obtained in different ways can provide different perspectives for studying rogue wave sea states. And it is not the first time that a phase-resolving model, such as HOSM, has been applied to assess rogue sea states (e.g., Trulsen et al., 2015; Fedele et al., 2016; Jiang et al., 2019). However, phase-resolving simulations can be performed only using several spectra selected from the full spectral evolution of the sea state because full simulations can be cumbersome. In this study, using the precalculated dataset and the extraction approach, phase-resolved non-Gaussianity indicators can be obtained conveniently. Thus, the "normal" or "abnormal" behaviors of sea state non-Gaussianity can be identified according to the statistics of the phase-resolved indicators, which are derived throughout the simulated sea states. Therefore, this study provides a new perspective for studying non-Gaussian sea states and the associated rogue wave sea states, and such studies can be performed based on the phase-resolved non-Gaussianity throughout the evolution of sea states of interest, not just at a few specific times.

With respect to real ocean conditions, our model is certainly an oversimplification; for example, it focuses purely on the nonlinearities among waves, ignoring other physical mechanisms that might influence non-Gaussianity, e.g., wind forcing, wave breaking. Furthermore, the interaction between two wave systems were also ignored, regarding that for two systems do not interact much, the combined wave system appears to be more Gaussian than the corresponding partitioned wave systems (Støle-Hentschel et al., 2020). Therefore, the model was established based only on unimodal spectral shapes and bimodal shapes were ignored. And in the application, the quantiles are just obtained based on samples simulated in a few months, involving more samples might provide more suitable criterion. Further improvement of the new estimation approach can be done in future studies.

5 CONCLUSION

In this study, we established a systematic approach for estimating the evolution of sea state non-Gaussianity. The newly established approach includes: i) a set of precalculated indicators representing the relation between sea state non-Gaussianity and various combinations of the SP, BW, and DS geometries, and ii) an approach that can extract the skewness and kurtosis indicators from the precalculated dataset for given 2D spectra. Because the precalculated dataset was obtained based on HOSM (phase-resolving) simulations, the indicators can be applied to real sea states without calibration of spectral shapes. With the newly developed extraction approach, estimating the phase-resolved non-Gaussianity of a given sea state becomes convenient, and the phase-resolved non-Gaussianity can be estimated with the evolution of sea state. Validation of the estimation approach was performed based on analysis of sea states in which real rogue waves occurred. The sea states were reproduced using the spectral wave model WWⅢ, and additional HOSM simulations were conducted with simulated wave spectra at the time of occurrence of the rogue wave events. The acceptable goodness of fit of the extracted indicators to the simulated results verified the feasibility of the estimation approach for use in real wave environments.

With the estimation approach, the evolution of the phase-resolved non-Gaussianity in the selected sea states was illustrated. It is found that the selected rogue wave events all occurred in non-Gaussian sea states with relatively large skewness, but normal kurtosis; and that large values of both the skewness and kurtosis indicators can be found at times before the occurrence of the selected events, suggesting high probability of rogue wave generation in such sea states. And it should be noted that the "large" and "normal" judgements given above are based on the statistics of the phase-resolved non-Gaussiniaty indicators, which are obtained throughout the evolution of the simulated sea state.

In comparison with the evolution trends of certain operational indicators derived from theoretical formulae, different behaviors of non-Gaussianity of the same selected the sea states are found. Although, the results of the newly proposed approach seem more reasonable, it is difficult to discuss which kind of estimation is better; and the newly proposed approach still has defects. However, this study provides a new perspective for studying the non-Gaussian rogue wave sea states, and further research could be performed on this basis.

6 DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available on request from the first author Xingjie JIANG (jiangxj@fio.org.cn) or the corresponding author Tingting ZHANG (zhangtt@fio.org.cn).

7 ACKNOWLEDGMENT

We thank Guillaume Ducrozet and Yves Perignon from the Research Laboratory in Hydrodynamics, Energetics and Atmospheric Environment (LHEEA) of the École Centrale de Nantes and The French National Centre for Scientific Research (CNRS) for their assistance in helping us understand the HOSM method and the use of HOS-ocean.

Electronic supplementary material

Supplementary material (Appendix) is available in the online version of this article at https://doi.org/10.1007/s00343-021-1236-1.