Large waves and navigation hazards of the Eastern Mediterranean Sea.

We present a statistical analysis of deep-water buoy measurements of large waves generated during two major storms of the Eastern Mediterranean in 2017 and 2018, respectively. The largest waves observed do display similar characteristics to those of the Draupner, Andrea, and El Faro rogue waves in that second order bound nonlinearities enhance the linear dispersive focusing of extreme waves. We also present a novel analysis of waves in space-time to predict potential risks posed by such large waves to navigation. In particular, we consider the scenario of two types of vessels of the Israeli Navy fleet navigating during the most intense stages of the two storms considered here and provide predictions for the largest waves likely to be encountered.

Introduction

In light of global warming, weather events are likely to become more severe and hazardous to coastal and offshore structures and maritime navigation. Not only does this relate to vast oceans, but also to large basins such as the Black Sea and the Mediterranean Sea. In particular, the Eastern Mediterranean deserves special attention because it lacks studies on the physical characteristics of extreme sea states and large waves. A climatology hot-spot and natural treasures area, the Eastern Mediterranean is a basin of increasing interest, both ecologically and economically[1]. Thus, a scientific investigation and prediction of wave extremes are of paramount importance, especially for rogue wave predictions and maritime navigation in extreme weather.

Conventional thresholds[2] have been adopted to distinguish if a wave is rogue or not. In particular, a wave is rogue if its crest height  or the crest-to-trough (wave) height , where is the significant wave height of the sea state in which the rogue wave occurred. For example, the first famous rogue wave ever recorded is the Draupner wave that occurred on the of January 1995 in an intense sea state with  m[3, 4]. The 26-m rogue wave hit the lower deck of the “Draupner S” oil platform at the North Sea. The Draupner wave has and satisfying the conventional criteria for being rogue[2]. On November 2007, the powerful Andrea storm crossed the “Ekofisk” oil platform. In a sea state with  m four wave gauge lasers mounted on the rig measured the Andrea rogue wave[5], whose wave height and crest height . In November 2020 a rogue wave has been detected off the Vancouver Island in the Atlantic Ocean[6], whose wave height and crest height are well over the conventional rogue thresholds[2]. Nevertheless, the crest height of that rogue wave reached 11.96 m in amplitude, and thus smaller than both the Draupner and Andrea waves.

There is a lack of scientific knowledge about the metocean characteristics of the eastern part of the Mediterranean Sea, especially for rogue wave predictions and vessel navigation in extreme weather. Reports of rogue wave occurrences involving vessel navigation have been mostly observed in the western part of the Mediterranean by eyewitnesses on cruise ships that have been impacted by extreme waves resulting in structural damage and loss of lives. Two cases have been recently reported: the “Voyager” accident that occurred in 2005 and the “Louis Majesty” accident in 2010. Both accidents triggered scientific interest to confirm, or rebut the possibility of cruise ships encountering rogue waves[7, 8]. A rogue wave event has been caught in the Northern Adriatic Sea by a web-camera placed on the oceanographic tower ’Acqua Alta’ off the Venice coast, Italy. The estimated wave crest height is in the range of 5.1–6.4 m[9] and exceeded the conventional rogue thresholds[2]. All the above-mentioned extreme waves have been observed in the western part of the Mediterranean Sea.

Recently, Knobler et al.[10] provided the first observations of several rogue wave events in the Eastern part of the Mediterranean Sea from measurements acquired by AXYS Technologies “3 metre” metadata buoy in deep waters (see “Methods” section for details). Surface elevations, wind and wave parameters were also collected by a meteorological station[10, 11]. As expected for this region, the strongest storm events occurred throughout October–April period[10]. They detected a sample of 109 rogue waves, whose wave heights exceeded the rogue threshold , and only nine waves out of the sample had crests exceeding in height[2]. However, no wave crest exceeded 4 m in height.

In the present work, we further examine the statistics of large waves in deep waters observed during the two most powerful storm events of the Eastern Mediterranean recorded during the last five years. In particular, the first storm occurred during January 16–22 in 2018, and hereafter referred to as Gaia. The second storm occurred during January 25–31 in 2017, and hereafter referred to as Jan 2017. We note that the largest waves observed in both storms do not exceed the rogue wave thresholds[2], in contrast to the Draupner and Andrea rogue waves[3–5, 12, 13] that occurred in intermediate water depths and did reach crest heights of  well over the rogue threshold. In our studies, we also exploit statistical distributions for non-stationary time series[14] to analyze the statistics of surface elevations sampled from the two storms. Finally, we carry out a novel space-time analysis of the most intense sea states as a basis for predicting potential rogue hazards for ships navigating through the rough waters of the storms.

Results

