Large nearshore storm waves off the Irish coast.

We present a statistical analysis of nearshore waves observed during two major North-East Atlantic storms in 2015 and 2017. Surface elevations were measured with a 5-beam acoustic Doppler current profiler (ADCP) at relatively shallow waters off the west coast of Ireland. To compensate for the significant variability of both sea states in time, we consider a novel approach for analyzing the non-stationary surface-elevation series and compare the distributions of crest and wave heights observed with theoretical predictions based on the Forristall, Tayfun and Boccotti models. In particular, the latter two models have been largely applied to and validated for deep-water waves. We show here that they also describe well the characteristics of waves observed in relatively shallow waters. The largest nearshore waves observed during the two storms do not exceed the rogue thresholds as the Draupner, Andrea, Killard or El Faro rogue waves do in intermediate or deep-water depths. Nevertheless, our analysis reveals that modulational instabilities are ineffective, third-order resonances negligible and the largest waves observed here have characteristics quite similar to those displayed by rogue waves for which second order bound nonlinearities are the principal factor that enhances the linear dispersive focusing of extreme waves.

Introduction

Recent studies[1,2] reveal that rogue waves can arise from a combination of the process of constructive interference and nonlinear effects specific to the complex dynamics of ocean waves. Under relatively rare conditions, waves locally propagate in an organized way or nearly in phase, resulting in an unusual case of constructive interference that generates waves with large amplitudes. However, this mechanism still cannot fully explain the sizes of rogue waves observed under actual oceanic conditions. Various discrepancies observed between theoretical models and actual observations can be attributed to the nonlinear nature of waves: they are not sinusoidal but vertically asymmetric, displaying shallower more rounded troughs, and higher sharper crests that result from the water surface being pushed upward against the pull of gravity. Thus, the nonlinearity of the ocean surface manifest in the lack of symmetry between wave crests and troughs needs to be accounted for[3-5]. Such nonlinearities do contribute to the effects of constructive interference noticeably. Indeed, recent studies[1,4] suggest that nonlinear effects due to second-order bound harmonics play a predominate role in this process and can cause an increase of 15 to 20 percent in crest height, i.e. the vertical distance from the mean sea level to the top of the wave.

The formation of a rogue wave at a given point of the ocean is simply a random or chance event[1]. Several cases of extreme wave occurrences of practical and theoretical interest such as the Andrea, Draupner and Killard waves[1] and the sinking of El Faro[2] have been studied in detail by way of higher order spectral wave simulations and validated with probabilistic wave models. These studies have shown that second-order statistical distributions of crests, in particular those often referred to as Tayfun[3,4,6,7] and Forristall[8] models, both describe rogue statistics reasonably well in intermediate to deep waters.

In this work, we will show that the Tayfun and Boccotti[9,10] models for wave heights, previously validated for both simple[4,10] and mixed seas[11-13] in deep water, describe the statistics of large waves in intermediate to relatively shallow waters reasonably well also. For comparison, we also consider the Forristall’s Weibull regression model[8] because of its frequent application and popularity in engineering design[14]. Our results here and several others elsewhere indicate that it does work quite well in describing the distribution of observed data. For example, Gibson et al.[14] use the Forristall’s model to explain their statistics of oceanic wave crests, but they do not use the Tayfun[3] or Tayfun-Fedele[4] models on the grounds that they require the calculation of a key parameter from the time trace of water surface elevations, not readily available from hindcast models. However, they were able to use the Boccotti[9] wave-height model with two parameters specifically dependent on the frequency spectrum or the time trace of surface elevations. The Forristall crest-height model likewise requires two parameters that depend on the spectral moments. Similarly, recent work by Katsardi et al.[15] indicates that various regression models fitted to observed data are not universally applicable nor do they provide an adequate description of waves in shallow water, as large waves are overestimated. However, only unidirectional laboratory waves propagating on rather mild impermeable slopes are explored, and neither the Tayfun crest-height model[3,4] nor the Boccotti wave-height models[9,10] are tested on the grounds that second-order models cannot describe highly nonlinear shallow-water waves affected by intense wave-breaking. Nonetheless, they consider the Forristall crest-height model in their comparisons. That model is also second-order, and all second-order models break down in shallow water where waves are highly nonlinear, prone to intense breaking and affected by various dissipative effects of the seabed. Nevertheless, all these models are equally applicable to some of their data representative of the relatively shallower waters of the transitional water depths.

In the present study, we consider all the aforementioned models and test them against directional waves observed in relatively shallow waters within the shoaling zone. In particular, we consider two wave data sets from ADCP measurements taken off the west-coast of Ireland near Killard Point in 2015 and near the Aran Islands in 2017 (see Fig. 1). In particular, the two locations are nearshore at a water depth of approximately 37 meters (Killard) and 45 meters (Aran). They are well-known high-energy coastlines where storm waves overtop cliffs, fracture bedrock, and move large rocks weighing 100 tons or more[16-19]. We then analyze wave statistics in relatively shallow waters during storm events. In particular, we examine data observed in two storms, namely the storm of 25–27 Feb 2015, hereafter referred to as Feb 2015, and Doris of 21–26 Feb 2017. Wave measurements carried out during these storms are described in the Methods section. Monochromatic waves propagating on a water depth d and characterized with a wave number k feel the presence of the bottom, and start being modified whenever the dimensionless depth parameter kd < π. The wave regime is classified as deep water if kd > π, as intermediate or transitional depth for π/10 < kd < π, and as shallow water if kd < π/10. This classification has significance both practically and theoretically in establishing how wave characteristics are modified and what processes need be included and modeled in their theoretical predictions. Obviously, defining the depth regime of a wind-wave field as a whole in a similarly precise fashion is impossible due to the wide range of wave numbers observed. As a compromise, we will define the depth regime for the two storms based on the dominant wavenumber k at the spectral peak, thus focusing our attention on the most energetic components with wave numbers at and near the spectral peak. On this basis, both storms are in the transitional water-depth regimes since 0.5 < kd < 2.5, as seen in the right panel of Fig. (2). Further, characteristics of both storms vary considerably over their durations of 70 hours, approximately. As a consequence, we propose novel probability models appropriate to non-stationary processes so as to be able to analyze the surface-elevation time series gathered during the two storms.

