<figure><img src="image_1.png" /></figure>

HAL
open science

# Statistics of long-crested extreme waves in single and mixed sea states

Lei Wang, Jinxuan Li, Shuxue Liu, Guillaume Ducrozet

► To cite this version:

Lei Wang, Jinxuan Li, Shuxue Liu, Guillaume Ducrozet. Statistics of long-crested extreme waves in single and mixed sea states. Ocean Dynamics, 2021, 71 (1), pp.21-42. 10.1007/s10236-020-01418-9.
hal-03124379

HAL Id: hal-03124379

https://hal.science/hal-03124379

Submitted on 11 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Statistics of long-crested extreme waves
in single and mixed sea states

Lei Wang<fnref n="1" />,<fnref n="2" />, Jinxuan Li<fnref n="1" />, Shuxue Liu<fnref n="1" />, Guillaume Ducrozet<fnref n="2" />

Most of the studies on extreme waves are focused on the systems with single-peak wave spectra. However, according to the statistics of occurrence, the bimodal spectral system is also frequent in real sea conditions. In order to summarize the statistics of extreme waves, irregular wave trains under single-peak and bimodal spectra for long durations are simulated in this paper, based on a two-dimensional High Order Spectral (HOS) numerical wave tank. A large number of configurations have been tested under unimodal and bimodal spectra. The investigation on the wave trains under single-peak spectrum indicates that although in conditions often referred as deep water ($k_p h > \pi$), the relative water depth has a significant influence on the probabilities of occurrence of extreme waves. A detailed analysis of the combined effect of Benjamin-Feir Index (BFI) and relative water depth is provided. However, the situation is more complex in real sea conditions, which may exhibit multimodal spectra. We focus in this study on long-crested bimodal spectra characterized by the same significant wave height $H_s$ and mean zero-crossing period $T_z$ of the sea states as the single-peak spectrum. The wave conditions under bimodal spectrum present milder extreme wave statistics than those under single-peak spectrum. In addition, mixed ocean systems with equivalent energy distribution (i.e., *Sea-Swell Energy Ratio* (SSER) is close to 1.0) and larger separation between partitions (i.e., *Intermodal Distance* (ID) > 0.10) are the less prominent to extreme waves appearance. The comparison of the mixed sea states and the corresponding single independent systems demonstrates that the complexity of the underlying physics of a given sea state (for instance the presence of modulational instability or other nonlinear process) cannot be deduced by an analysis limited to the statistical content of the combined sea state. The wave energy being distributed among frequencies plays a major role. Additionally, Gram-Charlier distribution can accurately predict the probability of large waves ($1.5 < H/H_s < 2.0$) compared to the MER distribution, but it underestimates the statistics of the wave height distribution when $H/H_s$ is larger than 2.0 for both single-peak and bimodal states.

**Keywords** Bimodal spectrum, Long-crested extreme waves, Kurtosis, Wave height distribution, HOS numerical wave tank

## 1 Introduction

As frequent human activities (such as the exploration and development of marine resources) move towards the deep ocean, the probability to encounter extreme waves (also known as rogue or freak waves) is becoming larger. The analysis and understanding of the formation of such events is of

great practical significance and has attracted more and more attention recently. Extensive studies have been performed in the past decades, focused mainly on the evaluation of the probability distribution of the wave elevation and height and hence to investigate the statistics of extreme waves in a large number of random waves from both theoretical and applied points of view (reviewed by Forristall 1984, 2000; Tayfun 1990, 2006).

Great progress in probabilistic characteristics of wave trains under single-peak spectrum have been made. In a stationary, Gaussian and extremely narrow-banded process, the wave height may be regarded as twice the envelop amplitude and that these are distributed according to a Rayleigh probability distribution (Longuet-Higgins 1952). Naess (1985) derived a model for the crest-to-trough wave height in a wide-band process, which has received wide popularity. In practice, the behavior of real sea waves is usually nonlinear and broad

<figure><img src="image_3.png" /></figure>

<footnote>Jinxuan Li
lijx@dlut.edu.cn</footnote>

<footnote n="1">¹ State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, Dalian 116024, China</footnote>

<footnote n="2">² Ecole Centrale de Nantes, LHEEA Lab (ECN and CNRS), Nantes, France</footnote>
banded. Thus, the linear model becomes invalid, especially for large wave height distribution (Tayfun 1981; Stansell 2004). Several studies based on theoretical developments as well as physical experiments have shown that the third-order interactions, also known as modulation instability, plays an important role in the appearance of large steepness waves (Janssen 2003; Onorato et al. 2004; Mori and Janssen 2006). Janssen (2003) has proposed a parameter known as the Benjamin-Feir Index (BFI), which is the ratio between the wave steepness $\epsilon$ and the relative spectral bandwidth $\Delta f/f_p$, as a measure of the strength of those third-order interactions in random wave trains. Onorato et al. (2004) have investigated the evolution of irregular sea states considering various BFI values in a large wave flume. Their results have shown that for a small value of BFI, the Rayleigh distribution predicts the experimental data rather well. However, for a large value of BFI, the Rayleigh distribution underestimates the tail of the probability density functions of free surface elevation and wave heights.

The third-order interactions can change the statistical properties of surface waves under specific conditions. These effects are mainly represented by the variation of the statistical parameter known as the kurtosis, the fourth-order moment of the probability density function of the free surface elevation. The kurtosis is considered a statistical parameter that can give an indication of the presence of extreme events in the time series, being a measure of the importance of the tail of the distributions (Janssen 2003). Mori and Yasuda (2002), followed by Mori and Janssen (2006), have discussed the formal relation between the kurtosis and the exceedance probability for wave heights. The kurtosis enters the distribution function as a correction to the Rayleigh distribution. When the kurtosis is 3.0, the distribution function is close to the Rayleigh distribution for narrow-banded wave fields. Considering the effect of the kurtosis, Mori and Janssen (2006) developed a distribution called a modified Edgeworth-Rayleigh distribution (MER) under the assumption of weak nonlinearity, narrow spectrum, and wave height twice of wave amplitude. Mori et al. (2007) have confirmed that the nonlinear correction to the maximum wave height depends on the kurtosis. Moreover, the tail of the wave height distribution increases along with the increase of the kurtosis.

Most of the previous studies are limited to sea states defined with a unimodal spectrum. However, 15–25% of the ocean state at different locations around the world is observed as a bimodal spectrum with wind-sea and the swell systems mixed (Guedes Soares 1984, 1991). Several studies have been performed on the two coexisting systems in crossing sea states. Onorato et al. (2006) have theoretically investigated the crossing seas through coupled nonlinear Schrödinger (CNLS) equations, and have found larger regions of instability and larger growth rate for coupled systems compared to single systems. Regev et al. (2008) have pointed out that a wind-sea and a weak swell at right angle can give rise to wave modulations, leading to increased

occurrence of freak waves. Similarly, Gramstad and Trulsen (2010) have computed the modification of the probability of freak waves in a wind-sea perturbed by a weak swell oriented at various angles to the wind-sea. They have found that the swell can enhance the occurrence of freak waves in the wind-sea slightly. Based on the investigation of Toffoli et al. (2011), the kurtosis is concluded to be increased for two systems with an interaction angle between 40° and 60°.

Contrary to the investigation in crossing sea states, the statistics of the unidirectional co-propagating mixed sea states have not yet been fully summarized. Some studies exist with sea states of equal significant wave height (Petrova and Guedes Soares 2009, 2011) or with the same wind-sea spectrum (Støle-Hentschel et al. 2020). Petrova and Guedes Soares (2009, 2011) only investigated the effects of the relative energy ratio (not considering the peak frequency separation between the low- and high-frequency wave fields on the wave height distribution) and had found that for swell-dominated seas, the kurtosis is lower than that in wind-sea-dominated seas. Støle-Hentschel et al. (2020) had noted that the mixed sea has milder extreme wave statistics than the pure wind-sea. They suggest the necessity to partition wind-sea and swell to better interpret the combined extreme wave statistics.

Real sea states are known to exhibit directional properties that have a significant influence on the wave statistics (Onorato et al. 2009). However, it is common practice in ocean engineering to consider only long-crested seas. In the design of structures at sea, this simplification enables the use of simplified numerical procedure in the evaluation of wave loads for instance. As a result, the unidirectional configuration is still the subject of many investigations, including the present one. In addition, the definition of the metocean conditions often describe the energy content of the sea states by simple integral parameters such as the significant wave height ($H_s$) or the mean zero-crossing period ($T_z$). The possible limitations inherent to the latter choice are investigated through the use of bimodal spectra (compared to unimodal ones).

However, the extensive study of the extreme wave occurrence in this simple physical context is not yet available. The objective of the present work is consequently to bridge this gap in order to provide more comprehensive understanding of the extreme wave statistics in different ocean environments. The first part introduces the different sea states studied during the analysis. It is followed by a second section dedicated to the presentation of the numerical wave tank and data analysis method at use in the study. Then, unidirectional single-peak spectra are investigated in detail, summarizing the main effects of BFI and relative water depth on the statistical properties of sea states. Finally, the case of the bimodal seas is studied extensively, assuming the wave field is described with a given set of integral parameters (the significant wave height $H_s$ and the mean zero-crossing period of the sea states $T_z$) but different values of *Sea-Swell Energy Ratio* (SSER) and *Intermodal Distance* (ID).
# 2 Irregular sea state conditions

## 2.1 Bimodal wave spectrum

In the following analysis, for the frequency spectrum $S(\omega)$, the six-parameter spectrum proposed by Ochi and Hubble (1976) is employed to represent both the single-peak spectrum and the bimodal spectrum:

$$S(\omega) = \frac{1}{4} \sum_j \frac{\left(\frac{4\lambda_j+1}{4}\right)^{\lambda_j}}{\Gamma(\lambda_j)} \left(\frac{\omega_{pj}}{\omega}\right)^{4\lambda_j} \frac{H_{sj}^2}{\omega} \exp\left[-\frac{4\lambda_j+1}{4} \left(\frac{\omega_{pj}}{\omega}\right)^4\right] \quad j=1,2 \quad (1)$$

$\Gamma$ is the gamma function, and $j=1, 2$ represent the low-frequency and high-frequency parts, respectively, which in turn correspond to the swell and wind-sea in the mixed sea system, respectively. Each part has three parameters, respectively, viz. the significant wave height $H_s$, the spectral peak angular frequency $\omega_p$ (=2$\pi f_p$), and the peak enhancement factor $\lambda$.

According to the Rice theory (Rice 1944), the significant wave height $H_s$ of the mixed waves can be determined by the method of the energy superposition

$$H_s = \sqrt{H_{s1}^2 + H_{s2}^2} \quad (2)$$

where $H_{s1}$ and $H_{s2}$ stand for the significant wave heights of the low-frequency and high-frequency parts, respectively.

## 2.2 Relative spectral bandwidth