This section is structured as follows. First, we discuss the metocean characteristics of the sea states generated by the two observed storms as they passed through the Eastern Mediterranean. Relevant wave parameters and statistical models are defined in the “Methods” section. Then, we make prediction on the extreme waves during the two storms. To do so, we draw on Fedele et al.[14] and use their formulation of the probability structure of non-stationary time series as a model for the sequence of sea states in a storm. In particular, we study the occurrence frequencies of rogue waves during Jan 2017 and Gaia storms encountered by an observer at a given point of the sea surface. The wave extremes observed at the peak of the storms and their characteristics are then compared to those of the Andrea and Draupner rogue waves, observed at different oil platforms in the North Sea in 1995 and 2007, respectively[12] as well as the simulated El Faro rogue wave[13]. The metocean parameters of these five sea states are summarized in Table 1. Finally, we present a space-time wave analysis[13, 15, 16] of the most intense sea states to predict rogue hazards for vessel navigation within an area.

Table 1

Wave parameters and various statistics of Jan 2017 and Gaia at the storm peak in comparison to the El Faro, Andrea and Draupner rogue sea states[12].

	El Faro	Andrea	Draupner	Jan 2017	Gaia
Significant wave height H_s [m]	9.0	10.0	11.2	6.86	6.16
Dominant wave period T_p [s]	10.2	14.3	15.0	11.9	14.1
Mean zero-crossing wave period T_0 [s]	9.2	11.6	11.9	8.6	9.2
Mean wavelength L_0 [m]	131	209	219	115	131
Depth d [m], k_0d with k_0=2π/L_0	4700, 2.25	74, 2.23	70, 2.01	258,14.11	258, 12.36
Spectral bandwidth ν	0.49	0.35	0.36	0.43	0.46
Angular spreading σ_θ[rad]	0.79	0.37	0.39	1.08	1.11
Parameter R=σθ2/2ν^2	1.34	0.56	0.59	3.18	2.88
Benjamin Feir Index BFI in deep water[17]	0.36	0.24	0.23	0.16	0.09
Tayfun  NB skewness λ_3,NB[18]	0.262	0.159	0.165	0.222	0.199
Observed skewness λ_3	0.162	0.141	0.146	0.103	0.105
Maximum NB dynamic excess kurtosis λ40,maxd[19]	10^-3	2.3×10^-3	2.1×10^-3	-4.8×10^-3	-1.6×10^-3
Janssen NB bound excess kurtosis λ40,NBd [20]	0.049	0.065	0.074	0.097	0.087
Observed excess kurtosis λ_40	0.042	0.041	0.032	0.100	0.087
Actual maximum crest height h/H_s	1.68	1.63	1.55	1.04	1.11
Actual maximum crest-to-trough (wave) height H/H_s	2.6	2.30	2.10	1.98	1.96

We refer to the “Methods” section for the definitions of wave parameters.

Our extreme wave analysis focused on the study of the time sequence of changing sea states during the two storms. While the analysis of Fedele et al.[14] yields an optimal sea state duration of 40–50 min to minimize the variation between waves of consecutive sea states of storms, the duration T of our actual measurements was limited to 18 min at the beginning of each hour to minimize data storage. So, the analyses here were unavoidably based on  min. Nevertheless, we observed that our wave statistics are very robust to changes in up to  min, provided that surface elevations of each sea state in the storm sequence are normalized by the respective significant wave height.

Metocean parameters

The metocean parameters history of Jan 2017 storm are shown in Fig. 1. In particular, the top panel depicts the hourly variation of the significant wave height. This is estimated as (black line). For comparison, is also shown (red line), which is smaller than . The mean zero-up-crossing wave period is shown in the second panel of the same figure, and the hourly variations of are reported in the third panel, where is the peak wavenumber. The bottom panel on the left depicts the frequency spectrum measured at the storm peak. The directional wave spectrum is also shown, in the bottom panel on the right. Similarly, the histories of metocean parameters of storm Gaia are presented in Fig. 2. The two storms generated waves in deep waters depth (, see[14]). Moreover, is well above the critical threshold 1.363, at which modulation instability disappears. Nevertheless, third-order nonlinear modulational effects are negligible because the energy of the two sea states spreads both in all directions and frequencies[12, 19] as clearly seen in the bottom panels of Figs. 1, 2. Furthermore, we depict in Fig. 3 the Boccotti parameter , which measures the spectral bandwidth of the sea states of the two storms. The values of the observed vary in the range [0.5–0.7] indicating sea states with broadband spectra. In Table 1 we report the metocean parameters of the two storm events at their peak stages in comparison to those of the El Faro, Draupner and Andrea rogue sea states[12, 13]. Note that the five sea states all have similar metocean characteristics. Both Jan 2017 and Gaia are in deep waters and developed comparable steep sea states as the other three sea states, even if their peak significant wave height is approximately smaller than that of the other rogue sea states. Indeed, the observed wave parameters have comparable values in the same range typical of oceanic sea states in storms. Our analysis indicates that the largest observed waves in Gaia and Jan 2017 storms have characteristics quite similar to those displayed by the El Faro, Andrea and Draupner rogue waves[12, 13], where second order bound nonlinearities are the dominant factor that enhances the linear dispersive focusing of extreme waves.