Figure 1

(Left) Map of Ireland: ADCP location off the Aran Islands in 2017 (upper Inset) and off Killard Point in 2015 (lower Inset). (Right) ADCP deployment with (blue-capped) instruments in a protective frame. The map was generated from data via OpenStreetMap and its contributors. Photo by F. Dias.

Figure 2

(Left) Variations of difference between consecutive sea states of Doris and Feb 2015 storms, measured in terms of the standard deviation of observed values of σ2/σ1 − 1 between two consecutive sea states in a storm sequence. (Right) Hourly variations of the depth factor kd and (dashed line) theoretical threshold above which plane waves are modulationally unstable in unidirectional seas[30,31].

Results

This section is structured as follows. First, we discuss the characteristics of the sea states generated by the two storms as they pass by the west coast of Ireland. The descriptions of various wavefield characteristics, associated principal statistical parameters and probability models employed in the analyses are described in the Methods section. Subsequently, we present the analysis of extreme waves. In particular, in order to be able to predict rationally the occurrence of extreme waves in a time-dependent storm, we first present the theoretical formulation of a non-stationary model for describing the sequences of sea states in a storm. On this basis, an optimal sea state duration is determined based on the rationale that variation of key statistics is minimal between two consecutive sea-state sequences. We then explore the occurrence frequency of rogue waves observed at a fixed point at sea. The largest waves observed at the peak of the storms and their characteristics are then compared to those of the Draupner and Andrea rogue waves, observed at different oil platforms in the North Sea in 1995 and 2007, respectively, the Killard rogue wave observed off the coast of Ireland in 2015[1] and the simulated El Faro rogue wave[2]. The metocean parameters of the six sea states are summarized in Table 1.

Table 1

Wave parameters and various statistics of Doris and Feb 2015 at the storm peak in comparison to the El Faro[2], Andrea, Draupner and Killard rogue sea states[1].

	Doris	Feb 2015	El Faro	Andrea	Draupner	Killard
Significant wave height Hs [m]	6.4	12.6	9.0	10.0	11.2	11.4
Dominant wave period Tp [s]	10.0	14.1	10.2	14.3	15.0	17.2
Mean zero-crossing wave period T0 [s]	8.8	13.3	9.2	11.6	12.1	14.0
Mean wavelength L0 [m]	120	216	131	209	219	268
Depth d [m], k0d with k0 = 2π/L0	45, 2.35	37, 1.08	4700, 225	74, 2.23	70, 2.01	39, 0.91
Spectral bandwidth ν	0.46	0.51	0.49	0.35	0.36	0.37
Angular spreading σθ [rad]	0.90	1˜.07	0.79	0.37	0.39	0.34
Parameter R = σθ2/2ν2 [37]	1.90	2˜.2	1.34	0.56	0.59	0.42
Benjamin Feir Index BFI[38]	0.18	0.22	0.36	0.24	0.23	0.18
Narrow-band (NB) skewness λ3,NB[49]	0.190	0.221	0.262	0.159	0.165	0.145
Observed skewness λ3	0.144	0.441	0.162	0.141	0.146	0.142
Maximum NB dynamic excess kurtosis λ40,maxd [39]	−10−3	−10−1	10−3	2.3 · 10−3	2.1 · 10−3	2.7 · 10−4
NB bound excess kurtosis λ40,NBd [41]	0.094	0.229	0.049	0.065	0.074	0.076
Observed excess kurtosis λ40	0.098	0.263	0.042	0.041	0.032	−0.011
Actual maximum crest height h/Hs	1.11	1.23	1.68	1.63	1.55	1.44
Actual maximum crest-to-trough (wave) height H/Hs	2.06	1.93	2.6	2.30	2.10	2.00

Note that the Killard rogue wave occurred on a water depth of 39 m, however the hincast input spectrum used by Fedele et al.[1] could only be computed at an averaged water depth of 58 m. We refer to the Methods section for the definitions of wave parameters.

Our statistical analysis of large waves focused on the study of the time sequence of changing sea states during the two storms. An optimal sea state duration T of 50 minutes for both storms was determined so as to minimize the degree of difference between waves of consecutive sea states. Drawing on Boccotti[9], this is measured by the standard deviation of the random variable V = σ2/σ1 − 1, where σ2 are the variances of two successive sea states in the storm sequence as shown in Fig. (2). Sampled values of V are obtained by dividing the non-stationary time series into N = D/T successive sea states of the same duration T and variances σ12,σ22, ..., σ2,σ + 12, ..., where D is the storm duration (70 hours). Then, N − 1 sampled values of V follow as V = σ + 1/σ − 1 for j = 1, ..., N − 1, from which the standard deviation of V can be estimated. Obviously, the mean of V tends to zero as T approaches smaller values.

This process ensured that resulting statistics are robust to variations in T up to ±20 min once the total population of surface elevations from each sea state in a storm sequence are normalized by the respective significant wave height.

Metocean parameters

Both storms generated directional sea states in transitional water depths. This is clearly seen in the right panel of Fig. (2), displaying the hourly variations of the depth coefficient kd during the two storms. Surface spectra observed were broadbanded so that ω ≈ 0.8ω, while k ≈ 0.7k.

Metocean parameters of both storms and how these vary are shown below in the top and bottom panels of Fig. (3) respectively. In particular, the left panels of Fig. (3) depict the hourly variation of the significant wave height H = 4σ. For comparison, the variation of actual significant height H1/3 representing the mean of the highest 1/3 of wave heights observed is also shown in the same panels. It is seen that it is about 5% smaller than H. The actual mean zero-up-crossing wave period T0 is also shown in the center panels of the same figure, whereas the right panels depict the wave spectra measured at the storm peaks. We observe that the high-frequency behavior in both cases is described by a logarithmic f−4 decay in conformity with Zakharov’s wave turbulence[20].