Irrespective of the peak frequency, the value of the relative spectral bandwidth $\Delta f/f_p$ (here $\Delta f$ is defined as half-width at half-spectrum height) becomes static once the value of the peak enhancement factor $\lambda$ in the Ochi-Hubble spectrum is defined. For reference, the relationship between the enhancement factor $\gamma$ in the widely used single-peak JONSWAP spectrum and the present parameter $\lambda$ in bimodal spectrum is calculated and listed in Table 1.

## 2.3 Definition of Benjamin-Feir Index

There are different possible definitions for Benjamin-Feir Index (BFI), depending on the context of the study. Here, BFI is defined as the ratio between the wave steepness $\epsilon = k_p H_s / 2$ and the relative spectral bandwidth $\Delta f/f_p$ for deep-water waves (Onorato et al. 2004):

$$\text{BFI} = \frac{\sqrt{2\epsilon}}{2\Delta f/f_p} \quad (3)$$

Complementary, the extension to shallower water depth can be achieved by means of a modified Benjamin Feir Index, namely $B_S$, expressed as a

function of the deep-water BFI (Janssen and Onorato 2007; Janssen and Bidlot 2009):

$$B_S^2 = -\text{BFI}^2 \times \left(\frac{c_g}{c_0}\right)^2 \frac{gX_{nl}}{k_p \omega_p \omega_p''} \quad (4)$$

where $c_g$ and $c_0$ are the group and phase velocities, respectively. $\omega_p''$ is the second derivative of angular frequency $\omega_p$ with respect to wavenumber $k_p$. $X_{nl}$ is interaction coefficient. The specific expressions can be referenced to Janssen and Onorato (2007) and Janssen and Bidlot (2009). Note that in the real deep water (i.e., $k_p h$ tends to infinity), the second part of the expression for $B_S^2$ becomes -1 and the modified Benjamin Feir Index $B_S$ reduces to the deep-water definition BFI.

## 2.4 Definition of the different sea states

To investigate the statistical properties of the irregular wave trains and enhance the understanding of the extreme wave occurrence, different sea states under single-peak and bimodal spectra are simulated and analyzed in this study. The different conditions tested are detailed hereafter.

### 2.4.1 Single-peak spectrum

A significant number of studies have been dedicated to characterizing the features of the wave trains under single-peak spectrum. As an example, nonlinear simulations, experiments, and theoretical works have been conducted to try to predict the evolution of kurtosis in different contexts (Onorato et al. 2004, 2005, 2006; Toffoli et al. 2009; Shemer and Sergeeva 2009; Fernandez et al. 2014, 2016; Janssen 2003; Janssen and Bidlot 2009; Fedele 2015). However, as far as we know, existing works in a fully nonlinear context only consider one influencing factor, either BFI (Onorato et al. 2004, 2005, 2006; Shemer and Sergeeva 2009) or relative water depth $k_p h$ (Toffoli et al. 2009; Fernandez et al. 2014, 2016), and exhibit a limited range of variations. There is consequently a need in providing a thorough study of the effects of both parameters over a large range for unidirectional irregular sea states.

Five cases listed in Table 2 are designed to study the occurrence probability of the irregular wave trains under single-peak spectrum. Among them, Case 1, Case 2, and Case 3 are with the same peak enhancement factor but different peak frequencies (i.e., relative water depth $k_p h$), whereas Case 2, Case 4, and Case 5 have the same peak frequency but different peak enhancement factors. In each case, there are five different values of BFI, viz. 0.2, 0.4, 0.6, 0.8, and 1.0. Detailed parameters of various cases under the single-peak spectrum are listed in Table 2.

As an example, Fig. 1 displays the input amplitude spectra of different cases listed in Table 2. Figure 1a compares the
<table><caption>**Table 1** The relationship between $\gamma$ in JONSWAP spectrum and $\lambda$ in bimodal spectrum</caption><thead><tr><th>Bimodal spectrum<br>λ</th><th>JONSWAP<br>γ</th><th>Δf/f<sub>p</sub></th><th>Bimodal spectrum<br>λ</th><th>JONSWAP<br>γ</th><th>Δf/f<sub>p</sub></th></tr></thead><tbody><tr><td>1.0</td><td>1.0</td><td>0.2900</td><td>9.0</td><td>3.1</td><td>0.0983</td></tr><tr><td>1.5</td><td>1.2</td><td>0.2375</td><td>9.3</td><td>3.3</td><td>0.0967</td></tr><tr><td>2.0</td><td>1.4</td><td>0.2054</td><td>9.5</td><td>3.4</td><td>0.0958</td></tr><tr><td>2.5</td><td>1.5</td><td>0.1825</td><td>10.0</td><td>3.5</td><td>0.0933</td></tr><tr><td>3.0</td><td>1.6</td><td>0.1683</td><td>10.5</td><td>3.6</td><td>0.0917</td></tr><tr><td>3.5</td><td>1.7</td><td>0.1558</td><td>11.0</td><td>3.9</td><td>0.0892</td></tr><tr><td>4.0</td><td>1.8</td><td>0.1475</td><td>11.5</td><td>4.1</td><td>0.0867</td></tr><tr><td>4.5</td><td>1.9</td><td>0.1383</td><td>12.0</td><td>4.3</td><td>0.0850</td></tr><tr><td>5.0</td><td>2.1</td><td>0.1308</td><td>12.5</td><td>4.4</td><td>0.0842</td></tr><tr><td>5.5</td><td>2.2</td><td>0.1258</td><td>13.0</td><td>4.6</td><td>0.0825</td></tr><tr><td>6.0</td><td>2.3</td><td>0.1200</td><td>13.5</td><td>4.9</td><td>0.0800</td></tr><tr><td>6.5</td><td>2.4</td><td>0.1150</td><td>14.0</td><td>5.1</td><td>0.0792</td></tr><tr><td>7.0</td><td>2.5</td><td>0.1108</td><td>14.5</td><td>5.6</td><td>0.0768</td></tr><tr><td>7.5</td><td>2.7</td><td>0.1075</td><td>15.0</td><td>5.7</td><td>0.0767</td></tr><tr><td>8.0</td><td>2.9</td><td>0.1033</td><td>15.5</td><td>6.3</td><td>0.0737</td></tr><tr><td>8.5</td><td>3.0</td><td>0.1008</td><td>16.0</td><td>7.0</td><td>0.0712</td></tr></tbody></table>

<table><caption>**Table 2** Detailed parameters of various cases under single-peak spectrum</caption><thead><tr><th>Case</th><th></th><th>f<sub>p</sub><br>(Hz)</th><th>k<sub>p</sub><br>(m<sup>-1</sup>)</th><th>h/L<sub>p</sub></th><th>k<sub>p</sub>h</th><th>H<sub>s</sub><br>(m)</th><th>ε = k<sub>p</sub>H<sub>s</sub>/2</th><th>λ</th><th>Δf/f<sub>p</sub></th><th>BFI</th></tr></thead><tbody><tr><td rowspan="5">Case 1</td><td>Case 1_a</td><td rowspan="5">0.6</td><td rowspan="5">1.44</td><td rowspan="5">0.92</td><td rowspan="5">5.80</td><td>0.0278</td><td>0.0201</td><td rowspan="5">16</td><td rowspan="5">0.071</td><td>0.2</td></tr><tr><td>Case 1_b</td><td>0.0556</td><td>0.0403</td><td>0.4</td></tr><tr><td>Case 1_c</td><td>0.0834</td><td>0.0604</td><td>0.6</td></tr><tr><td>Case 1_d</td><td>0.1112</td><td>0.0805</td><td>0.8</td></tr><tr><td>Case 1_e</td><td>0.1390</td><td>0.1007</td><td>1.0</td></tr><tr><td rowspan="5">Case 2</td><td>Case 2_a</td><td rowspan="5">0.8</td><td rowspan="5">2.58</td><td rowspan="5">1.64</td><td rowspan="5">10.30</td><td>0.0156</td><td>0.0201</td><td rowspan="5">16</td><td rowspan="5">0.071</td><td>0.2</td></tr><tr><td>Case 2_b</td><td>0.0313</td><td>0.0403</td><td>0.4</td></tr><tr><td>Case 2_c</td><td>0.0469</td><td>0.0604</td><td>0.6</td></tr><tr><td>Case 2_d</td><td>0.0625</td><td>0.0805</td><td>0.8</td></tr><tr><td>Case 2_e</td><td>0.0782</td><td>0.1007</td><td>1.0</td></tr><tr><td rowspan="5">Case 3</td><td>Case 3_a</td><td rowspan="5">1.0</td><td rowspan="5">4.02</td><td rowspan="5">2.56</td><td rowspan="5">16.10</td><td>0.010</td><td>0.0201</td><td rowspan="5">16</td><td rowspan="5">0.071</td><td>0.2</td></tr><tr><td>Case 3_b</td><td>0.020</td><td>0.0403</td><td>0.4</td></tr><tr><td>Case 3_c</td><td>0.030</td><td>0.0604</td><td>0.6</td></tr><tr><td>Case 3_d</td><td>0.040</td><td>0.0805</td><td>0.8</td></tr><tr><td>Case 3_e</td><td>0.050</td><td>0.1007</td><td>1.0</td></tr><tr><td rowspan="5">Case 4</td><td>Case 4_a</td><td rowspan="5">0.8</td><td rowspan="5">2.58</td><td rowspan="5">1.64</td><td rowspan="5">10.30</td><td>0.0187</td><td>0.0240</td><td rowspan="5">12</td><td rowspan="5">0.085</td><td>0.2</td></tr><tr><td>Case 4_b</td><td>0.0373</td><td>0.0481</td><td>0.4</td></tr><tr><td>Case 4_c</td><td>0.0560</td><td>0.0721</td><td>0.6</td></tr><tr><td>Case 4_d</td><td>0.0747</td><td>0.0962</td><td>0.8</td></tr><tr><td>Case 4_e</td><td>0.0934</td><td>0.1202</td><td>1.0</td></tr><tr><td rowspan="5">Case 5</td><td>Case 5_a</td><td rowspan="5">0.8</td><td rowspan="5">2.58</td><td rowspan="5">1.64</td><td rowspan="5">10.30</td><td>0.0227</td><td>0.0292</td><td rowspan="5">8</td><td rowspan="5">0.103</td><td>0.2</td></tr><tr><td>Case 5_b</td><td>0.0454</td><td>0.0585</td><td>0.4</td></tr><tr><td>Case 5_c</td><td>0.0681</td><td>0.0877</td><td>0.6</td></tr><tr><td>Case 5_d</td><td>0.0908</td><td>0.1169</td><td>0.8</td></tr><tr><td>Case 5_e</td><td>0.1135</td><td>0.1462</td><td>1.0</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>Fig. 1 Initial amplitude spectra for different cases under single-peak spectrum (BFI = 0.6). a Different peak frequencies b Different peak enhancement factors</figcaption></figure>

cases with the same peak enhancement factor (λ) and BFI value (BFI = 0.6), but different peak frequencies. Figure 1b compares the cases with the same peak frequency and BFI value (BFI = 0.6), but various peak enhancement factors.