Figure 1

Metocean parameters history during January 2017 storm: (top panels) hourly variation of the significant wave height estimated as (black) and (red), mean zero-up-crossing wave period , and depth factor ; (bottom-left) normalized frequency spectrum and (bottom-right) the directional spectrum  () at the storm peak. Here,  () and  () are the spectral amplitude and frequency at the storm peak.

Figure 2

Metocean parameters history during Gaia storm: (top panels) hourly variation of the significant wave height estimated as (black) and (red), mean zero-up-crossing wave period , and depth factor ; (bottom-left) normalized frequency spectrum and (bottom-right) the directional spectrum  () at the storm peak. Here,  () and  () are the spectral amplitude and frequency at the storm peak.

Figure 3

Metocean parameters history during storms  Jan 2017 (left panel, start date Jan 25, 2017, 00:00) and Gaia (right panel, start date Jan 16, 2017, 00:00): hourly variation of Boccotti parameter and the significant wave height estimated as .

Statistical measures of ocean nonlinearities

We use integral statistics such as the Tayfun steepness and the excess kurtosis of the zero-mean surface elevation to measure ocean wave nonlinearities. Wave skewness describes the effects of second-order bound nonlinear effects, which alter the geometry and statistics of the sea surface with higher sharper crests and shallower more rounded troughs[23-25]. The excess kurtosis indicates deviations of the distribution of surface elevations from Gaussianity. It comprises a dynamic component[19, 26] measuring third-order quasi-resonant wave-wave interactions and a bound contribution induced by both second- and third-order bound nonlinearities[23-28].

Figure 4 depicts the statistical parameters history for the two analyzed storms. The top panels display the hourly variations of the observed Tayfun steepness , the observed excess kurtosis and the theoretical maximum Narrow Band (NB) dynamic excess kurtosis[19] for storm Jan 2017. The theoretical narrow-band predictions are also shown. The same comparisons for storm Gaia are also reported in the bottom panels. In Table 1 we compare the statistical parameters of Gaia and Jan 2017 strongest sea states and the Draupner and Andrea rogue sea states (from Fedele et al.[14]). All sea states are short-crested as characterized by large directional spreading, especially at the peak of the storms as . In open ocean, wave energy can spread directionally and third-order quasi-resonant nonlinearities are essentially insignificant[12, 19, 29–32] in such realistic oceanic conditions[12, 33, 34]. Indeed, the right panels of Fig. 4 indicate that the maximum dynamic excess kurtosis is of  and negative for both Jan 2017 and Gaia, and negligible in comparison to the associated bound component. The values of the observed excess kurtosis from the hourly 18 min samples are erratic due to the short record duration. Estimates of the third order moment, or wave skewness from the buoy’s time records of surface elevations converged for time series 70–100 wave periods long, or 12–15 min records. In contrast, estimates of the fourth-order moment, or excess kurtosis were at the verge of convergence for the available 18-min record. Nevertheless, our analysis indicates that third-order quasi-resonant interactions, including NLS-type modulational instabilities play an insignificant role in the formation of large waves[12, 19] in agreement with oceanic observations available so far[24, 35, 36].

Figure 4

Statistical parameters history during storms (top panels) Jan 2017 and (bottom panels) Gaia: (left) hourly variation of the theoretical NB (bold line) and observed Tayfun steepness  (thin blue line), (center) observed excess kurtosis (thin blue line) and NB estimate (bold line) and (right) theoretical maximum NB dynamic excess kurtosis[19] (blue thin line). The significant wave height history is also reported (black line).

Next we analyze the local steepness of the individual waves observed. To do so, we consider the two definitions and based on the crest height h and the half of the crest-to-trough height H, respectively. Here, k is the wavenumber estimated from the measured local wave period using the nonlinear dispersion relation[21, 37]. Toffoli et al.[21] used to analyze the statistics of individual wave steepnesses from both laboratory and open sea measurements. They concluded that waves can reach values above the Stokes limit 0.448 and up to the limit 0.55, hereafter referred to as TBOW. Clearly, H/2 is smaller than the crest height h on average because of the crest-trough asymmetry due to bound nonlinearities[23]. As a consequence is smaller than . Indeed, in a numerical study of breaking wave groups Barthelemy et al.[22] observed values that did not exceed the limit 0.72, which is based on the fifth order Stokes wave solution[38]. Recently, Knobler et al.[39] observed in their experiments waves of average amplitude with steepness above 0.8. The steepest waves appear as forms of wave dislocations due to crest pairing, or splitting[40]. However, the largest waves in their dataset did not exceed the theoretical limit[22] 0.72. Figure 5 depicts the local steepness of the individual waves measured during Jan 2017 and Gaia storms. Note that the observed values do not exceed the TBOW limit 0.55. Similarly, the observed values do not exceed the 0.72 limit proposed by Barthelemy et al.[22].