Figure 3

History of metocean parameters in (top panels) Doris and (bottom panels) Feb 2015 storms: (left) hourly variations of significant wave height H = 4σ compared to actually observed H1/3, (center) variation of mean zero-up-crossing period T0, and (right) surface spectra at peak stage of storms and their high-frequency saturation compared to a logarithmic f−4 type decay.

The two states analyzed here do not present any characteristics typical of mixed or crossing seas such as swell waves overlapping with wind seas because the frequency spectra S(f) displays a unimodal structure, as depicted in the right panels of Fig. (3). In particular, an examination of the directional spectrum S(ω, θ)/σ2 estimated using the Bayesian direct method (BDM) at the peak of Doris storm and shown in the left panel of Fig. (4) clearly displays a unimodal broad-banded wind-wave field, also confirmed by the attendant unimodal directional spreading function D(θ) in the right panel of the same figure.

Figure 4

Doris: (left panel) Estimated normalized directional spectrum S(ω,θ)/σ2 [Hz−1rad−1] at the storm peak using the Bayesian direct method (BDM) and (right panel) associated angular spreading function .

Figure (5) displays the scatter diagrams of crest heights h/H versus corresponding wave periods T/T0 for both storms. Large crest heights (and similarly wave heights, not reported here) do not violate the Miche-Stokes limits. These are depicted in the same figure by two bold red lines representing the Miche-Stokes limits for the most intense and weakest sea states of the storms. In seas generated by intense storms, nonlinear wave dispersion is effective in limiting wave growth as a precursor to breaking[21-23]. Thus, the onset of wave-breaking can occur well below the Miche-Stokes upper limit[22,24-27].

Figure 5

Scatter diagrams of crest heights h/H versus wave periods T/T0 observed in Doris and Feb 2015 storms compared with Miche-Stokes limits for the weakest (upper red curves) and the most intense (lower red curves) sea states, respectively.

In Table 1 we compare the metocean parameters observed during the peak states of the two storms with those of the El Faro, Draupner, Andrea and Killard rogue sea states[1]. Clearly, all six sea states have similar metocean characteristics. Killard, Doris and Feb 2015 are in shallower waters and the last two have a greater steepness than the other four sea states. Indeed, the observed values of skewness λ3 and excess kurtosis λ40 are larger than those observed in the other four cases (see also Fig. (6)). This suggests that the largest waves observed were near the onset of incipient breaking or already breaking, thus lessening the likelihood of occurrence of larger rogue events[21,22,28].

Figure 6

Histories of statistical parameters in (top panels) Doris and (bottom panels) Feb 2015 storms: (left) hourly variation of observed values of μ = λ3/3[6] compared to theoretical NB steepnesses[3,4] μ and μ; (right) hourly variation of observed Λ and Λ = 8λ40/3 implied by NB theory[41] compared to NB estimates Λ and Λ. Note that Λ is practically the same as Λappr. Subscripts m and p refer to definitions of parameters based on the mean and dominant wavenumbers k and k, respectively.

Modulational instability in intermediate or transitional water depths

In the Feb 2015 storm, the depth coefficient kd, depicted in the right panel of Fig. (2), was below the critical depth threshold 1.363 whereas not so in Doris. Above the threshold, plane waves are modulationally unstable[29] in the one-dimensional (1-D) wave dynamics described by the Nonlinear Schrödinger (NLS) equation[30,31]. However, the wave fields analyzed here are directional sea states, and according to the 2-D hyperbolic NLS equation, plane waves are modulationally unstable even at depths below that critical value, if they are perturbed by appropriate oblique disturbances[29,32-34]. Nevertheless, it is also recognized that instabilities ensuing from such disturbances are not likely to occur for values of kd < 0.5[32]. So, it is plausible that rogue waves could be generated by modulational instability, as in unidirectional seas[35,36] during both of the storms analyzed here. The kurtosis statistics is often used as an indicator if any rogue waves are present in a sea state. In sea states where third-order nonlinearities are significant, excess kurtosis λ40 = λ40 + λ40 comprises a dynamic component[37-40] λ40 due to nonlinear quasi-resonant wave-wave interactions and a Stokes bound harmonic contribution[40,41] λ40, given in the Methods section. As for the dynamic component, drawing on Janssen[40], Fedele’s[39] one-fold integral formulation is extended to narrowband (NB) waves in intermediate waters as

In the preceding expression, defines the Benjamin-Feir index in deep water at the spectral peak, ν the spectral bandwidth, and Im(x) denotes the imaginary part of x, ω and k the dominant spectral frequency and wavenumber. Depth effects on wave directionality, measured by R, are represented by R = βR by way of the factor β, and α is the depth factor. The latter two depend on the dimensionless depth kd (see Methods section). In the deep-water limit, both α and β become 1. The maximum of dynamic excess kurtosis is well approximated by[39]andwhere (which corrects a misprint in[39]) and b = 2.48.

In the focusing regime (0 < R < 1 and α > 0) the dynamic excess kurtosis of an initially homogeneous Gaussian wave field grows, reaches a maximum at the intrinsic time scale and then monotonously decreases and eventually vanishes over longer times. Such a regime is typical of unidirectional narrowband waves on water depths kd < 1.363. Indeed, in 1-D waves modulational instability disappears below that critical threshold and α < 0. As a result, wave dynamics becomes of defocusing type and excess kurtosis is negative. In particular, it reaches a minimum at t and then tends to zero over larger times. The kurtosis formulation in Eq. (1) extends Fedele’s[39] stochastic approach to NB waves at intermediate water depths, and it indicates that in directional seas such as those analyzed in this study, modulation instability can also occur at depths below the critical threshold 1.363 for α > 0 (see also Alber[42]). Further, for waves propagating over a broad range of directions in the open sea, Fedele et al.[1] show that such instabilities are ineffective in triggering rogue waves as excess kurtosis becomes negative, provided that R > 1. A rogue wave regime is a more likely occurrence only if the surface spectrum is sufficiently narrow-banded (R < 1) as well as characterized by a relatively large positive excess kurtosis.