## 2.4.2 Bimodal spectrum

Although different studies under bimodal spectrum already exist in the literature as mentioned in the introduction, most of them take place in a more complex configuration with different crossing directional sea states at various angles. However, as a special case, the unidirectional co-propagation mixed sea state with different peak frequencies may be the representative case for the wave properties and the nonlinear phenomena in unidirectional waves are obviously stronger. Further, extreme waves can be more easily generated in unidirectional wave trains. It is more reasonable to use unidirectional irregular waves as incident waves in practical engineering design. Therefore, an extensive study in a simplified unidirectional context with different parameters is worthy of analyzing in more details. By assuming that the significant wave height $H_s$ and the mean zero-crossing period $T_z$ of the sea states are fixed with that under single-peak spectrum, the present study is focused on investigating the effects of the type of wave spectra, ranging from single-peak spectrum to bimodal spectrum, for a given sea state parameter. Irregular sea states are often characterized by those integral parameters (for instance during the design of the structures at sea). The distribution of energy content is most of the time not taken into account. This study intends to investigate the influence of this feature on wave statistics, which is essential in the context of design for instance. It is different from and complementary to the purpose of the recent work of Støle-Hentschel et al. (2020), focused on the influence of the presence of the swell on a fixed wind-sea partition.

With respect to the study of the irregular wave trains defined by bimodal spectrum, two dimensionless parameters are commonly used to describe the characteristics of the mixed wave system. One is the Sea-Swell Energy Ratio (SSER) (Rodriguez et al. 2002), which can be defined as the ratio between the zero-order spectrum moment ($m_0$) of the wind-sea frequency band (high-frequency part) to that of the swell wave field (low-frequency part), given as:

$$SSER = \frac{m_{0_{wind-sea}}}{m_{0_{swell}}} \qquad (5)$$

Equivalently, this is the ratio of the significant wave heights squared. The value of SSER describes the energy distribution in the mixed wave system. Wave fields with SSER value smaller than 1.0 represent swell-dominated sea state, and wave conditions with SSER value greater than 1.0 represent wind-sea-dominated sea state. The value of SSER close to 1.0 corresponds to the sea-swell energy equivalent sea state.

Another parameter to describe the frequency separation between the spectral frequency peaks corresponding to the swell and wind-sea system is termed as Intermodal Distance (ID) (Rodriguez et al. 2002), and can be expressed as:

$$ID = \frac{f_{p2} - f_{p1}}{f_{p1} + f_{p2}} \qquad (6)$$

One wave condition with the parameter value of ID being close to 0 represents a sea state with swell and wind-sea peak frequencies in proximity. A wave condition with ID value being close to 0.10 or larger corresponds to a sea state with swell and wind-sea system located in different frequency zones and well separated.

Case 3_c in Table 2, with moderate value of BFI and deep-water depth $k_p h = 16.10$, is taken as a reference for a comparison with the results obtained from a single-peak spectrum. The adopted test cases under bimodal spectrum with the same
water depth ($h = 4.0$ m) and peak enhancement factor ($\lambda = 16$) with that under single-peak spectrum are listed in Table 3, in which different cases with various SSER and ID values are considered. For Cases A–J, the values of SSER are fixed as 1.0, and ID values are changed from 0.02 to 0.35. While for Cases AA–AE, AF–AJ, AK–AO, and AP–AT, the values of ID are fixed as 0.06, 0.08, 0.10, and 0.20, respectively, and SSER values are changed.

With fixed significant wave height of the sea states $H_s$ and the value of SSER, once the significant wave height of one partition $H_{s1}$ (or $H_{s2}$) is given, the other value $H_{s2}$ (or $H_{s1}$) can be determined according to the relation in Eqs. (2) and (5).

The mean zero-crossing period $T_z$, i.e., the characteristic period of the mixed system based on the method of energy weighting, can be expressed as:

$$ \sqrt{\frac{H_{s1}^2 + H_{s2}^2}{H_{s1}^2/T_{z1}^2 + H_{s2}^2/T_{z2}^2}} = T_z \quad (7) $$

$$ T_{z1} = 2\pi\sqrt{\frac{m_{01}}{m_{21}}} \quad T_{z2} = 2\pi\sqrt{\frac{m_{02}}{m_{22}}} \quad (8) $$

where $m_0$ and $m_2$ are the zero- and second-order spectrum moment, respectively.

Similarly, with fixed mean zero-crossing period $T_z$ and the value of ID, once the peak frequency of one partition $f_{p1}$ (or $f_{p2}$) is given, the other value $f_{p2}$ (or $f_{p1}$) can be determined according to the relation in Eqs. (7), (8), and (6). Detailed parameters are listed in Table 3.

Figure 2 illustrates the initial amplitude spectra for different cases under bimodal spectra compared with the single-peak spectrum. In Fig. 2a, SSER is fixed as 1.0 and ID is varying from 0.02 to 0.20. From the amplitude spectra, it can be observed that for smaller ID (narrower peak frequency width), there is no obvious bimodal phenomenon, but only an apparently wider spectral bandwidth. When ID is greater than 0.10 corresponding to Cases F–J, the two peaks of the low- and high-frequency part are completely unaffected by the other one and the spectra appear as the double-peak form, while in Fig. 2b, ID is fixed as 0.20 and SSER is various in each case.

# 3 Numerical setup and data analysis

## 3.1 Numerical model

High Order Spectral (HOS) method was proposed by Dommermuth and Yue (1987) and West et al. (1987). It was initially developed to solve the time evolution of a wave field specified throughout the domain $\eta(x, t=0)$ and $\Phi(x, z=\eta, t=0)$, assuming an infinite homogeneous spatial domain. Ducrozet

<table><caption>**Table 3** Detailed parameters for different cases under bimodal spectrum ($h = 4.0$ m, $\lambda = 16$)</caption><thead><tr><th>Case</th><th>f<sub>p</sub> (Hz)</th><th>H<sub>s</sub> (m)</th><th>ε = k<sub>p</sub>H<sub>s</sub>/2</th><th>SSER</th><th>ID</th></tr></thead><tbody><tr><td>Single</td><td>1.0</td><td>0.03</td><td>0.0604</td><td>—</td><td>—</td></tr><tr><td>Case A</td><td>0.98</td><td>0.0212</td><td>0.0410</td><td>1.00</td><td>0.02</td></tr><tr><td></td><td>1.02</td><td>0.0212</td><td>0.0444</td><td></td><td></td></tr><tr><td>Case B</td><td>0.96</td><td>0.0212</td><td>0.0393</td><td>1.00</td><td>0.04</td></tr><tr><td></td><td>1.04</td><td>0.0212</td><td>0.0461</td><td></td><td></td></tr><tr><td>Case C</td><td>0.94</td><td>0.0212</td><td>0.0376</td><td>1.00</td><td>0.06</td></tr><tr><td></td><td>1.06</td><td>0.0212</td><td>0.0478</td><td></td><td></td></tr><tr><td>Case D</td><td>0.92</td><td>0.0212</td><td>0.0359</td><td>1.00</td><td>0.08</td></tr><tr><td></td><td>1.08</td><td>0.0212</td><td>0.0495</td><td></td><td></td></tr><tr><td>Case E</td><td>0.90</td><td>0.0212</td><td>0.0342</td><td>1.00</td><td>0.10</td></tr><tr><td></td><td>1.09</td><td>0.0212</td><td>0.0512</td><td></td><td></td></tr><tr><td>Case F</td><td>0.84</td><td>0.0212</td><td>0.0302</td><td>1.00</td><td>0.15</td></tr><tr><td></td><td>1.14</td><td>0.0212</td><td>0.0552</td><td></td><td></td></tr><tr><td>Case G</td><td>0.78</td><td>0.0212</td><td>0.0262</td><td>1.00</td><td>0.20</td></tr><tr><td></td><td>1.18</td><td>0.0212</td><td>0.0592</td><td></td><td></td></tr><tr><td>Case H</td><td>0.73</td><td>0.0212</td><td>0.0228</td><td>1.00</td><td>0.25</td></tr><tr><td></td><td>1.21</td><td>0.0212</td><td>0.0627</td><td></td><td></td></tr><tr><td>Case I</td><td>0.67</td><td>0.0212</td><td>0.0192</td><td>1.00</td><td>0.30</td></tr><tr><td></td><td>1.25</td><td>0.0212</td><td>0.0662</td><td></td><td></td></tr><tr><td>Case J</td><td>0.62</td><td>0.0212</td><td>0.0163</td><td>1.00</td><td>0.35</td></tr><tr><td></td><td>1.27</td><td>0.0212</td><td>0.0691</td><td></td><td></td></tr><tr><td>Case AA</td><td>0.98</td><td>0.0276</td><td>0.0534</td><td>0.18</td><td>0.06</td></tr><tr><td></td><td>1.11</td><td>0.0117</td><td>0.0287</td><td></td><td></td></tr><tr><td>Case AB</td><td>0.96</td><td>0.0249</td><td>0.0461</td><td>0.46</td><td>0.06</td></tr><tr><td></td><td>1.08</td><td>0.0168</td><td>0.0396</td><td></td><td></td></tr><tr><td>Case AC</td><td>0.94</td><td>0.0212</td><td>0.0375</td><td>1.00</td><td>0.06</td></tr><tr><td></td><td>1.06</td><td>0.0212</td><td>0.0479</td><td></td><td></td></tr><tr><td>Case AD</td><td>0.92</td><td>0.0173</td><td>0.0294</td><td>2.06</td><td>0.06</td></tr><tr><td></td><td>1.04</td><td>0.0245</td><td>0.0531</td><td></td><td></td></tr><tr><td>Case AE</td><td>0.90</td><td>0.0111</td><td>0.0181</td><td>6.46</td><td>0.06</td></tr><tr><td></td><td>1.01</td><td>0.0279</td><td>0.0578</td><td></td><td></td></tr><tr><td>Case AF</td><td>0.98</td><td>0.0283</td><td>0.0548</td><td>0.12</td><td>0.08</td></tr><tr><td></td><td>1.15</td><td>0.0099</td><td>0.0264</td><td></td><td></td></tr><tr><td>Case AG</td><td>0.95</td><td>0.0254</td><td>0.0461</td><td>0.40</td><td>0.08</td></tr><tr><td></td><td>1.12</td><td>0.0160</td><td>0.0402</td><td></td><td></td></tr><tr><td>Case AH</td><td>0.92</td><td>0.0212</td><td>0.0359</td><td>1.00</td><td>0.08</td></tr><tr><td></td><td>1.08</td><td>0.0212</td><td>0.0495</td><td></td><td></td></tr><tr><td>Case AI</td><td>0.90</td><td>0.0185</td><td>0.0301</td><td>1.68</td><td>0.08</td></tr><tr><td></td><td>1.06</td><td>0.0236</td><td>0.0531</td><td></td><td></td></tr><tr><td>Case AJ</td><td>0.88</td><td>0.0144</td><td>0.0224</td><td>3.44</td><td>0.08</td></tr><tr><td></td><td>1.03</td><td>0.0263</td><td>0.0566</td><td></td><td></td></tr><tr><td>Case AK</td><td>0.98</td><td>0.0287</td><td>0.0555</td><td>0.09</td><td>0.10</td></tr><tr><td></td><td>1.20</td><td>0.0087</td><td>0.0250</td><td></td><td></td></tr><tr><td>Case AL</td><td>0.94</td><td>0.0257</td><td>0.0457</td><td>0.37</td><td>0.10</td></tr><tr><td></td><td>1.15</td><td>0.0155</td><td>0.0412</td><td></td><td></td></tr><tr><td>Case AM</td><td>0.90</td><td>0.0212</td><td>0.0344</td><td>1.00</td><td>0.10</td></tr><tr><td></td><td>1.10</td><td>0.0212</td><td>0.0511</td><td></td><td></td></tr><tr><td>Case AN</td><td>0.88</td><td>0.0192</td><td>0.0300</td><td>1.49</td><td>0.10</td></tr><tr><td></td><td>1.08</td><td>0.0231</td><td>0.0537</td><td></td><td></td></tr><tr><td>Case AO</td><td>0.85</td><td>0.0141</td><td>0.0206</td><td>3.66</td><td>0.10</td></tr><tr><td></td><td>1.04</td><td>0.0265</td><td>0.0575</td><td></td><td></td></tr><tr><td>Case AP</td><td>0.95</td><td>0.0287</td><td>0.0521</td><td>0.09</td><td>0.20</td></tr><tr><td></td><td>1.43</td><td>0.0088</td><td>0.0361</td><td></td><td></td></tr><tr><td>Case AQ</td><td>0.88</td><td>0.0263</td><td>0.0410</td><td>0.31</td><td>0.20</td></tr><tr><td></td><td>1.32</td><td>0.0145</td><td>0.0508</td><td></td><td></td></tr><tr><td>Case AR</td><td>0.79</td><td>0.0212</td><td>0.0264</td><td>1.00</td><td>0.20</td></tr><tr><td></td><td>1.18</td><td>0.0212</td><td>0.0591</td><td></td><td></td></tr><tr><td>Case AS</td><td>0.75</td><td>0.0184</td><td>0.0209</td><td>1.66</td><td>0.20</td></tr><tr><td></td><td>1.13</td><td>0.0237</td><td>0.0603</td><td></td><td></td></tr><tr><td>Case AT</td><td>0.70</td><td>0.0123</td><td>0.0121</td><td>4.99</td><td>0.20</td></tr><tr><td></td><td>1.05</td><td>0.0274</td><td>0.0608</td><td></td><td></td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Initial amplitude spectra for different cases under bimodal spectrum compared with the single-peak spectrum. **a** SSER = 1.0. **b** ID = 0.20</figcaption></figure>

