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.

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][17][18][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 p 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 < k p d < 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.

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  Table 1.
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 sea 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 σ j 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 s /T sea successive sea states of the same duration T sea and variances σ 1 2 ,σ 2 2 , ..., σ j 2 ,σ j + 1 2 , ..., where D s is the storm duration (70 hours). Then, N − 1 sampled values of V follow as V j = σ j + 1 /σ j − 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 sea approaches smaller values.
This process ensured that resulting statistics are robust to variations in T sea 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 k p d during the two storms. Surface spectra observed were broadbanded so that ω p ≈ 0.8ω m , while k p ≈ 0.7k m .
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 s = 4σ. For comparison, the variation of actual significant height H 1/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 s . The actual mean zero-up-crossing wave period T 0 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 .
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 d (ω, θ)/σ 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 (5) displays the scatter diagrams of crest heights h/H s versus corresponding wave periods T/T 0 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][22][23] . Thus, the onset of wave-breaking can occur well below the Miche-Stokes upper limit 22,[24][25][26][27] .  30,31 .
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 .

Modulational instability in intermediate or transitional water depths.
In the Feb 2015 storm, the depth coefficient k p d, 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][33][34] . Nevertheless, it is also recognized that instabilities ensuing from such disturbances are not likely to occur for values of k p d < 0.5 32 . So, it is plausible that rogue waves could be generated by modulational instability, as in unidirectional seas 35 , 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, BFI k 2 / p σ ν = defines the Benjamin-Feir index in deep water at the spectral peak, ν the spectral bandwidth, i 1 = − and Im(x) denotes the imaginary part of x, ω p and k p the dominant spectral frequency and wavenumber. Depth effects on wave directionality, measured by R, are represented by R S = β S R by way of the factor β S , and α S is the depth factor. The latter two depend on the dimensionless depth k p d (see Methods section). In the deep-water limit, both α S and β S become 1. The maximum of dynamic excess kurtosis is well approximated by 39   www.nature.com/scientificreports www.nature.com/scientificreports/  48. In the focusing regime (0 < R S < 1 and α S > 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 k p d < 1.363. Indeed, in 1-D waves modulational instability disappears below that critical threshold and α S < 0. As a result, wave dynamics becomes of defocusing type and excess kurtosis is negative. In particular, it reaches a minimum at t c 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 α S > 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 S > 1. A rogue wave regime is a more likely occurrence only if the surface spectrum is sufficiently narrow-banded (R S < 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 α S > 0. However, all their sea states are directionally spread and characterized with mostly negative excess kurtosis since R S > 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 m yield similar negligible estimates of the dynamic excess kurtosis.
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 λ ij 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 d measuring third-order quasi-resonant wave-wave interactions and a bound contribution λ 40 b 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 m , rather than k p , tend to yield more favorable comparisons with deep-water observations or theories 6,44 . Definitions based on k p lead to predictions that noticeably underestimate the observed and/or theoretical statistics of broadband waves since k p < k m 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 m or k p 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 μ m and μ p based on the mean and dominant wavenumbers k m and k p , respectively, and (right) the observed fourth-order Λ coefficient, its approximation Λ appr and the NB predictions 37,41 Λ m and Λ p based on k m and k p for Doris. The same comparisons are also reported in the bottom panels for Feb 2015. Clearly, the two NB predictions μ m overestimate and μ p 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 s . Drawing on 9,50-52 , the storm is modeled as a non-stationary continuous sequence of sea states of duration dt, and dt/T 0 (t) is the number of waves in the sea state, where T 0 (t) is the time-changing mean zero-crossing wave period. Consider now a wave of the storm and define the probability P ns (ξ) that its crest height C exceeds the threshold h = ξH s as observed at a fixed location where H s = 4σ. Equivalently, this is the probability of randomly picking a wave crest whose height C exceeds the threshold ξ = h/H s from the non-stationary time series observed at a fixed point of the storm. Then, The definition of P ns 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 s sea states of duration T sea = D s /N s . In each sea state S j centered at t = t j , waves are sampled and their crest elevations (h) are normalized with the local significant wave height as h/H s (t j ) and put all in the same population . Then, the empirical exceedance probability P ns (ξ) 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, N t dt EX t dt d t ( , ) () 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 s (t) in the same sea state. Then, P ns 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 ns is estimated from data as the weighted average where N w (t j ) is the number of waves sampled in the sea state S j . Equation (4) also implies that the threshold ξH s is exceeded on average once every N h (ξ) = 1/P ns (ξ) waves. Thus, N h 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 ns (ξ) that a wave of the storm has a crest height C larger than ξ = h/H s . However, we still need to find in what sea state (of significant wave height H s ) 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 as  Table 1) is observed in the sea state at the storm peak (H s = 12.6m). The pdf p(t|h) estimated for that crest amplitude, and shown in Fig. 8 www.nature.com/scientificreports www.nature.com/scientificreports/ 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 max /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.
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 ns (z = η/σ) of wave surface elevations η for storms is defined as where p(z = η/σ) is the Tayfun approximation 3,44 for nonlinear stationary sea states, described by ( ) 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 s 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 ns , which is practically the same as the stationary approximation estimated using the weighted average steepness μ based on Eq. (10). This . 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  www.nature.com/scientificreports www.nature.com/scientificreports/ 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 s . 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 R and wave height H/H R versus the number of waves N h . Here, h R and H R are the crest and crest-to-trough-thresholds exceeded with probability 1/N h in a Gaussian sea in accord with the Rayleigh law.
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.
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 f f /2 c p , where f p is the frequency at the spectral peak 53 . 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][55][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][12][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 (k p d = 0.6-2.5) in 2015 and off the Aran Islands at 45 m depth (k p d = 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 s = 1.25 and H/H s = 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 s (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 k p d 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][33][34] . This type of instabilities are not very likely to appear in theory 32 if k p d < 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 s 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.  Table 1 and Methods section for definitions).  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. denotes spectral moments. Further, S d (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 H 1/3 estimated from the actual time series.
The dominant wave period T p = 2π/ω p follows from the cyclic frequency ω p of the spectral peak and T 0 is the observed mean zero-crossing wave period. For Gaussian seas, this is equal to 2π/ω 0 , with m m /  〉μ m Δ = 2μ m Δ, where brackets denote statistical average. When Δ ST = 0, some algebra shows that λ 3,NB is the same as the original Marthinsen-Winterstein formulation 49,71 , developed nearly three decades ago in the form: The coefficients D 1 and D 2 arise from the frequency-difference and frequency-sum terms of second-order wave-wave interactions. Note that D 1 = 2Δ and D 2 = 2α exactly in unidirectional waves for which Δ ST = 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 αμ m ≤ 1 and βμ m /α ≤ 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 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][22][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 s . For weakly nonlinear nonlinear seas, this probability can be described by the third-order Tayfun  where ξ 0 follows from the quadratic equation ξ = ξ 0 + 2μξ 0 2 . 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 λ nm of surface elevation η and its Hilbert transform η 4 . In particular, . 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 TF in Eq. (19) leads to the Tayfun wave-crest distribution 3,6 ξ ξ = − P ( ) exp ( 8 ) , (20) T 0 2 where ξ = ξ 0 + 2μξ 0 2 . For Gaussian seas, ξ 0 = ξ since μ = 0 and Λ = 0, and P TF reduces to the Rayleigh distribution R 2 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  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   where r m = r(T m /2) is the value of the envelope r(t) of the covariance ψ(t) at t = T m /2. Finally we note that as spectral bandwidth ν tends to zero, all three parameters ψ * , ̈⁎ ψ and r m tend to unity, and the Boccotti and Tayfun distributions both reduce to the Rayleigh distribution given by R 2

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.