Both storms analyzed here are in transitional water depths and prone to potential rogue occurrences induced by modulational instability since α > 0. However, all their sea states are directionally spread and characterized with mostly negative excess kurtosis since R > 1. This can clearly be seen in Fig. (7). Thus, the potential recurrence or focusing of large waves as observed in unidirectional seas is largely attenuated or suppressed[1,2,7,39]. Indeed, our statistical analysis indicates that the effects of third-order resonance or modulational instabilities are negligible, and that second-order bound nonlinearities are the dominant factor in shaping the large waves observed. We also point out that the NB predictions based on the mean wavenumber k yield similar negligible estimates of the dynamic excess kurtosis.

Figure 7

Variations of hourly estimates in (top panel) Doris and (bottom panel) Feb 2015 storms of the maximum of dynamic excess kurtosis λ40 (thin line) and significant wave height H (thick line).

In summary, our analysis indicates that fourth order cumulants are essentially trivial to begin with, implying that the analyzed sea states are ordinary: nothing specially rogue about them. The present analyses of the two storm wave datasets discussed in the following section confirm all this and show that there is hardly anything beyond second-order nonlinearities to explain their statistics.

Nonlinear wave statistics

The relative importance of nonlinearities in a sea state can be assessed by way of various integral statistics. These include the observed values of the wave steepness μ = λ3/3 [6] and the coefficient Λ = λ40 + 2λ22 + λ04 [4], where λ3 is skewness and λ are fourth-order cumulants of the zero-mean surface elevation η(t) and its Hilbert transform. In particular, λ40 is the excess kurtosis of surface elevations. Skewness is a measure of vertical asymmetry, and it describes the effects of second-order bound nonlinearities on the geometry and statistics of the sea surface with higher sharper crests and shallower more rounded troughs[3,4,6]. Excess kurtosis indicates whether the tail of the distribution of surface elevations is heavier or lighter relative to a Gaussian distribution. It comprises a dynamic component λ40 measuring third-order quasi-resonant wave-wave interactions and a bound contribution λ40 induced by both second- and third-order bound nonlinearities[3-6,37,43].

In describing wave statistics, the theoretical NB predictions based on the mean wavenumber k, rather than k, tend to yield more favorable comparisons with deep-water observations or theories[6,44]. Definitions based on k lead to predictions that noticeably underestimate the observed and/or theoretical statistics of broadband waves since k < k[44]. Nonetheless, the sea states analyzed here are in intermediate water depths and characterized with broad-banded spectra both in frequency and direction. Describing the statistics in such cases based on either k or k is neither very reliable nor realistic.

In Table 1 we compare the statistical parameters of the most intense sea states of Doris and Feb 2015, and also the Draupner, Andrea and Killard rogue sea states[1,2]. The maximum dynamic excess kurtosis is negative and negligible. Thus, third-order quasi-resonant interactions, including NLS-type modulational instabilities, should not play any significant role in the formation of large waves in comparison to bound nonlinearities[1,39] especially as the wave spectrum broadens[21] in agreement with oceanic observations available so far[4,45,46]. The values of excess kurtosis λ40 and Λ are mostly due to bound nonlinearities[7,47,48].

In the top panels of Fig. (6), we compare the hourly variations of (left) the observed values of μ = λ3/3 [6] versus the theoretical NB approximations[3,49] μ and μ based on the mean and dominant wavenumbers k and k, respectively, and (right) the observed fourth–order Λ coefficient, its approximation Λappr and the NB predictions[37,41] Λ and Λ based on k and k for Doris. The same comparisons are also reported in the bottom panels for Feb 2015. Clearly, the two NB predictions μ overestimate and μ slightly underestimate the observed values of μ for Doris. In contrast, both NB predictions underestimate the observed steepness for Feb 2015. Similarly, the actual Λ is mostly overestimated by both the NB estimates. Moreover, Λ is practically the same as Λappr (see Methods section). In the final analysis, the sea states analyzed here are characterized by broadband spectra both in frequency and direction and the NB assumption is unrealistic. Thus, hereafter we use the observed values of μ and Λ to evaluate wave statistics.

Occurrence frequency of extreme waves in storms

We now describe a novel approach for the statistics of extreme waves encountered by an observer at a fixed point of the ocean surface during a storm of duration D. Drawing on[9,50-52], the storm is modeled as a non-stationary continuous sequence of sea states of duration dt, and dt/T0(t) is the number of waves in the sea state, where T0(t) is the time-changing mean zero-crossing wave period. Consider now a wave of the storm and define the probability P(ξ) that its crest height C exceeds the threshold h = ξH as observed at a fixed location where H = 4σ. Equivalently, this is the probability of randomly picking a wave crest whose height C exceeds the threshold ξ = h/H from the non-stationary time series observed at a fixed point of the storm. Then,where P(ξ, t) = Pr[C > ξH(t)] stands for the probability that a wave crest height C exceeds the threshold ξH(t) in the sea state occurring in the time interval [t, t + dt]. This probability depends on wave parameters around time t. The definition of P is consistent with the way wave crests are sampled from non-stationary wave measurements during storms. In particular, a storm is partitioned into a finite sequence of N sea states of duration T = D/N. In each sea state S centered at t = t, waves are sampled and their crest elevations (h) are normalized with the local significant wave height as h/H(t) and put all in the same population . Then, the empirical exceedance probability P(ξ) is estimated as the ratio of the number of waves that exceed ξ to the total number of waves in the population. This is consistent with the way Eq. (4) is formulated. Indeed, is the expected number of waves during a sea state in [t, t + dt] and P(ξ, t)EX(t)dt is the number of waves whose crest heights exceed the threshold h = ξH(t) in the same sea state. Then, P in Eq. (4) follows by cumulating (integrating over time) the number of waves of all the sea states whose crest heights exceed h. In practice, the non-stationary P is estimated from data as the weighted averagewhere N(t) is the number of waves sampled in the sea state S. Equation (4) also implies that the threshold ξH is exceeded on average once every N(ξ) = 1/P(ξ) waves. Thus, N can be interpreted as the conditional return period of a wave whose crest height exceeds ξ. For weakly nonlinear seas, the probability P(ξ, t) is hereafter described by the third-order Tayfun-Fedele[4] (TF), second-order Tayfun[3] (T), Forristall[8] (F) and the Rayleigh (R) distributions (see Methods section). Note that we are now able to estimate the probability P(ξ) that a wave of the storm has a crest height C larger than ξ = h/H. However, we still need to find in what sea state (of significant wave height H) such a wave most likely occurs.