et al. (2012) and Li and Liu (2015) have enhanced this model to represent a water wave tank, including a wave-maker and an absorbing beach. In these models, the initial conditions represent the fluid at rest and the waves are generated with a time-varying boundary condition (representing the wave-maker). The wave tank is of constant water depth *h* and horizontal length *L*<sub>x</sub>. The numerical solution is achieved with the velocity potential being split into the sum of a free surface spectral potential component Φ<sub>f</sub> and a prescribed non-periodic component Φ<sub>w</sub>.

The free surface spectral potential component Φ<sub>f</sub> satisfies the Laplace equation, the free surface boundary conditions, and the bottom condition. This velocity potential is expressed on a set of specific basis functions and described with a given number of modes N<sub>x</sub> (the number of points along the horizontal x direction). It can be solved using the traditional HOS method proposed by Dommermuth and Yue (1987), in which the velocity potential is written in a perturbation series up to an arbitrary order M. The free surface boundary conditions are modified with the non-periodic component acting as forcing terms. This additional velocity potential Φ<sub>w</sub> satisfies the Laplace equation, the wave-maker boundary condition, and the bottom condition. This velocity potential is expressed on another set of specific basis functions with regard to a given number of modes N<sub>z</sub> (the number of points along the vertical dimension, necessary to describe the wave-maker's motion).

The numerical problem can be solved with reference to Bonnefoy et al. (2010). For more details of the present implementation, the reader is referred to Li and Liu (2015).

## 3.2 Numerical setup

In a wave tank configuration, it is reminded that the wave statistics possibly evolve in space and that they are obtained from temporal signals at each location. To ensure that the random waves evolve over a sufficiently long extent in space, the numerical wave tank is at least 50 times of the corresponding peak wavelength for the different peak frequencies in the

numerical simulation. Therefore, the effective length of the numerical wave tank under single-peak spectrum is 220 m, 120 m, and 100 m corresponding to the peak frequency of 0.6 Hz (Case 1), 0.8 Hz (Cases 2, 4, 5), and 1.0 Hz (Case 3) and it is fixed to 220 m in all cases considering bimodal spectrum (Cases A–J and Cases AA–AT), with a water depth of 4.0 m.

Appendix A provides details on the generation of irregular waves. In this paper, the total number of wave components, describing the temporal wave-maker motion (N<sub>f</sub>), is chosen as 200 in all different cases. Irrespective of the case in the following analysis, the frequency range (f<sub>L</sub>, f<sub>H</sub>) is chosen as (0.2 Hz, 3.0 Hz).

With respect to the discretization in space, for cases under single-peak spectrum listed in Table 2, the number of points per peak wavelength N<sub>Lp</sub> is fixed to 30 based on the convergence analysis detailed in the Appendix B section. For bimodal spectra, the frequency content is distributed over a larger frequency range. Then N<sub>Lp</sub> is adopted uniformly as 30, based on the largest peak frequency of the bimodal spectrum for all cases. The discretization in time T<sub>p</sub>/Δt is adopted as 100 and the order of the HOS method M is chosen as 5 in all cases, referring to the convergence analysis in Appendix B.

To obtain sufficiently stable statistics at each fixed location, a large number of waves should be recorded. In order to prevent the reflection from doing too long simulations, the total necessary number of waves expected (for accurate statistics) should be split into different simulations, which have random phases. The effective statistical numerical time is around 6500 s at each location with 10 different random seeds. Then there are at least 5000 waves in total for each configuration. Furthermore, random phases in the same segment for different cases considering various peak frequencies are kept the same value for possible comparison in the time domain.
## 3.3 Data analysis

### 3.3.1 Kurtosis prediction

As an indication of the occurrence probability of the extreme waves, the kurtosis of the free surface elevation is also investigated in this paper. It is defined as:

$$ kurtosis = \lambda_4 = \frac{\langle (\eta - \langle \eta \rangle)^4 \rangle}{\sigma^4} \qquad (9) $$

where $\langle \cdot \rangle$ represents the average over time and $\sigma$ is the standard deviation of $\eta$. In a Gaussian wave process, the value of kurtosis is 3.0, whereas a larger value represents a larger probability of occurrence for the extreme waves.

A theoretical solution for the steady-state kurtosis in finite water depth $\lambda_{4,NB}$, assuming narrow-banded wave field, is introduced for comparison. This kurtosis $\lambda_{4,NB}$ can be written as the sum of two terms accounting for both the dynamic contributions $\lambda_{4,NB}^{dyn}$ and the static contributions $\lambda_{4,NB}^{b}$. The dynamic contributions represent the effect of resonant interactions, which are expected to be significant in the present configuration (Janssen and Bidlot 2009; Fedele 2015) The static contributions correspond to the effect of bound modes, which are the dominant process for the skewness and one of the contributions to the kurtosis. The specific calculation formulas can be found in Janssen and Bidlot (2009) and Ducrozet and Gouin (2017). The expression should represent, for unidirectional, narrow-banded and weakly nonlinear wave field, the asymptotic value of the kurtosis of the free surface elevation. Note that the theoretical spatial evolution can also be defined following Fedele (2015).

### 3.3.2 Probability density functions: Naess, MER, and GC

The statistical analysis of random wave height is an important parameter in the study of extreme wave characteristics. Various theoretical and empirical models have been proposed to characterize the wave height probability distribution.

In a stationary, Gaussian and extremely narrow-banded process, the wave height is regarded as twice the envelope amplitude and that these are distributed according a Rayleigh probability distribution (Longuet-Higgins 1952). It is given by:

$$ E(H) = \mathrm{EXP}\left(-\frac{H^2}{8}\right) \qquad (10) $$

in which $H$ is the non-dimensionalized wave height with the square root of the zero-order spectral momentum (i.e., $\sqrt{m_0}$). However, the assumption of the wave height as being twice the amplitude of the envelop amplitude is not totally exact due

to the modulated structure of a narrow-banded process and the time lag between a crest and the adjacent trough. Naess (1985) derived a linear distribution model for the crest-to-trough wave height in a narrow-banded wave trains given by:

$$ E(H) = \mathrm{EXP}\left[-\frac{H^2}{4(1-\rho(\tau/2))}\right] \qquad (11) $$

But Naess (1985) did not define the key parameter $\rho(\tau/2)$ in his distribution exactly. Its asymptotic approximation was given by Boccotti (2000), that $\rho(\tau/2)$ represents the value of the normalized autocorrelation function of the sea surface elevation at the time when it attains its first minimum.

In nature, most of the wave fields are nonlinear. For stronger nonlinear waves, using Rayleigh or Naess distribution will sometimes produce great deviation. Hence, the Rayleigh or Naess distribution is no longer applicable. Considering the effect of nonlinear interaction, Mori and Janssen (2006) developed a distribution called a modified Edgeworth-Rayleigh distribution (MER) under the assumption of weak nonlinearity, narrow spectrum, and wave height twice of wave amplitude. It is expressed as:

$$ E(H) = \mathrm{EXP}\left(-\frac{H^2}{8}\right) \left[1 + \left(\frac{\pi}{\sqrt{3}}\mathrm{BFI}^2\right) \frac{H^2}{384}(H^2-16)\right] \qquad (12) $$

where BFI is chosen as the deep-water BFI in this study.

And then, scholars found that the third-order nonlinearity of the wave surface elevation has a great influence on the statistical characteristics of random wave trains and meanwhile the fourth-order moment of wave surface elevation is directly related to the third-order nonlinearity from the perspective of statistical analysis. Taking kurtosis up to the third order into account, Tayfun and Fedele (2007) and Fedele et al. (2016) have proposed a Gram-Charlier (GC) distribution model, given as:

$$ E(H) = \mathrm{EXP}\left(-\frac{H^2}{8}\right) \left[1 + \frac{\Lambda}{1024} H^2 (H^2 - 16)\right] \qquad (13) $$

$$ \Lambda = \lambda_{40} + 2\lambda_{22} + \lambda_{04} \qquad (14) $$

$$ \lambda_{mn} = \langle \eta^m \hat{\eta}^n \rangle / \sigma^{m+n} + (-1)^{m/2} (m-1)(n-1); \quad m+n=4 \qquad (15) $$

where $\hat{\eta}$ is the Hilbert transform of the free surface elevation.

Note that the difference between MER and GC model is that in MER model the influence of third-order nonlinearity on wave surface elevation is obtained by theoretical analysis and can be directly calculated by the given BFI values of the wave trains. However, in GC model, although Eq. (13) is also derived by theoretical consideration, the nonlinear effect is based on the parameter $\Lambda$, which should be statistically determined from the analyzed wave series according to Eqs. (14) and (15). So, the latter slightly changes with the wave propagation and can only be used as a reference model for the wave height distribution analysis at some specific position for bimodal waves.
MER and GC nonlinear model will be used as a reference
distribution for the following wave height distribution analysis.

**4 Numerical results and discussions**

This section is dedicated to the analysis of the different sea
states presented in Section 2.4. Unidirectional single-peak
spectra are investigated in detail, summarizing the main ef-
fects of BFI and relative water depth. Then, the cases defined
by bimodal spectra, assuming the wave field is described with
a given set of integral parameters (H<sub>s</sub> and T<sub>z</sub>) but different
SSER and ID, are studied to investigate the effects of the
distribution of the energy content in bimodal spectra on sta-
tistical properties.

4.1 Single-peak spectrum

Figure 3 presents the evolution of the kurtosis along the wave
flume for the different cases with various BFI values in the
numerical investigations, in which the horizontal axis repre-
sents the distance away from the wave-maker normalized by
the corresponding peak wavelength. For the case with smaller
BFI (BFI = 0.2), the kurtosis is almost constant along the wave
flume. However, for the cases with larger BFI, the values of
kurtosis obviously grow along the wave flume and deviate
from the Gaussian value of 3.0. The cases with larger BFI
have faster growth rates and larger maximum values. When
the value of BFI is 1.0, the maximum value of kurtosis can
reach up to 4.75. These phenomena have been observed in
different physical experiments (Onorato et al. 2004; Shemer
and Sergeeva 2009), theoretical analyses (Janssen 2003;
Janssen and Bidlot 2009; Fedele 2015), and fully nonlinear
numerical simulations (Onorato et al. 2006; Regev et al. 2008;
Gramstad and Trulsen 2010). They had proposed that these
phenomena are caused by modulation instability that depends
on the wave steepness and spectral bandwidth (BFI).

Although the evolution of kurtosis is strongly dependent on
BFI value of the sea states, it can be observed from Fig. 3 that
the maxima of the kurtosis as well as its asymptotic value are
different from each other for the cases with the same BFI
value. Consequently, the evolutions of kurtosis along the
wave flume for different relative water depths and peak en-
hancement factors with the same BFI values are compared in
Fig. 4, in which a theoretical solution for the steady-state
kurtosis in finite water depth λ<sub>4</sub>, NB is introduced as dash-
dotted lines for comparison. From Fig. 4, it can be clearly seen
that these theoretical values are significantly smaller than the
numerical results from the fully nonlinear numerical model,
especially for the cases with larger BFI value and deeper water
depth. It indicates that the physics of the process requires to be
taken into account higher dispersive and nonlinear effects than
this simple theoretical model.

More detailed analysis is presented in Fig. 4, where the sea
states with the same BFI are studied, but with varying other
parameters. Figure 4a, c, and e investigate the influence of the
relative water depth *k*<sub>p</sub>*h* on the spatial evolution of the kurto-
sis, while Fig. 4b, d, and f study the effect of the wave steep-
ness ε (or eq. Δ*f*/f<sub>p</sub>, since BFI is fixed) on the same spatial
evolution. It can be noted that, for Fig. 4b, d, and f with
constant BFI value and the same water depth, the horizontal
axis is normalized by *L*<sub>p</sub>/*ε*<sup>2</sup>.

From the observation, it can be found that, with constant BFI
value, the influence of relative water depth is significant, even
for values of k_p h > 5.8, often considered deep-water regime,
while the influence of the wave steepness (or Δf/f_p) appears
limited to the growth rate of the kurtosis, the asymptotic values
being nearly identical. It is found that for a fixed BFI and k_p h, a
steeper and narrower sea state has a faster growth rate. This is
associated with the higher-order dispersion and nonlinear ef-
fects than the theoretical prediction available (Fedele 2015).
The latter is indeed based on narrow-banded third-order solu-
tion of the problem that does not predict such behavior, the
evolution being only a function of BFI and k_p h.

Considering that the modified Benjamin Feir Index B<sub>S</sub> in-
cludes the influence of finite water depth, to investigate the
influence factors on the variation of the kurtosis under single-
peak spectrum, it is more relevant to take into account these two
independent factors: the deep-water BFI and the relative water
depth k<sub>p</sub>h. Therefore, a series of simulations with a fixed peak
frequency f<sub>p</sub> = 0.8 Hz for deep-water BFI = 0.6, 0.8, and 1.0 are
carried out. Several water depths are tested with h = 4 m, 3 m,
2 m, 1.5 m, 1.2 m, 1 m, 0.8 m, and 0.6 m, respectively. It should
be noted that as the water depth becomes shallower, the wave-
number k<sub>p</sub> becomes larger due to the dispersion equation in-
volved the influence of water depth. Then the relative water
depth corresponds to k<sub>p</sub>h = 10.32, 7.74, 5.16, 3.87, 3.10, 2.60,
2.12, and 1.66, respectively. At the same time, in order to keep
the deep-water BFI fixed at 0.6, 0.8, and 1, the significant wave
height needs to be adjusted accordingly. Subsequently, the sta-
tistics of the maximum kurtosis along the numerical wave
flume in each condition are analyzed. The contour plot of the
maximum kurtosis with respect to different initial deep-water
BFI and 1/(k<sub>p</sub>h) is illustrated in Fig. 5, where cubic interpolation
is used to plot the contour. Figure 5 a gives the numerical result
simulated by HOS method and Fig. 5b gives the narrow-banded
weakly nonlinear theoretical values obtained from Janssen and
Bidlot (2009) and Ducrozet and Gouin (2017). It should be
pointed out that the effect of finite water depth has been con-
sidered in the theoretical predictions above (i.e., in the calcula-
tion of the theoretical value, the modified Benjamin Feir index
B<sub>S</sub> replaces the deep-water BFI). In Fig. 5, the dotted black line
level represents the classical limit of the deep-water conditions
corresponding to k<sub>p</sub>h = π. We note that large variations of the
maximum kurtosis are observed in this deep-water regime
(k<sub>p</sub>h > π, below the dotted line).
<figure><img src="image_6.png" /><figcaption>(a) Case 1 ($k_p h = 5.80, \lambda = 16$)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Case 2 ($k_p h = 10.30, \lambda = 16$)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) Case 3 ($k_p h = 16.10, \lambda = 16$)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) Case 4 ($k_p h = 10.30, \lambda = 12$)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(e) Case 5 ($k_p h = 10.30, \lambda = 8$)</figcaption></figure>

<caption>**Fig. 3** The comparisons of the evolution of kurtosis for various BFI values under single-peak spectrum. **a** Case 1 ($k_p h = 5.80, \lambda = 16$). **b** Case 2 ($k_p h = 10.30, \lambda = 16$). **c** Case 3 ($k_p h = 16.10, \lambda = 16$). **d** Case 4 ($k_p h = 10.30, \lambda = 12$). **e** Case 5 ($k_p h = 10.30, \lambda = 8$)</caption>

Even though the value of deep-water BFI is fixed, the maximum values of the kurtosis are different owing to the different relative water depths. It is observed that the deeper the water depth is, the smaller the influence on the maximum kurtosis is. However, it is worth mentioning that at the classical deep-water limit ($k_p h = \pi$), the deviation from Gaussian sea state (a.k.a. kurtosis enhancement, $\lambda_4=3$) is approximately 25% smaller than the deep-water solution (that is reached at a far larger relative water depth $k_p h$). The kurtosis becomes smaller with a shallower water depth. This implies that the shallower water depth can depress the wave instability and larger wave development. It is due to the fact that the reduced water depth decreases the effects of modulation instability (Janssen and Onorato 2007). Hence, the values of kurtosis are different from each other in the so-called deep-water conditions, as

defined with the peak frequency and the same value of deep-water BFI.

The effect of the water depth on the wave dynamics has been studied extensively in the literature. Modulational instability vanishes as the depth decreases below the critical value corresponding to $k_p h \approx 1.363$ and this can be interpreted by a stability analysis of the Nonlinear Schrödinger (NLS) equation (Slunyaev 2005). Analyzing the coefficients of this equation tells us that the modulational instability is nonexistent below this threshold. However, for larger $k_p h$ values, and for instance the ones we study here, the coefficient is not constant. As a result, the modulational instability is still influenced by the water depth, so are the wave statistics.

Comparing the numerical results with the theoretical values in Fig. 5, these two have similar trends, but the theoretical
<figure><img src="image_11.png" /><figcaption>**Fig. 4** The evolution of kurtosis along the wave flume for different relative water depths and peak enhancement factors under single-peak spectrum. **a** BFI = 0.6, $\lambda$ = 16. **b** BFI = 0.6, $k_p h$ = 10.30. **c** BFI = 0.8, $\lambda$ = 16. **d** BFI = 0.8, $k_p h$ = 10.30. **e** BFI = 1.0, $\lambda$ = 16. **f** BFI = 1.0, $k_p h$ = 10.30</figcaption></figure>

values underestimate the actual fully nonlinear numerical simulation results. The maximum kurtosis can reach more than 4.0 for larger deep-water BFI, much larger than the theoretical values. This implies that it is quite necessary to adopt a fully nonlinear numerical model to simulate sea state waves with high-order nonlinearity and to predict accurately the occurrence of large waves, especially in deep water.

Generally, the kurtosis is associated with the presence of large waves. Complementary to the evolution of the kurtosis, the probability density function should be further studied. Figure 6 displays the exceedance probabilities of wave height (counted by zero up-crossing method) at the location of the maximum value of kurtosis. In the figure, MER distribution, considering the nonlinear interaction, is plotted for comparison. It can be easily observed that for smaller BFI, the wave height distribution almost obeys the MER distribution, and for

larger BFI there is a little deviation from the MER distribution. MER distribution can better predict the probability of large waves ($1.5 < H/H_s < 2.0$), but underestimates the statistics of the wave height distribution when $H/H_s$ is larger than 2.0, which is the classical value to identify events are rogue or freak waves (Klinting and Sand 1987). Meanwhile, the exceedance probabilities of wave height at maximum values of kurtosis considering the influence of the relative water depth and peak enhancement factor are exhibited in Fig. 7. The evolution of the wave height distribution is consistent with that obtained from the kurtosis. From Fig. 7a, it can be observed that for deeper water depth, MER model can better predict the probability of large waves ($1.5 < H/H_s < 2.0$), but also underestimates the statistics of the wave height distribution when $H/H_s$ is larger than 2.0. From Fig. 7b, irrespective of the peak enhancement factor, their wave height distributions
<figure><img src="image_12.png" /><figcaption>**Fig. 5** Contour plot of the maximum kurtosis with regard to different initial deep-water BFI and $1/(k_p h)$ (the dotted black line level represents the limit of the deep-water condition). **a** Numerical results. **b** Theoretical values</figcaption></figure>

keep basically consistent. And when $1.5 < H/H_s < 2.0$, the wave height distribution is very close to the MER theoretical prediction, which still underpredicts the occurrence of extreme waves ($H/H_s > 2.0$).

## 4.2 Bimodal spectrum

In sing-peak spectrum, modulational instability plays an important role during extreme wave generation, and BFI value

<figure><img src="image_13.png" /><figcaption>**Fig. 6** The influence of BFI on exceedance probability of wave height under single-peak spectrum for Case 2 ($k_p h = 10.30$, $\lambda = 16$) (the solid lines represent MER distribution)</figcaption></figure>

can be an indicator, while, for bimodal spectrum, the energy is divided into two partitions. The nonlinear interaction process is more complex. It is not easy to define an index to estimate the extreme wave generation. So, in this section, two spectral shape parameters are considered. We fix the integral parameters (the significant wave height $H_s$ and the mean zero-crossing period $T_z$ of the sea states) and vary the set of parameters: the energy distribution SSER and the frequency separation ID. As detailed in the introduction, one of the present objectives is to assess the importance of the wave energy distribution in the bimodal sea states and not only through the abovementioned integral parameters.

Figure 8 shows the evolution of the kurtosis along the wave flume for different cases under bimodal spectrum, where the solid black lines represent the results for the single-peak spectrum for a reference. The horizontal axis represents the distance away from the wave-maker normalized by the wavelength $L_z$ corresponding to the mean zero-crossing period (obtained from linear dispersion relation). It should be noted that comparing to Fig. 3c, the evolution takes place on a longer spatial extent in order to ensure sufficient evolution for bimodal waves. In Fig. 8a, the results of Cases A–J with fixed SSER but various ID values have been compared. The variation under single-peak spectrum has the faster growth rate and the larger maximum kurtosis than that under the bimodal spectrum. As ID becomes larger, the maximum kurtosis becomes gradually smaller. For Cases F–J where there is obvious bimodal structure of the spectrum with ID larger than 0.10, the values of kurtosis are close to Gaussian value of 3.0 and the changes during the wave propagation are not significant. This is consistent with the findings of Stole-Hentschel et al. (2020). Their study indeed reveals that considering a bimodal spectrum, the non-Gaussian behavior of one partition of the spectrum may be hidden due to a simple superposition. However, additional phenomena can be
<figure><img src="image_14.png" /><figcaption>**Fig. 7** The influence of relative water depth and peak enhancement factor on exceedance probability of wave height under single-peak spectrum (the black solid lines represent MER distribution). **a** BFI = 0.8, $\lambda$ = 16. **b** BFI = 0.8, $k_p h$ = 10.30</figcaption></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 8** The variation of the kurtosis along the wave flume for different cases under bimodal spectra compared with that under the single-peak spectrum. **a** SSER = 1.0. **b** ID = 0.06. **c** ID = 0.20</figcaption></figure>
found in Fig. 8b and c in which the energy repartition between wind-sea and swell (SSER) is varied. In each subplot, as the value of SSER becomes larger, the growth rate and the maximum value of the kurtosis become smaller for values lower than 1.0. When SSER reaches 1.0, which represents the sea-swell energy equivalent sea state, the kurtosis gets the minimum value (represented by the blue lines). Subsequently, as SSER continues increasing, the kurtosis becomes gradually larger. The amount of energy in each partition has a significant effect on the combined sea state statistics.

Figure 9 displays the relationship between the maximum values of kurtosis and ID values under sea-swell energy equivalent sea state that helps to understand the effect of the ID value on the kurtosis. We find that when the ID value is smaller than 0.10, the maximum value of kurtosis becomes smaller with increasing ID. This can be explained that when the two peaks of the input spectra have not been completely separated, the effect of increasing ID is mainly to have broaden the spectrum, which leads to the decrease of the effective BFI at constant $H_s$. The maximum values of the kurtosis decrease consequently. However, when the ID values are larger than 0.10, the maximum values of the kurtosis are nearly constant around 3.3 and there is almost no change in statistics of the bimodal sea. This indicates that statistics of the combined (bimodal) sea states with two completely separated peaks in the input spectra deviate moderately from Gaussian state. Note that this value of 3.3 is the maximum of a space-evolving kurtosis. It is driven by nonlinear interactions (also known as dynamic kurtosis), the importance of bound modes being limited (just depend on steepness).

Similarly, the variation of the maximum values of the kurtosis with various SSER values but fixed ID value is demonstrated in Fig. 10. The black dotted line represents the maximum kurtosis of the wave trains under single-peak spectrum. For swell-dominated sea (i.e., SSER < 1.0), the maximum value of kurtosis becomes smaller with SSER increasing up to 1.0. Thereafter, as SSER becomes larger than 1.0 (i.e., wind-sea-dominated), the maximum value of kurtosis becomes larger with increasing SSER. Meanwhile, the line for

<figure><img src="image_17.png" /><figcaption>**Fig. 9** The relationship between the maximum values of kurtosis and ID values under sea-swell energy equivalent sea state (SSER = 1.0)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 10** The variation of the maximum values of the kurtosis with various SSER values</figcaption></figure>

smaller ID value always lies above the larger one. These phenomena indicate that under the swell-dominated conditions or the wind-sea-dominated conditions, all the kurtosis maxima are larger than those of the case under sea-swell energy equivalent sea state. It demonstrates that with the same value of ID, the wave trains under sea-swell energy equivalent sea state are the less prominent to extreme wave appearance.

According to the study of Støle-Hentschel et al. (2020), it is necessary to partition the mixed sea states and analysis should be performed on each partition to better interpret the combined sea state wave statistics. Cases D–G from Table 3 are chosen to investigate the statistics difference between the mixed system and the two independent systems. The BFI value of each partition (swell and wind-sea) is reported in Table 4. As a reminder, the reference single-peak spectrum with identical $H_s$ and $T_z$ exhibits BFI = 0.6. The swell is weakly nonlinear relative to the wind-sea. Figure 11 gives the initial amplitude spectra of selected cases for further comparison with single independent system conditions. The vertical dotted line marks the splitting frequency $f_m$ that is employed to separate the combined swell and wind-sea (mixed) into swell and wind-sea contributions.

From the definition of kurtosis (Eq. (9)), Støle-Hentschel et al. (2020) calculated the kurtosis of combined sea states

<table><caption>**Table 4** Detailed parameters of various cases to investigate the statistics difference between the mixed system and the independent system</caption><thead><tr><th>Case</th><th>f<sub>p</sub> (Hz)</th><th>H<sub>s</sub> (m)</th><th>ε = k<sub>p</sub>H<sub>s</sub>/2</th><th>BFI</th></tr></thead><tbody><tr><td rowspan="2">Case D</td><td>0.92</td><td>0.0212</td><td>0.0359</td><td>0.36</td></tr><tr><td>1.08</td><td>0.0212</td><td>0.0495</td><td>0.49</td></tr><tr><td rowspan="2">Case E</td><td>0.90</td><td>0.0212</td><td>0.0342</td><td>0.34</td></tr><tr><td>1.09</td><td>0.0212</td><td>0.0512</td><td>0.51</td></tr><tr><td rowspan="2">Case F</td><td>0.84</td><td>0.0212</td><td>0.0302</td><td>0.30</td></tr><tr><td>1.14</td><td>0.0212</td><td>0.0552</td><td>0.55</td></tr><tr><td rowspan="2">Case G</td><td>0.78</td><td>0.0212</td><td>0.0262</td><td>0.26</td></tr><tr><td>1.18</td><td>0.0212</td><td>0.0592</td><td>0.59</td></tr></tbody></table>
<figure><img src="image_19.png" /><figcaption>Fig. 11 Initial amplitude spectra of selected cases for further comparison with single independent system conditions</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 12 Comparison of the statistic of the mixed sea state with the results under single independent system conditions and the theoretical results obtained from Eq. (16). (The black solid lines represent the theoretical</figcaption></figure>

based on the assumption that the two systems are independent (i.e., the free surface is a simple superposition of the wind-sea and swell):

$$ \lambda_{4,mixed} = \frac{\lambda_{4,swell} \sigma_{swell}^4 + 6 \sigma_{swell}^2 \sigma_{wind-sea}^2 + \lambda_{4,wind-sea} \sigma_{wind-sea}^4}{(\sigma_{swell}^2 + \sigma_{wind-sea}^2)^2} \quad (16) $$

where $\sigma^2$ is the variance of the corresponding free surface elevation.

Figure 12 compares the kurtosis evolution of the mixed sea state with the results obtained under single independent system conditions. The theoretical results obtained from Eq. (16) are also included. Irrespective of the ID value, analysis of the entire wave system shows that the mixed sea has milder extreme wave statistics than the pure wind-sea. When one part of the sea state (wind-sea) exhibits strongly nonlinear behavior, put together with a swell with a mildly nonlinear feature, it looks like weakly nonlinear according to the superposition. However, the physics taking place in the mixed sea state are

<figure><img src="image_21.png" /><figcaption>results, and red solid lines represent the result of the mixed sea states.) a Case D ID = 0.08. b Case E ID = 0.10. c Case F ID = 0.15. d Case G ID = 0.20</figcaption></figure>
complex, especially in the wind-sea partition which exhibits
strong non-Gaussian properties. When the ID values are larger
than 0.10 (Fig. 12c and d), i.e., the swell partition and the
wind-sea partition are separated, the statistic agrees with the
theoretical value well. This means that there are no significant
interactions between the wind-sea and the swell in this com-
bined sea state. While as the ID values are less than 0.10
(Fig. 12a and b), the deviations between the numerical mixed
system and the theoretical values become larger. This demon-
strates that for the cases with small ID values, the mixed wave
system exhibits nonlinear interactions between partitions. It is
associated with the fact that the theoretical result does not
involve the influence of one partition (swell or wind-sea) on
the other one. It is a simple averaging that do not provide the
whole information about the physical processes at play. In this
configuration (ID < 0.10), it is consequently not advisable to
consider the combined sea state as a simple linear superposi-
tion of two independent partitions.

As a general point, it is of major importance to stress that
when facing multimodal systems, the complexity of the underly-
ing physics (typically resulting from nonlinear interactions) can-
not be deduced by an analysis limited to the statistical content of
the whole sea state (such as kurtosis or probability distribution)
without information about the energy distribution (i.e., the spec-
trum shape). It is indeed demonstrated here that wave fields
characterized by weakly non-Gaussian features may actually
contain strongly nonlinear phenomena (in the present case, the
wind-sea partition exhibits a strongly non-Gaussian behavior),
and nonlinear interactions between partitions of mixed sea also
play significant roles as ID < 0.10. This is essential in a numerical
point of view, assessing the importance of highly nonlinear
models, such as the HOS method at use in the present study,
compared to weakly nonlinear approaches.

Finally, exceedance probability distributions of the wave
height at the location with maximum kurtosis are given in
Fig. 13 for different choices of parameters. The MER distri-
bution is represented as the black solid line, and the case under
single-peak spectrum is represented as black scatters for ref-
erences in the figures. It can be observed that the wave height
distributions have certain degree of departure from the MER
distribution when H/Hs is larger than 1.5 in all subplots. In
Fig. 13a, for smaller ID values (Cases A and C), the wave
height distributions are close to MER distribution and for
larger ID values (Cases E-I) the departures become gradually
larger as ID values become larger. However, until the swell
partition and the wind-sea partition are completely separated
(Cases G and I), the departure remains unchanged due to the
global weak nonlinearity in the wave trains. In Fig. 13b, it can
be reported that MER distribution lies in the middle of the
analyzed results and the departure from the MER line is not
mono changing. The larger departure can be observed for the
single case and case with SSER = 1.0. All these results from
the exceedance probability distribution of wave height are
consistent with those obtained from the study of kurtosis as
shown in Figs. 9 and 10. It implies that the wave trains under
bimodal spectrum with equivalent energy distribution and
larger intermodal distance between the two peaks of the initial
spectrum are the less prominent to extreme wave appearance.

Figure 14 compares the wave height distribution at the
location with maximum kurtosis between numerical results
and the corresponding results of MER and GC model to in-
vestigate the applicability of the theoretical models under dif-
ferent spectral types. The results are derived for the case under
single-peak spectrum and typical cases under bimodal spec-
trum. It can be clearly seen that comparing to MER model, GC
model has better approximation for extreme wave statistics,

<figure><img src="image_22.png" /><figcaption>**Fig. 13** Exceedance probabilities of the wave height distribution at the location with maximum kurtosis for different cases under bimodal spectrum compared with that under the single-peak spectrum (the black solid lines represent MER distribution). **a** SSER = 1.0. **b** ID = 0.10</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Fig. 14** Exceedance probabilities of the wave height distribution at maximum kurtosis for some different cases compared with the corresponding MER and GC distribution</figcaption></figure>

which is resulted from the fact that during the wave propaga-
tion, the third-order nonlinearity is larger than the theoretical
value due to the wave-wave interaction. Meanwhile, there is
still some deviation from the theoretical prediction for the
rogue waves ($H/H_s > 2$), which can be explained by higher-
order effects (higher than third order).

**5 Conclusions**

In this paper, the statistical properties of unidirectional wave
trains under single-peak and bimodal spectra are analyzed
from a series of numerical simulations taking into accounts
the nonlinearities in the process. The analysis is limited to
long-crested waves, as a typical configuration during the de-
sign of structures at sea to define the environmental conditions
of operation.

For the wave trains under single-peak spectrum, as is
known, the kurtosis and large wave generation probability

are dependent on the BFI value significantly. The sea states
with larger BFI generally have larger kurtosis and are more
prominent to generate extreme waves. However, the relative
water depth has significant influence on the large wave occur-
rence probability. It is demonstrated that in conditions often
referred as deep water ($k_p h > \pi$), the influence of water depth
on the wave statistics can be significant. A detailed analysis of
the combined effect of deep-water BFI and relative water
depth is provided in the present manuscript. It is shown that
the maximum excess kurtosis associated $k_p h = \pi$ can be up to
25% lower than the one of exact deep-water configuration.

Furthermore, the cases of sea states defined with fixed sig-
nificant wave height ($H_s$) and mean zero-crossing period ($T_z$)
but different energy distributions (SSER) or peak frequency
widths (ID) are investigated. It covers single-peak and bimod-
al spectra and intends to assess the major importance of these
features on the physics at play. An extensive study is provid-
ed, with the influence on statistical properties of two parame-
ters SSER and ID considered. The limitations associated to the
description of a sea state with the only two integral parameters $H_s$ and $T_z$ are demonstrated. Compared with the results under single-peak spectrum, the kurtosis and wave height distribution under single-peak spectrum are obviously larger than those under bimodal spectrum, irrespective of the case. This indicates that the sea states described by bimodal spectra are less likely to contain extreme (freak) waves. From an ocean engineering design point of view, the classical unimodal configuration is consequently conservative. However, taking into account the fact that 15 to 25% of the sea states are mixed ones, a more accurate estimation of extreme events and responses can be achieved with the bimodal description.

With respect to the cases with fixed SSER value, when the ID value is smaller than 0.10, the maximum value of kurtosis becomes smaller with increasing ID. As the swell and the wind-sea partitions are completely separated (ID > 0.10), the maximum values of the kurtosis of the combined sea states are constant around 3.3 and the wave height distributions remain a certain degree of departure from the MER distribution. While with respect to the cases with ID fixed with 0.10, when SSER value is close to 1.0, the kurtosis gets the minimum value, and the wave height distribution has a larger departure from the MER distribution. It can be concluded that with fixed ID, the wave train under sea-swell energy equivalent state is the less prominent to extreme wave appearance.

Additionally, comparing to MER model, GC distribution can better predict the probability of large waves (1.5 < $H/H_s$ < 2.0), especially for unidirectional bimodal wave trains, but underestimates the statistics of the wave height distribution when $H/H_s$ is larger than 2.0 (caused by higher than third-order effects) for long-crested waves in both single-peak state and bimodal state.

Furthermore, we have compared the mixed sea states with the corresponding single independent systems. When one part of the sea state (wind-sea) exhibits strongly nonlinear behavior, and is mixed together with a swell with a mildly nonlinear feature, it behaves in a weakly nonlinear way, according to the superposition. Then, it has been demonstrated that the complexity of the underlying physics of a given sea state (mainly the effect of nonlinear interactions, resonant or not) cannot be deduced by an analysis limited to the statistical content of the combined sea state. The way energy is distributed among frequencies plays a major role. This is of huge importance in the context of wave modeling, pointing out the necessity to use highly nonlinear models, even if the statistics may appear weakly non-Gaussian.

**Funding** This research was supported by the National Natural Science Foundation of China (51739010, 51879037), the National Key Research and Development Plan (2016YFC1401405, 2016YFE0200100), the Fundamental Research Funds for the Central Universities of China, and the National Scholarship for Building High Level Universities, China Scholarship Council (201706060082).

# Appendix A: Generation of irregular waves

The wave-maker motion is determined by the target free surface elevation and the linear transfer function, which depends on its geometry. According to the linear superposition theory, the two-dimensional irregular free surface elevation can be represented as a superposition of regular wave components with different frequencies:

$$ \eta(x, t) = \sum_{i=1}^{N_f} a_i \cos(k_i x - \omega_i t + \varphi_i) \quad (17) $$

where the subscript i stands for the ith wave component; $N_f$ is the total number of wave components; $a_i$, $k_i$, and $\omega_i$ are the wave amplitude, wave number, and wave frequency of each component; and $\varphi_i$ is the initial random phase, which is uniformly distributed from 0 to $2\pi$. The value of $k_i$ can be obtained through the dispersion equation:

$$ \omega_i^2 = gk_i \thicksim k_i h \quad (18) $$

with $\thicksim$ the gravity acceleration. The wave amplitude of each component $a_i$ is defined from the wave spectrum:

$$ a_i = \sqrt{2S(\omega_i)\Delta\omega_i} \quad (19) $$

where $\Delta\omega_i$ is the division widths of the frequencies, $\Delta\omega_i = 2\pi(f_H - f_L)/(N_f - 1)$. ($f_L$, $f_H$) is the considered frequency range. To avoid the periodic time repetition, the actual frequency of the ith wave component $\omega_i$ is selected randomly in the frequency spin as

$$ \omega_i = (i - 1 + \text{rand})\Delta\omega_i \quad (20) $$

where rand is a random number distributed from 0 to 1.

# Appendix B: Convergence analysis

The data from the irregular wave experiment under single-peak spectrum by Li et al. (2013) are used to conduct the convergence analysis of the established numerical model. The experiment was carried out in the State Key Laboratory of Coastal and Offshore Engineering in Dalian University of Technology, China. The wave flume is 69.0 m long, 2.0 m wide, and 1.8 m deep. The water depth for the experiment is set at 1.2 m. Eleven resistive wave gauges are arranged along the length of the wave flume to record the free surface elevations. The experimental setup is displayed in Fig. 15.

Irregular waves are characterized by the JONSWAP spectrum with random initial phases. Details of the experimental parameters are listed in Table 5. Two cases with different BFI values of 0.51 and 0.87, which represent different stabilities of the wave group, have been chosen by varying the peak period $T_p$ and the peak enhancement factor $\gamma$ in the JONSWAP
<figure><img src="image_24.png" /><figcaption>**Fig. 15** The setup of the experiment</figcaption></figure>

spectrum. BFI is defined as deep-water BFI = $\frac{\sqrt{2\epsilon}}{2\Delta f/f_p}$, with $\epsilon = k_p H_s/2$ the wave steepness and $\Delta f/f_p$ the relative spectral bandwidth. For the two cases, the significant wave height $H_s$ is the same. However, case II has a deeper relative water depth $k_p h$, larger wave steepness $\epsilon$, and larger BFI value, which represents a stronger nonlinear wave condition than case I. Thus, the convergence analysis is performed with the irregular wave trains for case II.

The numerical wave tank replicas the experimental setup. Its length is enlarged to 80 m to ensure correct wave absorption. The error is measured thanks to the free surface elevation recorded at the middle of the wave tank on a time window fixed to 100$T_p$ and calculated as:

$$ \text{error} = \frac{\int_t |\eta_{\text{numerical}} - \eta_{\text{experiment}}|}{\int_t |\eta_{\text{experiment}}|} \quad (21) $$

The convergence analysis is carried out with the parameters $N_c = N_x/4$, $M = 5$, and the result of the convergence analysis with respect to the number of points per peak wavelength $N_{Lp}$ is represented in Fig. 16. The figure displays two lines representing the convergence rates of 2 and 3 as a reference. The numerical convergence rate is slightly larger than 2nd order, which is the theoretical expected value (Ducrozet et al. 2012). Considering the numerical simulation of the uni-directional irregular wave trains, the overall error is about 5% with 30 points per corresponding peak wavelength. This value of $N_{Lp} = 30$ is chosen as converged parameter for the rest of the study. It ensures the accuracy of the numerical simulation as well as a fast solution. Similar convergence analysis has been conducted for the discretization in time (time-step $\Delta t$) as well as the order of nonlinearity of the method $M$. The final numerical configuration is $N_{Lp} = 30$, $T_p/\Delta t = 100$, and $M = 5$.

The comparison of the input target spectrum and numerically generated one at $x = 5$ m is exhibited in Fig. 17.

<table><caption>**Table 5** Parameters of the physical experiment</caption><thead><tr><th>Case</th><th><i>T</i><sub>p</sub> (s)</th><th><i>k</i><sub>p</sub>h</th><th>γ</th><th><i>H</i><sub>s</sub> (m)</th><th>ε = <i>k</i><sub>p</sub><i>H</i><sub>s</sub>/2</th><th>BFI</th></tr></thead><tbody><tr><td>Case I</td><td>1.5</td><td>2.20</td><td>7.0</td><td>0.055</td><td>0.05</td><td>0.51</td></tr><tr><td>Case II</td><td>1.0</td><td>4.82</td><td>3.3</td><td>0.055</td><td>0.11</td><td>0.87</td></tr></tbody></table>

The comparison ensures the correct wave generation procedure. In addition, the results obtained with a linear solver ($M = 1$) is added as a reference. It can be observed that the generated spectrum with $M = 1$ is exactly the target spectrum, which demonstrates the effectiveness of the established numerical model in dealing with random wave trains. Consequently, the discrepancies between the input target spectrum and the generated one in a nonlinear context ($M = 5$) are only due to the nonlinear effects in the process of the wave propagation, the generation of waves being strictly equivalent.

## Appendix C: Numerical validation

Complementary to the previous convergence analysis, which demonstrated the accuracy of the HOS numerical wave tank compared to the experiments on the free surface elevation at a single location, this paragraph assesses its relevance for the systematic study presented in Sections 3 and 4. The experimental data listed in Table 5 are further used to validate the accuracy of the established numerical model for long-time simulation and statistical characteristics.

The total simulation duration is 2000 s, i.e., more than 1000 waves are involved in the simulated wave trains. 2D HOS numerical model is adopted to reproduce the two wave trains. Figure 18 compares the free surface elevations between the

<figure><img src="image_25.png" /><figcaption>**Fig. 16** Convergence of the numerical simulation with respect to the number of points per corresponding peak wavelength</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. 17** Comparison of amplitude spectra input target and numerically generated at x = 5 m</figcaption></figure>

experimental data and numerical results at some pre-setting locations along the wave tank. The horizontal axis represents the time corrected from the group velocity at the peak of the spectrum $c_g$. This allows to follow the wave groups in their evolution and possibly identify the large wave events in each

case. It can be noted that all the free surface elevations of the numerical simulation remain consistent with the experimental data along the wave tank, irrespective of the value of BFI.

The evolution of the kurtosis along the wave flume for both cases with different BFI values is presented in Fig. 19, where experiments and numerical results are compared. The horizontal axis represents the distance away from the wave-maker normalized by the corresponding peak wavelength. Complementary, the probability distribution of the wave height at different positions of the numerical waves is compared to the experimental data in Fig. 20, in which Naess distribution are plotted as a reference. It can be found that the kurtosis and the wave height distribution of the numerical results have good agreement with the experimental data, even for case II that exhibits the largest waves. These comparisons in Figs. 18, 19, and 20 demonstrate that the established HOS numerical model is able to accurately simulate the irregular wave trains for a long time, even in the presence of significant modulation instability. Furthermore, in Fig. 20, the wave height distributions show different behaviors for the two cases at different locations. This is an awaited behavior from literature: lower $k_p h$ and lower BFI induce weaker modulation instability (Janssen 2003; Janssen and Bidlot 2009; Fedele 2015) and consequently smaller occurrence of extreme

<figure><img src="image_27.png" /><figcaption>**Fig. 18** The comparisons of the free surface elevations between experimental data and numerical results along the wave tank (the black solid lines represent the experimental data, and red solid lines represent the numerical results). a Case I. b Case II</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>**Fig. 19** The comparisons of kurtosis between experimental data and numerical results along the wave tank ($L_p$ is the wavelength corresponding to the peak frequency). a Case I. b Case II</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 20** The comparisons of the probability distribution of the wave height between experimental data and numerical results at different locations ($L_p$ is the wavelength corresponding to the peak frequency). a Case I. b Case II</figcaption></figure>
waves (characterized by a smaller kurtosis). We are consequent-
ly confident in the accuracy of the numerical model to study this
phenomenon.

References

Boccotti P (2000) Wave mechanics for ocean engineering. Elsevier
Science, Oxford

Bonnefoy F, Ducrozet G, Le Touzé D, Ferrant P (2010) Time domain simulation of nonlinear water waves using spectral methods. Adv Numer Simul Nonlinear Water Waves:129–164

Dommermuth DG, Yue DKP (1987) A high-order spectral method for the study of nonlinear gravity waves. J Fluid Mech 184:267–288

Ducrozet G, Gouin M (2017) Influence of varying bathymetry in rogue wave occurrence within unidirectional and directional sea-states. J Ocean Eng Mar Energy 3(4):309–324

Ducrozet G, Bonnefoy F, Le Touzé D, Ferrant P (2012) A modified high-order spectral method for wavemaker modeling in a numerical wave tank. Eur J Mech-B/Fluids 34:19–34

Fedele F (2015) On the kurtosis of deep-water gravity waves. J Fluid Mech 782(7):25–36

Fedele F, Brennan J, Ponce de León S, Dudley J, Dias F (2016) Real world ocean rogue waves explained without the modulational insta-bility. Sci Rep 6:27715

Fernandez L, Onorato M, Monbaliu J, Toffoli A (2014) Modulational instability and wave amplification in finite water depth. Nat Hazards Earth Syst Sci 14(3):705–711

Fernandez L, Onorato M, Monbaliu J, Toffoli A (2016) Occurrence of extreme waves in finite water depth. Extreme Ocean Waves. Springer, Cham, pp 45–62

Forristall GZ (1984) The distribution of measured and simulated wave heights as a function of spectral shape. J Geophys Res 89(C6):547–552

Forristall GZ (2000) Wave crest distributions: observations and second-order theory. J Phys Oceanogr 30(8):1931–1943

Gramstad O, Trulsen K (2010) Can swell increase the number of freak waves in a wind-sea? J Fluid Mech 650:57–79

Guedes Soares C (1984) Representation of double-peaked sea wave spectra. Ocean Eng 11(2):185–207

Guedes Soares C (1991) On the occurrence of double peaked wave spectra. Ocean Eng 18(1–2):167–171

Janssen PAEM (2003) Nonlinear four-wave interactions and freak waves. J Phys Oceanogr 33(4):863–884

Janssen PAEM, Bidlot J-R (2009) On the extension of the freak wave warning system and its verification. Technical Memorandum, 588. ECMWF, Reading, pp 42

Janssen PAEM, Onorato M (2007) The intermediate water depth limit of the Zakharov equation and consequences for wave prediction. J Phys Oceanogr 37(10):2389–2400

Klinting P, Sand SE (1987) Analysis of prototype freak waves. Coastal Hydrodynamics, ASCE, 618–632

Li JX, Liu SX (2015) Focused wave properties based on a high order spectral method with a non-periodic boundary. China Ocean Eng 29(1):1–16

Li JX, Li PF, Liu SX (2013) Observations of freak waves in random wave field in 2D experimental wave flume. China Ocean Eng 27(5):659–670

Longuet-Higgins MS (1952) On the statistical distribution of the heights of sea waves. J Mar Res 11:245–266

Mori N, Janssen PAEM (2006) On kurtosis and occurrence probability of freak waves. J Phys Oceanogr 36(7):1471–1483

Mori N, Yasuda T (2002) A weakly non-Gaussian model of wave height distribution for random wave train. Ocean Eng 29(10):1219–1231

Mori N, Onorato M, Janssen PAEM, Osborne AR, Serio M (2007) On the extreme statistics of long-crested deep water waves: theory and experiments. J Geophys Res 112(C9):C09011

Naess A (1985) On the statistical distribution of crest to trough wave heights. Ocean Eng 12:221–234

Ochi MK, Hubble EN (1976) Six-parameter wave spectra. Proceeding of 15th Coastal Engineering Conference, 301-328

Onorato M, Osborne AR, Serio M, Cavaleri L, Brandini C, Stansberg CT (2004) Observation of strongly non-Gaussian statistics for random sea surface gravity waves in wave flume experiments. *Phys Rev E* 70(6):067302

Onorato M, Osborne AR, Serio M, Cavaleri L (2005) Modulation insta-bility and non Gaussian statistics in experimental random water-wave trains. Phys Fluids 17(7):1–4

Onorato M, Osborne AR, Serio M (2006) Modulational instability in crossing sea states: a possible mechanism for the formation of freak waves. Phys Rev Lett 96(1):014503

Onorato M, Cavaleri L, Fouques S, Gramstad O, Janssen PAEM, Monbaliu J, Osborne AR, Pakozdi C, Serio M, Stansberg C, Toffoli A, Trulsen K (2009) Statistical properties of mechanically generated surface gravity waves: a laboratory experiment in a three-dimensional wave basin. J Fluid Mech 627:235–257

Petrova PG, Guedes Soares C (2009) Probability distributions of wave heights in bimodal seas in an offshore basin. Appl Ocean Res 31(2): 90–100

Petrova PG, Guedes Soares C (2011) Wave height distributions in bimodal sea states from offshore basins. Ocean Eng 38(4):658–672

Regev A, Agnon Y, Stiassnie M, Gramstad O (2008) Sea-swell interaction as a mechanism for the generation of freak waves. Phys Fluids 20(11):112102

Rice SO (1944) Mathematical analysis of random noise. Bell Syst Tech J 23(3):282–332

Rodriguez G, Guedes Soares C, Pacheco M, Perez-Martell E (2002) Wave height distribution in mixed sea states. J Offshore Mech Arct Eng 124:34–40

Shemer L, Sergeeva A (2009) An experimental study of spatial evolution of statistical parameters in a unidirectional narrow-banded random wavefield. J Geophys Res 114(C1)

Slunyaev AV (2005) A high-order nonlinear envelope equation for gravity waves in finite-depth water. J Exp Theor Phys 101(5):926–941

Stansell P (2004) Distributions of freak wave heights measured in the North Sea. Appl Ocean Res 26(1–2):35–48

Støle-Hentschel S, Trulsen K, Nieto Borge JC, Olluri S (2020) Extreme wave statistics in combined and partitioned windsea and swell. Water Waves

Tayfun MA (1981) Distribution of crest-to-trough wave heights. *J Waterw Port Coast Ocean Eng* 107(3):149

Tayfun MA (1990) Distribution of large wave heights. J Waterw Port Coast Ocean Eng 116(6):686–707

Tayfun MA (2006) Statistics of nonlinear wave crests and groups. Ocean Eng 33(11–12):1589–1622

Tayfun MA, Fedele F (2007) Wave-height distributions and nonlinear effects. Ocean Eng 34(11–12):1631–1649

Toffoli A, Benoit M, Onorato M, Bitner-Gregersen EM (2009) The effect of third-order nonlinearity on statistical properties of random directional waves in finite depth. Nonlinear Process Geophys 16:131–139

Toffoli A, Bitner-Gregersen EM, Osborne AR, Serio M, Monbaliu J, Onorato M (2011) Extreme waves in random crossing seas: laboratory experiments and numerical simulations. Geophys Res Lett 38: L06605

West BJ, Brueckner KA, Janda RS, Milder DM, Milton RL (1987) A new numerical method for surface hydrodynamics. J Geophys Res 92(C11):11803–11824