Figure 5

Observed Local steepnesses and versus local period for (top panels) Jan 2017 and (bottom panels) Gaia storms. Stokes limit = 0.448, TWBO limit[21] = 0.55 and Barthelemy et al.[22] = 0.72. k = local wave number, h = crest height, H = crest-to-trough height and  = mean zero-crossing period. Because of the crest-trough asymmetry, on average and .

Occurrence frequency of a rogue wave during storms

We now describe the statistics of large waves encountered by an observer at a fixed point of the ocean surface during a storm of duration . To do so, we draw on Fedele et al.[14] and model the storm as a non-stationary continuous sequence of sea states of duration dt, and is the number of waves in the sea state, where is the time-dependent mean zero-crossing wave period. We want to estimate the probability, or occurrence frequency of a wave crest height exceeding a given threshold  as encountered by a fixed observer, where . That is the probability to randomly pick from the non-stationary time series observed at a fixed point a wave crest that exceeds the threshold , and it is given by[14]where is the stationary exceeding probability of the sea state occurring in the time interval t and , which depends on wave parameters around time t. The definition of is consistent with the way wave crests are sampled from non-stationary wave measurements during storms. Moreover, the threshold  is exceeded on average once every waves, where is also the conditional return period of a wave whose crest height exceeds . In practice, the non-stationary is estimated as the weighted average[14]where is the number of waves sampled in the sea state j and a storm is partitioned in a finite sequence of sea states of duration . For weakly nonlinear random sea states (stationary seas), the probability P is hereafter described by the third order Tayfun-Fedele[24] (TF), modified narrow-band[41] (MNB), Tayfun[23, 25] (T), Forristall[43] (F) and the linear Rayleigh (R) distributions (see “Methods” section). Similarly, the non-stationary occurrence frequency of a wave in a storm whose crest-to-trough (wave) height exceeds the threshold are described by the same Eqs. (1) by simply replacing P(h) with the exceedance probability P(H) proper for wave heights of stationary seas. This is hereafter described by the generalized Boccotti[42] (B), Tayfun[44] (T) and linear Rayleigh (R) distributions (see “Methods” section).

Figure 6 summarizes the wave statistics for storm Jan 2017. In particular, the left panel depicts the empirical distribution (squares) for crest heights plotted versus the number of waves . This is compared against theoretical predictions of the non-stationary second-order F, MNB and T models as well as the third-order TF and Rayleigh (R) distributions. T and TF models are based on the observed Tayfun steepness . Note that TF is practically the same as MNB as an indication that second-order effects are dominant. F slightly overestimates the return periods, whereas the linear R model underestimates them. This indicates the dominance of second-order nonlinearities in shaping the sea surface. Similarly, the right panel of the same Fig. 6 reports the empirical distribution for crest-to-trough wave heights for Jan 2017. The observed statistics is well described by both the generalized Boccotti (B) and Tayfun (T) models. Similar conclusions also hold for the wave statistics for Gaia presented in Fig. 7.

Figure 6

Wave statistics for storm Jan 2017: (left) crest heights and (right) crest-to-trough wave heights versus number of waves . The empirical distribution (squares) from the observed wave ensemble is compared against with theoretical models: T = Tayfun[23, 25], TF = Tayfun-Fedele[24], MNB = modified narrow-band[41], B = generalized Boccotti[42], F = Forristall[43] and R = Rayleigh distributions. Confidence bands are also shown (light dashes). is the inverse of the exceedance probability , and similarly is for crest-to-trough heights. Horizontal lines denote the rogue threshold  and  for crest and wave heights[2].

Figure 7

Wave statistics for storm Gaia: (left) crest heights and (right) crest-to-trough wave heights versus number of waves . The empirical distribution (squares) from the observed wave ensemble is compared against with theoretical models: T = Tayfun[23, 25], TF = Tayfun-Fedele[24], MNB = modified narrow-band[41], B = generalized Boccotti[42], F = Forristall[43] and R = Rayleigh distributions. Confidence bands are also shown (light dashes). is the inverse of the exceedance probability for crest heights, and similarly is for crest-to-trough heights. Horizontal lines denote the rogue threshold  and  for crest and wave heights[2].

The wave profile with the largest wave crest height observed during Gaia ( m,  m) is shown in the right panel of Fig. 8. For comparison, the El Faro, Draupner and Andrea rogue wave profiles are also shown[12, 13]. From Table 1, Gaia wave did not exceed the rogue threshold  as the other three rogue waves did. Nevertheless, the five sea states have very similar characteristics, suggesting a common generation mechanism due to linear dispersive focusing enhanced by second-order nonlinearities[12].

