Extreme Wave Statistics in Combined and Partitioned Windsea and Swell

We investigate how the extreme wave statistics of a combined windsea and swell appears to be different from the extreme wave statistics of the corresponding windsea and swell partitions. We consider the situation of following long-crested windsea and swell in laboratory experiments and in simulations using a high-order spectral method (HOSM). We also consider the situation of short-crested windsea and swell crossing at nearly right angle, corresponding to the sea state when the Prestige accident happened, in hindcast simulations combined with HOSM. For cases when the two wave systems do not interact much, the combined wave system appears to be more Gaussian than the corresponding partitioned wave systems, consistent with the central limit theorem. This result is found for kurtosis and exceedance probability of envelope and crest height.


Introduction
Multiple wave systems constitute 15-25% of sea states observed in different locations around the world [3,20,21]. They may be distinguished by their interaction angle that is in the range of [ [42]. The analysis of 270 ship accidents, Dedicated to Kristian B. Dysthe on the occasion of his 80th birthday. performed by Toffoli et al. [42], showed that most of them occurred in crossing seas and were associated with a rapid change of the sea state.
Examples of accidents that happened in crossing seas are the Louis Majesty incident [4] that occurred at a crossing angle of 40 • -50 • , the Prestige accident [45] at a crossing angle of 90 • , and the El-Faro accident [13] that occurred in a sea state multimodal both in direction and frequency. Crossing seas were also claimed as an important condition preceding the Suwa-Maru incident [39]. The latter happened in an 'extremely narrow' [39] unimodal swell that had developed from two swells and one windsea.
The hypothesis of increased freak wave probability due to mutual interaction of coexisting wave systems was studied for the special case of crossing seas with identical peak periods by Onorato et al. [33] and Toffoli et al. [43]. Onorato et al. [33] investigated crossing seas theoretically by coupled nonlinear Schrödinger (CNLS) equations and found larger regions of instability and large growth rate for coupled systems compared to single systems. Simulations of the CNLS were then compared with simulations of the higher order spectral method (HOSM) and laboratory measurements. The kurtosis was found to be increased for two systems with an interaction angle between 40 • and 60 • [43], indicating enhanced probability of freak waves according to Mori and Janssen [28]. Generalization into crossing wave systems that do not have the same peak period was done by Gramstad et al. [18].
There are two opposite regimes for the interaction between windsea and swell. For interactions over a long scale, there is evidence that a windsea may modify a swell [27] and may enhance the occurrence of freak swell waves [39] through resonant transfer of energy from the windsea to the swell, such that the spectrum may become narrower. This may happen when the separation in wave periods between the two systems is small, making it unfeasible to partition the combined system into two separate systems. Some ship accidents indeed seem to have occurred in conditions of narrowing wave spectra, e.g. Tamura et al. [39], Waseda et al. [47,48].
For the opposite regime of short interaction scale, the swell can modulate the amplitude and phase of the short waves. Regev et al. [37] found that a windsea and a weak swell crossing at right angle can suffer wave modulations that lead to freak waves. Similarly, Gramstad and Trulsen [16] computed the modification of the probability of freak waves in a windsea perturbed by a weak swell crossing at various angles. They found that the swell could enhance the occurrence of freak waves in the windsea slightly. The increase was found to be minimal in the case of right angle between the swell and the windsea. This regime corresponds to large separation between the two wave systems such that a partitioning is both feasible and desirable. The focus herein is the latter regime, in particular if the resulting modification of the nonlinear wave evolution can make extreme or freak waves more probable.
Following Longuet-Higgins and Stewart [26], a large number of studies were devoted to how longer waves modulate shorter waves, e.g. Longuet-Higgins [25], Craik [6], Grimshaw [19], Henyey et al. [23], Zhang and Melville [49], Zhang and Melville [50]. The amplitudes of the short waves are enlarged and their wavelength is shortened close to the crests of the underlying long waves. Near the troughs of the long waves, the short waves obtain longer wavelengths and smaller amplitudes. Zhang and Melville [50] concluded that the modulational instability of weakly nonlinear short waves was enhanced when they were riding on finite-amplitude long waves. Similar results for the nonlinear evolution of short gravity waves on long waves were also obtained by Naciri and Mei [29][30][31]. First they studied short gravity waves on long rotational Gerstner waves [29], then irrotational short waves on irrotational long waves [30] and then two-dimensional interaction of obliquely intersecting waves [31]. For the latter, the instability of the short wave due to oblique side bands was shown to be enhanced by the presence of the long wave. However, the obliqueness becomes important only when the steepness of the long wave is sufficiently large.
The work presented herein concerns the wave statistics of different bi-modal sea states that are combinations of swell and windsea. Mixed seas have been studied by Petrova and Guedes Soares [34], Petrova and Guedes Soares [35], Petrova and Guedes Soares [36] and Guedes Soares and Carvalho [22] with focus on wave height distributions. For swell dominated seas, they found that skewness and kurtosis are lower than in windsea dominated seas [34,36]. They also found that the exceedance probability of windsea dominated cases have thicker tails than those of swell dominated cases [34][35][36]. Within each of the modes, swell or windsea, the wave systems are usually multi-directional or short-crested. Important contributions on how to correctly account for kurtosis in short-crested sea was given by Fedele [9] as an improvement over the theory of Mori and Janssen [28], and applications were presented in Fedele and Tayfun [10], Fedele et al. [11], Fedele et al. [12].
We suggest the importance of partitioning a well-separated bi-modal sea state of swell and windsea into the swell alone and the windsea alone, and to carry out analysis of each of the partitioned systems as well as the combined system. In Sect. 2 we describe a laboratory experiment in which we generated a unidirectional (long-crested) bimodal sea state of following swell and windsea, and in Sect. 3 we describe numerical simulations of the same wave configuration. In Sect. 4 we summarize the characteristics of the sea states considered in the experiments and the numerical simulations, we present their kurtosis in Sect. 5, their envelope in Sect. 6, and their crest height in Sect. 7. In Sect. 8 we consider the multi-directional (short-crested) bi-modal sea state during the Prestige accident, with swell and windsea crossing at right angle [45].
Our main conclusion is that the extreme wave statistics of well-separated multimodal sea states can be better comprehended if analysis is done on the partitioned wave systems as well as the combined sea state.

Setup of the Laboratory Experiments
Experiments with unidirectional bi-modal waves were performed in a wave flume at the Department of Mathematics at the University of Oslo. The dimensions of the tank and the setup are sketched in Fig. 1. The wave tank is 24.6 m long, 0.5 m wide and was filled to a mean water depth h = 0.7 m. Unidirectional waves were mechanically generated at one end of the tank by a hydraulic piston and were subsequently damped at the other end of the tank by a 3 m-long beach to prevent reflection. The surface elevation was measured with ultrasonic probes at the four positions indicated in Fig. 1. Key parameters for the laboratory experiment are listed in Table 1. Visual inspection  during the experiments did not reveal wave breaking, although the measurements do indicate dissipation. Waves were generated according to a Torsethaugen spectrum given in its simplified version for windsea dominated spectra [46]: The values for the peak periods and significant wave heights are summarized in Table 2.
Parameters for windsea and swell are distinguished by the corresponding subindex W and S. We have renamed the shape factor γ in Veritas [46] to ν to avoid collision with the skewness γ . Similarly, we use α instead of σ in Veritas [46], since we will use σ for standard deviation. We believe the peak wavelength of the windsea is long enough to suppress crosstank modulations provoked by modulational instability, the parameter μ of Trulsen et al. [44] having the value μ = 0.4 0.1 for the peak period which should ensure we are outside the domain of transversal modulationally unstable modes.

Setup of the Numerical Wave Tank
Numerical simulations of unidirectional bi-modal waves were carried out with a higher order spectral method, the HOS-NWT numerical wave tank described by Bonnefoy et al. [2] and Ducrozet et al. [7]. We applied a 1D wave tank with the same dimensions as in our laboratory, spatially discretized at 1024 points. The water depth was set to five meters. A numerical wave paddle at one end of the tank was programmed to generate a two-peak Ochi and Hubble [32] spectrum This is a superposition of a windsea spectrum, with index W , and a swell spectrum, with index S. Both spectra are defined by the same formula depending on the significant wave height H s j , the peak period T p j , the shape factor ν j which was set to 3.3 for swell and windsea. Γ is the Gamma function. The Ochi and Hubble spectrum was preferred for the simulations since the swell and the windsea are better separated than in the Torsethaugen spectrum, thus facilitating the subsequent analysis. Unfortunately, since the experiments were carried out a long time before the numerical simulations and the subsequent analysis, we have employed different spectra in the experiments and the simulations. The performance characteristics of the paddle were set to the default values. The damping beach was simulated by the inherent numerical beach with an absorption coefficient of unity [2,7].
Key parameters for the numerical simulations are listed in Table 1.

Characterization of the Sea States
In the simulations the spectrum develops fast in the vicinity of the wavemaker, and then reaches an approximately steady shape beyond five peak wavelengths from the wavemaker. Figure 2 shows an example of the spectra at the position x = 23 λ p where λ p is the peak wavelength of the windsea. The vertical dotted lines mark the splitting frequencies ω m that were employed to partition the combined swell and windsea (SW ) into swell (S) and windsea (W ). Torsethaugen and Ochi and Hubble spectra (a) (b) Fig. 2 Scaled spectra for the different laboratory experiments and simulations. The spectra are distinguished by the sea swell energy ratio (SSER) as defined by Guedes Soares [20] require different splitting frequencies, resulting in ω m = 0.685 ω p for laboratory and ω m = 0.75 ω p for simulation, where ω p is the peak frequency of the windsea. The different sea states are labeled by the sea-swell energy ratio (SSER) as suggested by Guedes Soares [20], based on the zeroth spectral moments for windsea, m 0W , and swell m 0S .
The key technique employed in the following is the partitioning of the combined wave system into two non-overlapping systems. The windsea was extracted by highpass filter (HP) and the swell by lowpass filter (LP). A similar technique was employed by, e.g. Trulsen et al. [45] and Støle-Hentschel et al. [38].
Since the windsea spectrum is equal for all the cases, its spectral peak S p = S W at the peak frequency ω p = ω pW serves as scaling parameters and the indices for windsea are dropped for these two parameters.
The wave steepness and the Benjamin-Feir index (BFI) are characteristics of the sea states that can only be estimated for individual wave systems [38,45]. The wave steepness definition employed here is ε = k p a c , where a c = H s / √ 8 is a characteristic amplitude and where H s is the significant wave height defined as four times the standard deviation of the surface elevation. We only consider the BFI of the windsea BFI W = ε W ω p /Δ W , where Δ W is the 'half width at half peak'. Table 2 summarizes the characteristics of the wave system employed in this study.
The splitting frequency ω m and the steepnesses are given in terms of the peak frequency ω p and the peak wavelength λ p of the windsea (a) (b) Fig. 3 The correlation coefficient between highpass filtered windsea and highpass filtered mixed cases with the indicated SSER As a first estimate of the degree of interaction between swell and windsea we calculate the correlation coefficient of pure windsea and the windsea part of a combined windsea and swell sea state: where η W,HP and η SW,HP are the surface elevations of pure windsea and combined swell and windsea both filtered by a high-pass filter. The corresponding standard deviations are denoted σ W,HP and σ SW,HP . In the case of the simulations, the correlation coefficient shown in Fig. 3b was averaged over all simulations. Close to the wavemaker, the pure windsea and the windsea parts of the mixed seas are practically identical for all cases. In the development along the tank the correlation coefficient decays depending on the strength of the swell. Close to the damping beach, the correlation coefficient decays to a level of 0.85 for the laboratory and 0.75 for the simulation for the strongest swell. These observations suggest that the swell has only marginal influence on the phases of the windsea over the given evolution length. Figure 4 presents the laboratory measurements and the simulation results in terms of H s and kurtosis for the highpass filtered windsea. In the case of the laboratory experiments the local kurtosis, κ, was computed for each fixed position along the tank. The damping beach in the laboratory is too short for the longest waves leading to seiche modes. These have been filtered by a highpass filter for the entire data analysis. For the simulations, we have employed a lower time resolution of Δt = 0.14 s and the part of the simulation time subject to analysis is limited to 77 T p . To ensure sufficient surface elevation measurements for adequate estimates of kurtosis 32 neighboring grid points were considered as one evaluation point in space, and 10 simulations were merged; thus the mean values and the confidence intervals are based on 50 simulations arranged into five groups. The strong decay of H s in the laboratory is attributed to dissipation that has a strong effect on the experiment due to the short peak wavelength of the windsea spectrum. This reduction in wave steepness affects the kurtosis as reflected in the reference curve based on Mori and Janssen [28]. For the simulations the significant wave height is constant along the tank. Overall the difference in kurtosis for the different spectra is smaller than the variation of the values. In Fig. 5 we analyze the local kurtosis for the waves of the combined spectra. Both the laboratory and the simulation results show that the SSER is correlated to the kurtosis of combined seas.

Kurtosis
From the definition of kurtosis we calculate reference values for kurtosis of combined sea states based on the assumption that the two systems are independent. From the analysis of the correlation coefficient, we anticipate that this may be a good approximation. We find the formula with {σ 2 S , κ S } the variance and kurtosis of the swell partition, and {σ 2 W , κ W } the variance and kurtosis of the windsea partition. Figure 5 shows that the statistics of the combined wave system differs clearly from the statistics of the windsea shown in Fig. 4. The estimate in Eq. (5) is similar but not equal to the calculated values.
The kurtosis of the windsea is generally high as expected for long-crested waves. According to Gramstad and Trulsen [15], wave crests longer than ten times the wavelength are associated with high kurtosis, even for large spectral bandwidths that are not typically associated with modulational instability according to a low Benjamin-Feir index [24].

Exceedance Probability of Wave Envelope
In Fig. 6 we show the exceedance distribution for the Hilbert wave envelope computed according to Fig. 6 Exceedance probability of Hilbert envelopes scaled by the standard deviation of the surface elevation at one position given in terms of the peak wavelength of the windsea. The mixed sea is based on the maximum swell for laboratory (SSER = 3.9) and simulation (SSER = 0.9), respectively. The windsea has a SSER of 62.8 for the laboratory and 17.0 for the simulation whereη is the Hilbert transform of the surface elevation η. The envelopes were scaled by the standard deviation of the corresponding surface elevation. For reference, we compare the exceedance distributions with the Rayleigh distribution corresponding to a Gaussian process. In this and the next section we present pure windsea and the combined sea with the strongest swell for laboratory and simulation. Pure swell simulations were added for comparison. The sea states are analysed both filtered and unfiltered. The results in Fig. 6 are clearer for the simulations than for the laboratory experiments. The swell envelopes of all configurations coincide with the Rayleigh distribution in good correspondence with the fact that the swell is essentially linear. While the envelopes of the combined wave fields differ, their windsea contributions are similar and deviate clearly from Gaussian behavior.

Exceedance Probability of Crest Height
In Fig. 7 we show the exceedance probability of the crest height η c scaled by the standard deviation σ of the corresponding surface elevation. We only consider the exceedance probability of crests for simulations since the number of measured crests in the laboratory was too small for meaningful analysis. For better resolution of extreme waves, the number of simulations was increased to 180, compared to that of Table 1. The simulated data were interpolated by truncated Fourier series expansion in time before determining the crest height as maximum between two zero crossings. Below the horizontal dotted lines the number of crests is lower than one hundred which we consider too few for statistical significance.
The exceedance probability is characterized by local variation, similar to what we previously found for the kurtosis. To illustrate the amount of local variation, we show the exceedance probability at four locations along the numerical wavetank in Fig. 7. When comparing the waves after partitioning into swell and windsea, the statistics are approximately identical irrespective if the other wave system is present or not.
We have employed the Rayleigh, Tayfun [40] and Tayfun-Fedele [12,41] distributions as reference. The crests of narrowbanded linear waves are Rayleigh distributed. Although the simulated swell is not narrowbanded, its crests are found to be approximately Rayleigh distributed. The Tayfun distribution, here based on the standard deviation of the windsea, is a reference distribution for narrowbanded second-order wave crests. The windsea distributions clearly exceed Tayfun even though they are clearly not narrowbanded. Hence we anticipate that nonlinearities of third order must play an important role for the windsea. As an approximation of a third-order distribution of crest heights we employ the Tayfun-Fedele distribution. It is a continuation of the second-order Tayfun distribution by the Gram-Charlier series up to the fourth cumulant [12,41]. The distribution changes its form with the wave steepness by its dependence on the skewness and the occurrence of extreme values by the incorporation of the kurtosis. The Tayfun-Fedele distribution shows good agreement with the simulated windsea exceedance probability of crest heights. Figure 8 reveals that the combined swell and windsea results in an exceedance probability for crest heights between that of pure swell and that of pure windsea, rather similar to narrowbanded waves of second order. However, the combined exceedance probability is based on a linear swell and a third-order windsea. Linear and secondorder nonlinear [14] simulations were added.

The Short-Crested Crossing Sea State of the Prestige Accident
In Trulsen et al. [45] we considered the bi-modal and multi-directional (short-crested) sea state in which the oil tanker MV Prestige accident occurred on 13 November 2002. The tanker suffered damage on the starboard side of its hull while sailing close to Cape Finisterre, off the northwestern coast of Spain.
In the previous hours of the accident, the windsea was evolving towards higher periods and the swell in the opposite way. The final situation after the accident was a swell with higher mean period and a windsea with lower mean period, with a similar value of dominant mean period T m about 9 s for both wave systems, both systems with a comparable significant waveheight (about 4.5 m) and crossing at about 90 • .
In Trulsen et al. [45] the global wind-forcing model ARPEGE [1,5] was used to force a ECWAM-based wind generation model [8] to provide hourly directional hindcast spectra. Those spectra were used as input of HOSM to simulate the temporal evolution in terms of sea surface wave fields. For each initial hindcast spectrum, a total of 50 realizations with the initial random phases between [0, 2π) were performed. The temporal evolution of the wave fields covered a total time of 100 T m , with a sampling interval of 10 T m . We here reconsider those simulations in order to check the performance of Eq. (5) on the combined and the partitioned swell and windsea. Figure 9 shows the kurtosis of the combined sea state in black, of the windsea partition in blue, of the swell partition in green, and of the reconstruction of the kurtosis according to Eq. (5) in red. The mismatch between the kurtosis of the combined sea state and the reconstruction according to Eq. (5) is likely due to nonlinear interaction between the swell and the windsea. During the first 4 h of the simulation, the red and the black curve are close and we may anticipate little interaction between the systems. During this period of time, the kurtosis of the combined sea state remains clearly below the windsea partition. At later times the increasing kurtosis of the swell appears to cause an interaction between the swell and the windsea increasing the extreme wave statistics of the combined sea state.

Discussion
It is known that describing a unified steepness for multi-modal sea states is difficult [16,17,29,30,38,45]. Here we have shown that the partitioning of the spectrum may not only be beneficial for estimating the various partial steepnesses, but also for better understanding of extreme wave statistics.
Equation (5) shows that the kurtosis of a combined bi-modal sea, where the modes are statistically independent, is a weighted combination of the kurtosis of the partitioned systems. It can be seen that Eq. (5) is in agreement with the central limit theorem; it predicts a combined kurtosis closer to a Gaussian than the partitioned systems. However, if the two systems interact, we can expect higher values of combined kurtosis as seen in the Prestige data. In general, a consideration of only the combined kurtosis could conceal the non-Gaussian behavior of a nonlinear windsea.
A similar behavior seems to occur for the exceedance probabilities. In our unidirectional bi-modal cases the swell affects the windsea so little that the two systems may be regarded as independent; however, a consideration of only the extreme wave statistics for the combined sea may conceal how dangerous the windsea partition really is.

Conclusion
We have shown that extreme wave statistics of a long-crested windsea can be nearly unaffected by the presence of a following swell; at the same time as the extreme wave statistics of the combined following long-crested windsea and swell appears to be more Gaussian and thus less extreme. We have also shown evidence that analysis of partitioned wave systems can give additional insight, reveal nonlinear interaction between the wave systems, and suggest if one of the systems should be studied on its own right. Consistent with the central limit theorem, superposition between independent partial wave systems should converge to a Gaussian sea, thus possibly concealing strong non-Gaussian behavior in one of the partitions. Consideration of a combined sea alone may undermine how dangerous a windsea partition really is. Therefore, we believe that mixed sea states of windsea and swell should be partitioned, and analysis be performed on each partition, in order to better interpret the combined extreme wave statistics.