To do so, we draw on[9,51,52] and express the probability density function (pdf) describing the time at which a wave crest C exceeds a specified or given level h in the interval [t, t + dt] as

The preceding pdf is estimated from data aswhere T is the sea state duration. As an example, consider the Feb 2015 storm. The largest wave with also the largest crest height of h = h = 1.23H = 14 m (see Table 1) is observed in the sea state at the storm peak (H = 12.6m). The pdf p(t|h) estimated for that crest amplitude, and shown in Fig. 8, is very narrow and concentrated around its absolute maximum coincident with the storm peak in agreement with what is expected intuitively. Instead, waves whose crest height exceeds the smaller threshold h = h/2 = 7 m have a pdf that still has its maximum at the storm peak, but it is broader indicating that crest heights exceeding 7 m likely occur also before or after the storm peak.

Figure 8

Probability density function p(t|h) [1/min] of a wave whose crest elevation exceeds the maximum observed height h = h = 1.23H = 14 m in the Feb 2015 storm. Hourly variations (bold line) of significant wave height H are also shown. The vertical dashed line indicates the sea state in which the largest crest was observed. For comparison, the pdf for h = h/2 = 7 m is also shown.

Similarly, the nonstationary occurrence frequency of a wave of the storm whose crest-to-trough wave height exceeds the threshold H as well as the pdf p(t|H) of the sea state in which such waves occur can be both described by the same Eqs ((4), (6)) by simply replacing P(h) with the exceedance probability P(H) appropriate for wave heights of stationary seas. This is hereafter described by the generalized Boccotti (B), Tayfun (T) and linear Rayleigh (R) distributions (see Methods section).

Finally, the second-order nonstationary pdf p(z = η/σ) of wave surface elevations η for storms is defined aswhere p(z = η/σ) is the Tayfun approximation[3,44] for nonlinear stationary sea states, described bywhere and μ = λ3/3, valid for the observed values of skewness coefficient λ3 < 0.6, approximately.

Note that the probability structure of storm-wave characteristics depends upon the time history of wave parameters, say α(t), such as significant wave height, skewness and excess kurtosis. The analyses of the data sets here indicate that the non-stationary distributions are well approximated by the corresponding stationary models of an equivalent sea state with duration equal to that of the actual storm (as if H does not vary in time) evaluated using the weighted average parameters

This follows from modeling the actual storm as if it had an ‘equivalent’ rectangular shape and assuming that the average parameters of the Tayfun and Boccotti models are the same in both the actual and equivalent storms. However, such approximations do not have general validity, and they may not work for other models or data sets. Thus, hereafer we will only consider the non-stationary models based on (4).

The empirical distributions of surface wave elevations for both storms Doris and Feb 2015 are shown in Fig. (9). These are for the most part well described by the non-stationary Tayfun pdf p, which is practically the same as the stationary approximation estimated using the weighted average steepness μ based on Eq. (10). This indicates the dominance of second-order nonlinearities in shaping the sea surface, especially for the more intense Feb 2015 storm.

Figure 9

Probability density function of scaled surface wave elevations η/σ in (left) Doris and (right) Feb 2015 storms. The empirical distributions (□) derived from the total wave population are compared with the non-stationary second-order Tayfun (T) distribution for storms. Dashed curves (G) describe the probability density for a standard unit Gaussian variable.

Figure (10) summarizes the wave statistics for Doris. In particular, the left panel of the figure depicts the empirical distribution for crest heights h/H observed in the total wave population plotted versus the number of waves N(ξ). This is compared against the theoretical predictions of the nonstationary second-order Tayfun (T), third-order Tayfun-Fedele (TF), Forristall (F) as well as the Rayleigh (R) distributions. Although the associated confidence bands on the empirical probabilities noticeably widen over the large waves, the theoretical predictions nonetheless still lie mostly within the same confidence bands. Note that TF is practically the same as T and F as an indication that second-order effects are dominant, whereas the linear R model underestimates the return periods. Similarly, the right panel of the same figure shows the empirical distribution for crest-to-trough wave heights H/H. The observed statistics is well described by both the non-stationary generalized Boccotti (B) and Tayfun (T) models. The small differences among the various models are magnified in Fig. (11), which shows the plots of the normalized crest height h/h and wave height H/H versus the number of waves N. Here, h and H are the crest and crest-to-trough-thresholds exceeded with probability 1/N in a Gaussian sea in accord with the Rayleigh law.

Figure 10

Doris: (left) crest heights h/H and (right) crest-to-trough wave heights H/H versus number of waves N(ξ). Empirical distributions (□) of total population of waves observed in comparison with non-stationary models (T = Tayfun, TF = Tayfun-Fedele, F = Forristall, B = generalized Boccotti and R = Rayleigh distributions). (Light gray lines) approximate confidence bands on observational estimates and (horizontal dashed lines) rogue thresholds[58] for crest (1.25H) and wave heights 2.2H, respectively.

Figure 11