Figure 8

Simulated wave profiles (solid) and mean sea levels (MSL) (dashed) versus for (from left to right) El Faro, Andrea and Draupner sea states. Gaia largest wave observed is shown in the right panel. Actual measurements (thin solid) and MSLs (thin solid) are also shown for Draupner. Note that the Gaia MSL is insignificant and the Andrea MSL is not available. is the maximum crest height and is the dominant wave period (see Table 1 and “Methods” section for definitions).

Wave statistics in space-time

The maximum surface wave height observed over a given area during a time interval, or space-time extreme is much larger than the maximum observed at a given point over the same time period[15, 16, 48]. Indeed, in the short-crested directional seas typically generated by storms it is very unlikely that an observed large crest at a given point in time is actually the largest crest of a group of waves propagating in space-time. In contrast, it is most likely that the sea surface was much higher somewhere near the measurement point, as Boccotti’s (2000) QD theory predicts.

Euler Characteristics of random fields[49-51] provide a theoretical framework for modelling space-time wave extremes stochastically[15, 16, 52]. Drawing on Fedele et al.[13] consider a 3-D non-Gaussian field in space-time over an area A for a time period of D. The area cannot be too large since the wave field may not be statistically homogeneous. The duration should be short so that the sea state can be assumed as stationary and spectral changes insignificant. Then, drawing on the linear Gaussian space-time model[49-51] the third-order nonlinear probability that the maximum surface elevation over the area A and during the time interval D exceeds the generic threshold is described by[13, 49–51]wheredenotes the third-order nonlinear Tayfun-Fedele model[24] for the probability of exceedance of crest heights at a fixed point. The parameter measures third-order nonlinear effects and can be estimated from the excess kurtosis as  (see “Methods” section). Note that our analysis indicates that third-order nonlinearities are negligible. The amplitude relates to the Rayleigh-distributed envelope via the Tayfun quadratic equation[11, 23, 25]

The space-time parameters and are the average number of 1-D and 2-D waves that can occur on the edges and boundaries of the volume , and is the average number of 3-D waves that can occur within the volume[15]. These all depend on the directional wave spectrum and its spectral moments defined in the “Methods” section.

Given the probability structure of the linear surface wave maximum in Eq. (3), the nonlinear mean maximum surface or crest height attained over the area A during a time interval D is given, according to Gumbel (1958), by[13]where is the inverse of the derivative of with respect , the Euler-Mascheroni constant , and the most probable surface elevation value follows from Eq. (3) as

Fedele et al.[13] observed that Gaussian space-time models[53-57] overestimate the maximum surface height over large areas and time intervals because they rely on Gaussianity. Indeed, Gaussian seas do not account for the wave physics of nonlinear dispersion[29], which limits the wave growth as a precursor to breaking[29, 58]. As a result, water surface elevations of a Gaussian sea can reach any height as there are no physical limits induced by wave breaking. Thus, the larger the area A or the time interval D, the greater the number of waves sampled in space-time, and unrealistically large amplitudes are likely to be sampled. Recently, Benetazzo et al.[59] proposed an approach to limit the effects of the unboundedness of rare Gaussian space-time events. In particular, they bound the range of wave crest amplitudes of a given sea state in the finite interval as waves inevitably break and limit crest heights. The threshold is chosen based on wave observations. For Draupner, Andrea and El Faro sea states Benetazzo et al.[59] set , but both Andrea and El Faro attained larger amplitudes (see Table 1). Such an approach was pioneered by Collins[45] and Battjes[46] in the early 70s (see Thornton and Guza[47] for a comprehensive review). Hereafter, we refer to it as the Collins-Battjes method. In simple words, it relies on defining a sea state of N waves whose sample space consists of non-breaking waves and breaking waves. Here, is the probability, or fraction of waves that do not break. All the mass probability of breaking waves is concentrated at and the associated probability density function (pdf) of crest heights has a Dirac delta spike of mass at so that the probability of the entire sample space is 1. The pdfs reported by Benetazzo et al.[59] in their Eqs. (13) and (16) miss the Dirac delta spike to concentrate the probability mass at . Following Collins-Battjes method[45-47], we bound the Tayfun-Fedele exceedance probability in Eq. (4) as

The correct pdf must have a Dirac delta spike at aswhere and the probability of the entire sample space adds up to 1. The space-time probability of exceedance can be estimated by the associated Gumbel model[15]limited in the range . The expected maximum surface height follows by integrating the exceedance probability over the finite range

The nonlinear mean maximum surface or crest height expected at a point during the time interval D follows from Eq. (3) by setting and . We note in passing that for Gaia and Jan 2017 storms the unbounded TF model well describes the statistics of large crests at a point (see Figs. 6 and 7), and so for El-Faro[13], Andrea and Draupner[12] rogue waves. The second-order form of the ST model above () has been implemented in WAVEWATCH III[60]. The linear limit follows from Eq. (3) by setting and .

To characterize rare space-time occurrences in third-order nonlinear random seas we consider the threshold exceeded with probability q by the maximum surface height over an area A during a sea state of duration D. This satisfies

In simple words, the maximum surface height observed within the area A during D exceeds the threshold only in realizations of an ensemble of N identical sea states. Consider the ST ratio , where is the mean maximum surface height expected over the area during a sea state of duration hour and is the mean maximum surface height expected at a point during the same time interval D. Figure 9 displays the unbounded ST ratios (dashed lines) as a function of the area width for the El Faro, Draupner and Andrea sea states[13] as well as the sea states at the peaks of Gaia and Jan 2017 storms. We also depict the bounded ratios (solid lines) using Collins-Battjes method and the threshold as in Benetazzo et al.[59]. The ST ratios for Draupner and Andrea are estimated using the European Reanalysis (ERA)-interim data[52]. For comparisons, we also report the empirical ST ratio from the HOS-simulated El Faro sea state[13] and the experimental observations at the Acqua Alta tower[16]. Fedele et al.[13] noted that the unbounded predictions fairly agree with the El Faro HOS simulations and observations for small areas (), but they overestimate over larger areas. Note that Gaia and Jan 2017 give similar predictions and slightly smaller than the other sea states. As expected, the bounded predictions limit the growth of the statistical ensemble of waves in space-time. In particular, the maximum surface height expected over an area is bounded by the threshold , and it does not exceed 1.6 times the surface height expected at a point.

Figure 9

Space-time extremes: theoretical ratio as a function of the area width for El Faro, Draupner and Andrea rogue sea states as well as the sea states at the peaks of Gaia and Jan 2017 storms. is the mean maximum surface height expected over the area during a sea state of duration hours and is the mean maximum surface height expected at a point. For comparisons, the empirical ST ratio from the El Faro HOS simulations (green solid line) together with the experimental observations at the Acqua Alta tower (squares) are also shown[16]. is the mean wavelength. Dashed lines denote the unbounded Gaussian predictions, whereas solid lines denote the bounded predictions using Collins-Battjes method[45–47] with .

Prediction of potential hazards for vessel navigation

Consider a vessel that navigates through a wave field of surface elevations  at a constant speed V along a straight path at an angle with respect to the x axis of a fixed Cartesian frame (x, y, t) over a time interval D. Define the occurrence frequency, or exceedance probability that a rogue wave exceeding a crest height h is encountered by the vessel during its navigation path. can be estimated as the probability that the maximum surface height exceeds h over the space-time volume spanned by the vessel in motion. To do so, we draw on Fedele et al.[13] and consider the encountered wave field  as seen by an observer in the frame  moving with the vessel, whereand (x, y, t) is the fixed Cartesian frame. The wave field encountered by the cruising vessel, or equivalently the field of surface fluctuations measured from the moving frame is

For an observer in the fixed frame (x, y, t), is a random wave field homogeneous in space and stationary in time, and so is the encountered field as seen from an observer on the vessel in the moving frame . However, will be non-stationary for a generic navigation route and variable vessel’s speed. The associated space-time covariance is given bywhere is the directional wave spectrum of the sea state, and , where k is the wavenumber related to the frequency f via the linear dispersion relation. The Doppler-shifted encountered, or apparent frequency is[61-63]

When the vessel moves faster than the waves coming from a direction , the apparent frequency and waves appear to move away from an observer on the vessel. In this case, the direction of those waves should be reversed[61], i.e. , and set as positive. The spectral moments of the encountered random field are given by[13]

In the moving frame , the vessel is stationary. As a consequence, the vessel spans the same area equal to its size L x W m over the duration D, i.e. it covers the space-time volume  m x D min. Here, L and W are the length and width of the vessel, respectively. The probability that the maximum surface height over exceeds the threshold h can be estimated from the spectral moments of the encountered wave field using the space-time model[13, 15, 52] in Eq. (3). In our analysis, we consider the worst-case scenario of the vessel cruising against the dominant wave direction of the sea state. In such conditions the moving vessel encounters more waves than if it were anchored at a fixed point because the apparent frequency of the encountered waves increases.