Doris: (top) crest heights h/h and (bottom) crest-to-trough wave heights H/H versus number of waves N(ξ). Empirical distributions (□) of total population of waves observed in comparison with non-stationary models (T = Tayfun, TF = Tayfun-Fedele, F = Forristall, B = generalized Boccotti and R = Rayleigh distributions). (Light gray lines) approximate confidence bands on observational estimates. Amplitude h (H) refers to Rayleigh-distributed crest (crest-to-trough) heights exceeded with probability P(ξ) in Gaussian seas.

Similar conclusions also hold for the wave statistics for the Feb 2015 storm, summarized in Figs (12) and (13). As regard to crests, TF slightly exceeds T, again as an indication that second-order effects are dominant.

Figure 12

Feb 2015: (left) crest heights h/H and (right) crest-to-trough wave heights H/H versus number of waves N(ξ). Empirical distributions (□) of total population of waves observed in comparison with non-stationary models (T = Tayfun, TF = Tayfun-Fedele, F = Forristall, B = generalized Boccotti and R = Rayleigh distributions). (Light gray lines) approximate confidence bands on observational estimates and (horizontal dashed lines) rogue thresholds[58] for crest (1.25H) and wave heights 2.2H, respectively.

Figure 13

Feb 2015: (top) crest heights h/h and (bottom) crest-to-trough wave heights H/H versus number of waves N(ξ). Empirical distributions (□) of total population of waves observed in comparison with non-stationary models (T = Tayfun, TF = Tayfun-Fedele, F = Forristall, B = generalized Boccotti and R = Rayleigh distributions). (Light gray lines) approximate confidence bands on observational estimates. Amplitude h (H) refers to Rayleigh-distributed crest (crest-to-trough) heights exceeded with probability P(ξ) in Gaussian seas.

The wave profiles η with the largest wave crest height observed during Doris and Feb 2015 are shown in the left panel of Fig. (14). In the other panels, we display the El Faro, Draupner, Andrea and Killard rogue wave profiles for comparison[1,2]. In the same figure, the mean sea level (MSL) below the crests is also shown. The estimation of the MSL follows by low-pass filtering the time series of zero-mean surface elevations with a frequency cutoff , where f is the frequency at the spectral peak[53].

Figure 14

Observed wave profiles η/η (thick solid) and mean sea levels (MSL) (thick dashed) versus t/T for (first panel from the left) Doris and Feb 2015, and (following panels from left to right) El Faro, Andrea, Draupner and Killard simulated waves (thin solid) and MSL (thick dashed). Actual measurements (thick solid) and MSLs (thick solid) are also shown for Andrea, Draupner and Killard. Note that the Killard MSL is insignificant and the Andrea MSL is not available. η is the maximum crest height and T is the dominant wave period (see Table 1 and Methods section for definitions).

All six wave profiles are similar, suggesting a common generation mechanism for rogue events. In particular, all cases have sharper crests and rounded troughs and they do not display any secondary maxima or minima. They appear more regular and behave as narrow-banded waves do[45]. In essence, this means that the temporal profile of a relatively large wave observed over a complete phase cycle of 2π displays a single dominant crest or a ‘global’ maximum with no local maxima or minima. In other words, the wave phase monotonously increases over the cycle without any reversals associated with local minima and maxima[45]. These characteristics typical of truly narrow-band waves in every wave cycle are similarly observed but locally in the largest group of waves in a wind-wave field although they are not in the least described by a narrow-banded spectrum. In narrow-band waves, the constructive interference is the primary mechanism for the generation of large displacements in the underlying first-order linear field. The second-order corrections are phase-locked to the linear field such that they always tend to enhance wave crests and flatten the troughs, leading to the basic vertical asymmetry observed in oceanic waves. The process is similar for relatively large waves in a wind-wave field[45].

Further, Doris and Feb2015 both display the characteristics of a dominant wind wave field and show no evident characteristics typical of mixed or crossing seas such as swell overlapping with local wind waves (see e.g. Fig. 4). That may explain the minor set-down observed below the largest waves observed. On the contrary, a set-up below the simulated El Faro and actual Draupner rogue waves is observed, most likely due to the multidirectionality of the respective sea states[2]. Indeed, recent studies showed that Draupner occurred in mixed seas consisting of swell waves propagating at approximately 80 degrees to wind seas[54-56]. Instead, the El Faro sea state showed a very broad directional spreading of energy typical of strong hurricane conditions. The multi-directionality of the two sea states may explain the set-up observed under the large wave[53] instead of the second-order set-down normally expected[57].

Discussion

There is at present no consistent crest-height or wave-height model that works effectively at shallow water depths where kd < π/10 and recent comparisons[15] simply serve to demonstrate this. The wave regime in such shallow waters can only be described by stochastic formulations of highly nonlinear shallow-water equations. Second-order theories or approximations tend to become ineffective at such depths. Our work does not overlap with or extend to such water depths where the second-order theory breaks down. Indeed, we have shown that the theoretical Tayfun[3,6] and Boccotti[9,10] models for crest and wave heights, largely applied to and validated for deep-water waves[4,10] and more recently for mixed/crossing seas[11-13], describe waves reasonably well in intermediate shallow waters also (π/10 < kd < π). So does the second order Forristall model[8].

In particular, we have analyzed actual wave data from ADCP measurements gathered during the passages of two major storms nearshore off Killard Point at the intermediate water depth of approximately 37 m (kd = 0.6–2.5) in 2015 and off the Aran Islands at 45 m depth (kd = 1.36–2.2) in 2017 (see Fig. 1). The observed sea states at the storm peak present the characteristics of a main dominant wind wave field. No evident crossing sea characteristics of overlapping swell and wind components are observed. We have observed time-dependent wave statistics and proposed a novel approach to rationally analyze the non-stationary surface series.