As an application, we hypothesize the scenarios of two vessels of the Israeli Navy fleet navigating at their cruise speeds against the dominant waves of the peak sea state of Gaia over a time interval of 10 min. We consider the Super Dvora Mk III-class patrol boat, a vessel of small size (length 20 m and width 6 m) cruising at the average speed of  50 km/h (27 knots) and maximum speed of 93 km/h (50 knots). We also consider the Sa’ar 6-class corvette, a large vessel (length 90 m and width 14 m) cruising at the average speed of 30 km/h (17 knots) and maximum speed of 50 km/h (27 knots). Figure 10 depicts the bounded predictions (solid lines) of the nonlinear threshold exceeded with probability 1/n by the maximum surface height seen by the observer on the two vessels over the duration  min. We used the Collins-Battjes method[45-47] described in the preceding section and adopt , which is the crest height of the largest simulated wave of the El Faro sea state[13]. The unbounded Gaussian predictions are also reported (dashed lines). Consider an observer on the large vessel navigating against the dominant waves of the peak sea state of Gaia at the average cruise speed of 30 km/h. The observer sees the vessel as stationary and the spanned space-time volume is  x 14  x10 min. The probability that the maximum surface height exceeds  m (here,  m, see Table 1) when the vessel cruises along the dominant direction of waves is close to , as shown in Fig. 10. An observer on the small vessel navigating against the dominant waves at the average cruise speed of 50 km/h has a smaller probability to encounter the rogue wave, because the vessel’s size is smaller than that of the large vessel. If the two vessels navigate along with the dominant waves reduces to 1/2400 (large vessel) and 1/10,500 (small vessel), because they will encounter less waves as the apparent wave frequency reduces.

Figure 10

Theoretical bounded predictions (Collins-Battjes method[45–47], ) for the maximum surface height exceeded with probability 1/n over the space-time volume  m x D min spanned by the El Faro vessel (black lines), the Super Dvora Mk III-class patrol boat (red lines, small vessel) and the Sa’ar 6-class corvette (blue lines, large vessel). The vessels navigate against the dominant waves in the peak sea state of storm Gaia for a duration D. L and W are the length and width of the vessel. The numerical HOS predictions (squares) for El Faro are also reported[13]. Dashed lines denote the unbounded Gaussian predictions. Dashed horizontal line denotes the threshold.

For an observer fixed at a point of the sea, reduces to  (see Fig. 8). We note in passing that if the two vessels navigate at their maximum speeds against the dominant waves, increases to 1/3160 (small vessel) and 1/1000 (large vessel).

Figure 10 also depicts the predictions for the El Faro vessel during Hurricane Joaquin based on HOS simulations[13]. The data suggests that the El Faro vessel was drifting at an average speed of approximately 2.5 m/s prior to its sinking[13]. Over a time interval of min prior sinking the probability to encounter a rogue wave that exceeds the threshold  m is , in agreement with the predictions by Fedele et al.[13]. Thus, we find the same probabilities for the hypothetical scenario of the El Faro vessel drifting through the sea state at the storm peak of Gaia. However, in this case the rogue wave would exceed a smaller height, i.e.  m since the significant wave height at the peak of Gaia ( m) is smaller than that of the El Faro sea state ( m, see Table 1).

Conclusions

The Eastern Mediterranean is a basin of increasing interest being a climatology hot-spot and a natural treasures area, both ecologically and economically[1]. However, there is a lack of studies on the metocean characteristics of its extreme sea states and large waves. Such studies are of paramount importance, especially in light of global warming where weather events are likely to become more severe and more hazardous to offshore and coastal structures as well as to vessel navigation.

To fill the gap of knowledge on metocean properties of the Eastern Mediterrenean sea storms, we have presented a scientific investigation and prediction of wave extremes of that area. Our statistical analysis indicates that the largest observed waves during the two major storms recorded in the last five years have characteristics quite similar to those displayed by the El Faro, Andrea, and Draupner rogue waves[12, 13]. In particular, linear dispersive focusing enhanced by second order bound nonlinearities are the dominant factor in the generation of large waves. Third-order resonances and NLS-type modulational instabilities are ineffective in shaping large waves[12, 19]. Finally, we have carried out a novel space-time analysis of the most intense sea states and demonstrated the ability to predict potential rogue hazards for ships of various sizes and cruise speed navigating within the rough waters of the analyzed storms.

Methods

Buoy measurements

The ‘3 metre’ buoy manufactured by AXYS Technologies is located on the south part of the Exclusive Economic Zone of Israel. The buoy is moored at a depth of 260 m at the coordinates ( Lat., Long.). This buoy operates by recording 18 minuets of heave displacements at the beginning of each hour with sampling frequency of  Hz. Additional outputs are also elaborated by the internal buoy’s software, i.e. the directional energy density spectrum spectrum with a resolution of 3° covering 360°, and typical wave parameters such as the significant wave height, wave period and steepness. The same wave parameters were estimated independently from the the time series of the surface wave fluctuations provided by the buoy, and they were practically identical to the buoy’s outputs. For a cross-site validation, the buoy’s wave parameters have been compared to data acquired by the nearby onshore (i) Ashkelon port meteorological station[64] and (ii) the DATAWELL Wave-Rider MKIII buoy[65] 40-km away from our buoy, and moored 2 km offshore at a 24-m depth. In both cases, the buoy’s wave parameters showed similar trends to those acquired by the two instruments, while differing slightly in values. This is expected in light of the distance between the instruments’ locations and the physical characteristics of shallow waters and onshore environments.