The large wave characteristics measured do not exceed the conventional rogue thresholds[58] h/H = 1.25 and H/H = 2.2 observed in laboratory experiments[15]. In contrast, Draupner, Andrea or Killard rogue waves[1,2], all observed in intermediate water depths, did attain crest heights of approximately 1.6H (see Table 1). Nevertheless, our analysis reveals that the largest waves observed here have characteristics quite similar to those displayed by the El Faro, Andrea, Draupner and Killard rogue waves[1,2] for which second order bound nonlinearities constitute the dominant factor enhancing the linear dispersive focusing of extreme waves.

Moreover, most observed values of the dimensionless depth kd were slightly above (Doris) and below (Feb 2015) the threshold 1.363 above which unidirectional waves are expected to become modulationally unstable[30,31]. The sea states analyzed here were multidirectional, and a carrier wave is modulationally unstable even at depths below that critical value if they are perturbed by appropriate oblique disturbances[32-34]. This type of instabilities are not very likely to appear in theory[32] if kd < 0.5. Nonetheless, rogue waves can be generated by modulational instability, as in unidirectional seas[35,36]. However, in directional seas such as the two considered here, energy can spread directionally and the recurrence of large waves as observed in unidirectional seas is largely attenuated or suppressed[1,2]. Indeed, our statistical analysis indicates that modulational instabilities are ineffective, third-order resonant effects are negligible and second order bound nonlinearities are the dominant factor in shaping the large waves observed.

Our results here also indicate that in shallower water depths, nonlinear dispersion effects intensify[2,21], inducing waves to break more rapidly than in deep waters. As a result, waves cannot breathe as they do not have time to grow and reach higher amplitudes above 1.25H as in deep water. Therefore, whereas the standard rogue thresholds are based on the Rayleigh law appropriate to linear non-breaking Gaussian seas, it makes sense to consider more realistic thresholds and models that account for wave breaking since the latter limits wave growth and impedes the occurrence of rogue waves[2,21,22].

Finally, large waves with higher and sharper crests do not display any secondary maxima or minima. They appear more regular or “narrow banded” than relatively low waves, and their heights and crests do not often violate the Miche–Stokes type upper limits[59]. Our results also suggest that third-order resonant nonlinearities do not affect the surface statistics in any discernable way, in agremeent with recent rogue wave studies[1,2]. Indeed, our analysis reveals that fourth order cumulants are negligible. As a consequence, the sea states analyzed here have nothing specially rogue about them.

Methods

ADCP measurements

A Teledyne Sentinel V acoustic Doppler current profiler (ADCP) was deployed off Killard Point, Ireland (upper inset of Fig. 1) during Spring 2015 and off the Aran Islands, Ireland (lower inset of Fig. 1) during Spring 2017, to measure wave events. The instrument itself was secured in a frame to ensure it stayed in position and to prevent damage (see right panel of Fig. 1). The frame and instrument were placed at rest at the sea bottom, at an average depth of 37 m (Killard Point) and 45 m (Aran Islands). The four slant beams made a 25° angle with the vertical, so that at the surface, the maximum distance between beams was approximately 35 m (Killard Point) and 42 m (Aran). This ADCP operates by emitted sound pulses in five beams (four slanted and one vertical) and using the Doppler effect to measure the movement of sound scatterers such as plankton and small particulates, within these beams[60]. Each beam divides the water column into 38 bins, separated by 122 cm. Because of hardware limitations, data were sampled at 2 Hz similarly to standard wave measurements gathered at oil platforms[14]. Drawing on[61], the sampling error on estimating crest and wave heights, the so-called quantization error[62,63], is approximately 1–2%. This is mitigated by correcting for crest amplitudes by quadratically interpolating the sampled crests, as in recent stereo-measurements of the ocean surface[62,63].

We correct the resulting data sets of echo intensity and velocity measurements for pitch, roll, and heading of the instrument in the water, and convert from instrument coordinates (a radial set) to geographical coordinates (North, East, Up) using standard transformations[64,65]. We interpolate the data to find the position of maximum intensity, corresponding to the location of the free surface[64].

Directional spectrum

We estimate the directional spectrum from the free-surface profiles of the four slanted beams using the DIWASP toolbox[66]. Consequently, we are able to determine angular spreading and directionality of a sea state. DIWASP uses a number of methods to estimate the directional spectrum from the cross-power spectrum of the data: direct Fourier transform method (DFTM), extended maximum likelihood method (EMLM), iterated maximum likelihood method (IMLM), extended maximum entropy principle (EMEP), and Bayesian direct method (BDM)[66].

However, the estimations are not perfect due to limited information and unknown factors. The EMLM spectrum is often more directionally-diffused with a lower peak. The EMEP spectrum produces a good directional spreading. However, although the peak is higher than EMLM, it is below the desired result[67]. EMEP and BDM can give very similar spreading results, but their peak values often differ significantly[68,69]. EMEP can calculate bi-directionality, while BDM is less sensitive to probe layout and more robust against errors[68]. In our analysis, we consider the BDM spectrum.

Wave parameters

The significant wave height H is defined as the mean value H1/3 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 omnidirectional surface spectrum as H ≈ 4σ, where is the standard deviation of surface elevations, denotes spectral moments. Further, S(f,θ) is the directional wave spectrum with θ as the direction of waves of frequency f, and the cyclic frequency is ω = 2πf. In this paper, we use the spectral-based estimate, which is 5–10% larger than the actual H1/3 estimated from the actual time series.

The dominant wave period T = 2π/ω follows from the cyclic frequency ω of the spectral peak and T0 is the observed mean zero-crossing wave period. For Gaussian seas, this is equal to 2π/ω0, with . The associated wavelength L0 = 2π/k0 follows from the linear dispersion relation , with d the water depth. The ‘mean’ or central frequency ω of the spectrum is defined as ω = m1/m0 [3] and the associated mean period T is 2π/ω. Theoretical NB steepness[3,49] is defined as μ = kσ, where k is the mean wavenumber corresponding to ω via the linear dispersion relation

The group velocitywhere ω' is the first derivative of the mean frequency with respect to the wavenumber k.

The spectral bandwidth ν = (m0m2/m12 − 1)1/2 gives a measure of the frequency spreading. For unimodal directional spectra, as those analyzed in this study, the angular spreading , where the angular spreading function and is the mean direction. In general, . Furthermore, we define R = σ2/2ν2 as a dimensionless measure of the directionality of a sea state[37,70]. In terms of ω, k and q = kd based on the spectral centroid, the directional factor is given by[40]and the depth factor[40]whereand

The nonlinear interaction coefficient is[40]where is the phase velocity in shallow waters,accounts for the wave-induced setdown[40], and

Statistical parameters

The normalized covariance function of zero-mean surface displacement η(t) is defined as . An alternative measure for the spectral bandwidth is given by the Boccotti parameter ψ* = |ψ(τ*)|, which is the absolute value of the first minimum of ψ at τ = τ* [9] and the corresponding second derivative with respect to τ.

Skewness coefficient λ3 and excess kurtosis λ40 of the zero-mean surface elevation η(t) are given by

Here, overbars imply statistical averages and σ is the standard deviation of surface elevations. Clearly, the wave steepness μ = λ3/3 [6] relates to the skewness coefficient λ3 of surface elevations. For third-order nonlinear (NB) random seas the excess kurtosis[37,38]comprises a dynamic component λ40 due to nonlinear quasi-resonant wave-wave interactions[39] and a Stokes bound harmonic contribution λ40 [41]. Drawing on[40] and using parameters based on ω, k and q = kd, wave skewness and bound excess kurtosis for narrowband (NB) waves in intermediate water are given bywherewhere

contributes positively to the wave induced setdown[40] due to the directional spread of waves. Assuming that the linear crest heights (ξ0) scaled with surface rms are Rayleigh-distributed, the mean wave-induced setdown in the still water level is given by 〈ξ02〉μΔ = 2μΔ, where brackets denote statistical average.

When Δ = 0, some algebra shows that λ3, is the same as the original Marthinsen-Winterstein formulation[49,71], developed nearly three decades ago in the form:where

The coefficients D1 and D2 arise from the frequency-difference and frequency-sum terms of second-order wave-wave interactions. Note that D1 = 2Δ and D2 = 2α exactly in unidirectional waves for which Δ = 0. Unfortunately, Eq. (13) are not valid in relatively shallow water depths as second and third-order terms of the associated Stokes expansion can be larger than the linear term (see Eq. (A18) in[41]) because of the divergent nature of α and β. Thus, the relative validity of the preceding results essentially assumes the constraints αμ ≤ 1 and βμ/α ≤ 1. These are satisfied for all seas states of both storms studied here.

Miche-Stokes upper limit

In Gaussian seas, surface displacements and thus wave and crest heights have unbounded ranges. In reality, surface elevations are neither exactly Gaussian nor unbounded. And, the crest-to-trough height of a large wave whose steepness approaches the Stokes limiting steepness is unlikely to exceed an upper bound. For long-crested waves in transitional water depths, Miche[59] approximated this upper bound as

where σ is the standard deviation of the sea state, k is the wavenumber and d the water depth. Following Tayfun[45] the corresponding Miche-Stokes limit for crest heights is estimated as

The Miche-Stokes limit can be rewritten as a function of wave period T via the linear dispersion relation. Finally, note that in realistic oceanic seas, nonlinear wave dispersion is effective in limiting the wave growth as a precursor to breaking[21-23]. Thus, in wave fields generated by intense storms, the onset of breaking can occur well below the preceding Miche-Stokes type upper bounds[22,24,25].

The Tayfun-Fedele model for crest heights

We define P(ξ) as the probability that a wave crest observed at a fixed point of the ocean in time exceeds the threshold ξH. For weakly nonlinear nonlinear seas, this probability can be described by the third-order Tayfun-Fedele model[4],where ξ0 follows from the quadratic equation ξ = ξ0 + 2μξ02. Here, μ = λ3/3 is the Tayfun steepness: it represents an integral measure of wave steepness and relates to second-order bound nonlinearities. The parameter Λ = λ40 + 2λ22 + λ04 is a relative measure of third-order nonlinearities expressed in terms of the fourth-order cumulants λ of surface elevation η and its Hilbert transform [4]. In particular, and . In this study, Λ is approximated solely in terms of the excess kurtosis as Λappr = 8λ40/3. This approximation follows from the NB relations between cumulants[43,72] λ22 = λ40/3 and λ04 = λ40, valid as the spectral bandwidth ν tends to zero. Numerical computations[1] indicate that Λ ≈ Λappr with an error of about 3% in wave fields where second-order nonlinearities are dominant, in agreement with observations[35,73].

For second-order seas, Λ = 0 and P in Eq. (19) leads to the Tayfun wave-crest distribution[3,6]where ξ = ξ0 + 2μξ02. For Gaussian seas, ξ0 = ξ since μ = 0 and Λ = 0, and P reduces to the Rayleigh distribution

Note that the Tayfun distribution represents an exact theoretical result for large second-order wave crest heights and it depends solely on the steepness parameter defined as μ = λ3/3[6].

The Forristall’s Weibull model for crest heights

The exceedance probability for crest heights is given by[8]where α = 0.3536 + 0.2561S1 + 0.0800U, β = 2 − 1.7912S1 − 0.5302U + 0.284U2 for multi-directional (short-crested) seas. Here, S1 = 2πH/(gT2) is a characteristic wave steepness and the Ursell number U = H/(k2d3), where k is the wavenumber associated with the mean period T.

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[10]

and the Boccotti parameters[9] ψ* and are defined above in the section where statistical parameters are described. For Gaussian seas (Λ = 0), the original Boccotti[9] model is recovered

The Tayfun model for wave heights is given by[4,74]where r = r(T/2) is the value of the envelope r(t) of the covariance ψ(t) at t = T/2. Finally we note that as spectral bandwidth ν tends to zero, all three parameters ψ*, and r tend to unity, and the Boccotti and Tayfun distributions both reduce to the Rayleigh distribution given by