Wave parameters

The significant wave height is defined as the mean value of the highest one-third of wave heights. It can be estimated either from a zero-crossing analysis or more easily but approximately from the wave omnidirectional spectrum as , where is the standard deviation of surface elevations, are spectral moments. Further, is the directional wave spectrum with as the direction of waves at frequency f, and the cyclic frequency is . In this paper, we use the spectral-based estimate, which according to our wave measurements is larger than estimated from the measured time series.

The dominant wave period refers to the cyclic frequency of the spectral peak. The mean zero-crossing wave period is equal to , with . The associated wavelength follows from the linear dispersion relation , with d the water depth. The mean spectral frequency is defined as [23] and the associated mean period is equal to . A characteristic wave steepness is defined as , where is the wavenumber corresponding to the mean spectral frequency [23]. The following quantitites are also introduced: , the phase velocity , the group velocity .

The spectral bandwidth gives a measure of the frequency spreading. The angular spreading , where and is the mean direction. Note that . An alternative measure of spectral bandwidth is given by the Boccotti parameter , which is the absolute value of the first minimum of the normalized covariance function of the zero-mean surface displacement , which is attained at [66] and the corresponding second derivative with respect to [66].

The wave skewness and the excess kurtosis of the zero-mean surface elevation are given by

Here, overbars imply statistical averages and is the standard deviation of surface wave elevations. Note that was first derived by[67]. For second-order waves in deep water[25]and[18]

Here, is the spectral bandwidth defined above and the characteristic wave steepness , where is the wavenumber corresponding to the mean spectral frequency [23]. For narrow-band waves, tends to zero and the associated skewness [23-25]. The excess kurtosis of weakly nonlinear random seascomprises a dynamic component due to nonlinear quasi-resonant wave-wave interactions[17, 26] and a Stokes bound harmonic contribution [20]. In deep water it reduces to the simple form [20, 26, 68] where is the skewness of narrowband waves[23]. In deep water, the dynamic component is given in terms of a six-fold integral[17, 19] that depends on the Benjamin-Feir index and the parameter , which is a dimensionless measure of the multidirectionality of dominant waves[19, 26, 69].

The Tayfun-Fedele model for crest heights

The probability that a wave crest observed at a fixed point of the ocean in time exceeds the threshold can be described by the third-order Tayfun-Fedele model[24],where follows from the quadratic equation [23]. Here, the Tayfun wave steepness  is of and it is a measure of second-order bound nonlinearities as it relates to the skewness of surface elevations[25]. The parameter  is a measure of third-order nonlinearities and is a function of the fourth order cumulants of the water surface elevation and its Hilbert transform [24]. In particular,  and . In our studies is approximated solely in terms of the excess kurtosis as  by assuming the relations between cumulants[70] and . These, to date, have been proven to hold for linear and second-order narrowband waves only[27]. For third-order nonlinear seas, our numerical studies indicate that within a relative error in agreement with observations[71, 72].

For second-order seas, referred to as Tayfun sea states[73], only and in Eq. (22) yields the Tayfun (T) distribution[23]

For Gaussian seas, and and reduces to the Rayleigh (R) distribution

Fedele and Tayfun[25] proved that the Tayfun distribution represents an exact result for large second order wave crest heights and it depends solely on the steepness parameter defined as .

The modified narrowband (MNB) model

Tayfun and Alkhalidi[41] introduced the MNB model as an improvement of the conventional weakly-nonlinear second-order finite-depth narrowband model so as to avoid certain unrealistic features of that model in relatively shallow waters and to extend its applicability to highly nonlinear waves propagating from deep water to the shoaling and surf zones. The MNB exceedance probability for crest heights is given bywhere follows from the quadratic equation[11, 41]whereand the associated excess kurtosis

The MNB model is valid over the range of skewness values .

The Forristall model

The exceedance probability is given by[43]where , for multi-directional (short-crested) seas. Here, is a characteristic wave steepness and the Ursell number , where is the wavenumber associated with the mean period and d is the water depth.

The generalized Boccotti and Tayfun models for crest-to-trough (wave) heights

The third-order nonlinear statistics for crest-to-trough wave heights is described in terms of the generalized Boccotti distribution[42, 66].where the Boccotti parameter  and are defined above in the section where statistical parameters are described. For Gaussian seas (), the original Boccotti model[66] is recovered

The Tayfun model for wave heights is given by[24, 44]where is the value of the envelope r(t) of the covariance at .

Space-time statistical parameters

For space-time extremes, the coefficients in Eq. (3) are given by[15, 48]whereare the average number of waves occurring during the time interval D and along the x and y sides of length and respectively. They all depend on the mean period , mean wavelengths and in x and y directions:and

Here,are the moments of the directional spectrum and