The sinking of the El Faro: predicting real world rogue waves during Hurricane Joaquin

We present a study on the prediction of rogue waves during the 1-hour sea state of Hurricane Joaquin when the Merchant Vessel El Faro sank east of the Bahamas on October 1, 2015. High-resolution hindcast of hurricane-generated sea states and wave simulations are combined with novel probabilistic models to quantify the likelihood of rogue wave conditions. The data suggests that the El Faro vessel was drifting at an average speed of approximately 2.5 m/s prior to its sinking. As a result, we estimated that the probability that El Faro encounters a rogue wave whose crest height exceeds 14 meters while drifting over a time interval of 10 (50) minutes is ~1/400 (1/130). The largest simulated wave is generated by the constructive interference of elementary spectral components (linear dispersive focusing) enhanced by bound nonlinearities. Not surprisingly then, its characteristics are quite similar to those displayed by the Andrea, Draupner and Killard rogue waves.

stability 10-m wind speed U 10 and direction are shown in the bottom-panels, respectively. The red vertical lines delimit the 1-hour interval during which the El Faro vessel sank.
The 1-hour sea state experienced by El Faro at and around the time and location of sinking had a significant wave height of H s ≈ 9 m and the maximum wind speed was U 10,max = 51 m/s. Waves were multidirectional (short-crested) as indicated by the large values of both the spectral bandwidth ν and angular spreading σ θ as shown in Fig. 2.
In Table 1 we report the metocean parameters of the El Faro sea state in comparison to those of the Draupner, Andrea and Killard rogue sea states 7 . Note that the four sea states have similar metocean characteristics. However, El Faro is a steeper sea state as the mean wavelengh L 0 is shorter than those observed in the other three cases.
Statistical properties of Hurricane Joaquin-generated seas. The relative importance of ocean nonlinearities can be measured by integral statistics such as the coefficients of skewness λ 3 and excess kurtosis λ 40 of the zero-mean surface elevation η(t). The skewness is a measure of 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 [8][9][10] . The excess kurtosis indicates whether the tails of the distribution of surface elevations is heavy-or light-tailed relative to a Gaussian distribution. It comprises a dynamic component λ d 40 measuring third-order quasi-resonant wave-wave interactions and a bound contribution λ b 40 induced by both second-and third-order bound nonlinearities [8][9][10][11][12][13] .
In deep waters, the dynamic kurtosis 14 depends on the Benjamin-Feir index BFI and the parameter R, a dimensionless measure of the multidirectionality of dominant waves 11,14,15 . For unidirectional (1D) waves R = 0. The bottom panel of Fig. 2 displays the hourly variations of the directional factor R during Hurricane Joaquin near the location where El Faro sank. Around the peak of the hurricane, the generated sea states are quite multidirectional (short-crested) as R > 1. As wave energy also spreads directionally, nonlinear focusing due to modulational instability effects diminishes 14,[16][17][18] and becomes essentially insignificant under such realistic oceanic conditions 7,14,19,20 .
The top panel of Fig. 3 displays the hourly variation of the Tayfun steepness μ (solid line) with associated bounds (dashed lines). The coefficient of excess kurtosis λ 40 mostly due to bound nonlinearities is shown in the center panel and the associated Λ parameter at the bottom. The red vertical lines delimit the 1-hour interval during which the El Faro vessel sank.
In Table 1 we compare the statistical parameters of the El Faro sea state and the Draupner, Andrea and Killard rogue sea states (from ref. 7). Note that the El Faro sea state has the largest directional spreading. Moreover, for all the four sea states the associated BFI are less than unity and the maximum dynamic excess kurtosis is of O(10 −3 ) and thus negligible in comparison to the associated bound component. Thus, third-order quasi-resonant interactions, including NLS-type modulational instabilities play an insignificant role in the formation of large waves 7

Higher Order Spectral (HOS) simulations of the El Faro sea state. We have performed
Higher-Order pseudo-Spectral (HOS) simulations 25,26 of the El Faro sea state over an area of 4 km × 4 km for a duration of 1 hour (see Methods section for a description of the numerical method). The initial wave field conditions are defined by the WAVEWATCH III hindcast directional spectrum S(f, θ) around the time and region of the El Faro sinking as shown in Fig. 4. This is the result of a balance of the energy fluxes due to wind input (S in ), exact four-wave resonance nonlinearities (S nl ) and dissipation due to wave breaking (S ds ). Wind gustiness and currents are not modeled. Our WW3 hindcast indicates that the flux S in is balanced out by S ds . In particular, around the spectral peak 60% of wind input is lost to dissipation. This offset increases away from the peak. Any wave growth associated with S in + S ds and S nl is accounted for in the WW3 model. It is the wave growth associated with quasi-resonant and bound harmonics nonlinear effects that is not modeled. In our study, we exploit the HOS wave solver to simulate the El-Faro sea state by accounting for quasi-resonant and bound nonlinearities up to fourth order in wave steepness. An estimate of the most likely rogue wave amplitude is then provided as discussed below. Note that both wind input and wave breaking are somewhat modeled in our HOS simulations as these are initialized with the WW3 spectrum. Clearly, our analysis suggests future studies on the relative importance of possible effects such wind gustiness 27 and wave breaking 28, 29 on the HOS model results.
The wavenumber-frequency spectrum S(k, ω) estimated from the HOS simulations is shown in Fig. 5. Here, dashed lines indicate the theoretical dispersion curves related to the first-order (1 st ) free waves as well as the second (2 nd ) and third-order (3 rd ) bound harmonic waves. The HOS predictions indicate that second-order nonlinearities are dominant with a weak effect of third-order nonlinear bound interactions, in agreement with recent studies of rogue sea states 7 . It appears that fourth-order effects are insignificant.
The wave skewness and kurtosis rapidly reach a steady state after a few (mean) wave periods as an indication that third-order quasi-resonant wave-wave interactions are negligible in agreement with theoretical predictions 14 and simulations 7 . Note that the theoretical narrowband (NB) predictions slightly overestimate the simulated values for skewness and excess kurtosis (see Table 1). The same trend is also observed in recent studies on rogue waves 7 . This is simply because NB approximations do not account for the directionality and finite spectral bandwidth of the El Faro wave spectrum.
Occurrence frequency of a rogue wave by a fixed observer: the return period of a wave whose crest height exceeds a given threshold. To describe the statistics of rogue waves encountered by an observer at a fixed point of the ocean surface, we consider the conditional return period N h (ξ) of a wave whose crest height exceeds the threshold h = ξH s , namely where P(ξ) is the probability or occurrence frequency of a wave crest height exceeding ξH s as encountered by a fixed observer. In other words, P(ξ) is the probability to randomly pick from a time series observed at a fixed point of the ocean a wave crest that exceeds the threshold ξH s . Equation (1) also implies that the threshold ξH s , with H s = 4σ, is exceeded on average once every N h (ξ) waves. For weakly nonlinear random seas, the probability P is hereafter described by the third-order Tayfun-Fedele 9 (TF), second-order Tayfun 8 (T), second-order Forristall 30 (F) and the linear Rayleigh (R) distributions (see Methods section).
Our statistical analysis of HOS wave data suggests that second-order effects are the dominant factors in shaping the probability structure of the El Faro sea state with a minor contribution of excess kurtosis effects. Such dominance is seen in Fig. 6, where the HOS numerical predictions of the conditional return period N h (ξ) of a crest exceeding the threshold ξH s are compared against the theoretical predictions based on the linear Rayleigh (R), second-order Tayfun (T) and third-order (TF) models from Eq. (17). It is noted that the HOS predictions are based on a sample population of 10 6 crests. In particular, N h (ξ) follows from Eq. (1) as the inverse 1/P(ξ) of the empirical probabilities of a crest height exceeding the threshold ξH s . An excellent agreement is observed between simulations and the third-order TF model up to crest amplitudes h/H s ~ 1.5. For larger amplitudes, the associated confidence bands of the estimated empirical probabilities widen, but TF is still within the bands. Donelan and Magnusson 31 suggest that the TF model agrees with the Andrea rogue wave measurements up to h/H s ~ 1.1, concluding that TF is not suitable to predict larger rogue crest extremes (see their Fig. 7 in ref. 31). Unfortunately, their analysis is based on a much smaller sampled population of ~10 4 crest heights and they do not report the confidence bands associated with their probability estimates, nor they provide any parameter values to validate their data analysis. The deviation of their data from the TF model is most likely due to the relatively smaller population of crests observed. Note also that TF slightly exceeds both the T and F models as an indication that second-order effects are dominant, whereas the linear R model underestimates the return periods.
For both third-and fourth-order nonlinearities, the return period N r of a wave whose crest height exceeds the rogue threshold 1.25H s ≈ 11 m 32 is nearly N r ~ 10 4 for the El Faro sea state and for the simulated Andrea, Draupner and Killard rogue sea states 7 . This is in agreement with oceanic rogue wave measurements 23 , which yield roughly the same return period. Similarly, recent measurements off the west coast of Ireland 33 yield N r ~ 6 · 10 4 . In contrast, N r ~ 3 · 10 5 in a Gaussian sea.
Note that the largest simulated wave crest height exceeds the threshold 1.6H s ≈ 14 m (see Table 1). This is exceeded on average once every 10 6 waves in a time series extracted at a point in third-and fourth-order seas and extremely rarely in Gaussian seas, i.e. on average once every 10 9 waves. This implies that rogue waves observed at a fixed point of the ocean are likely to be rare occurrences of weakly random seas, or Tayfun sea states 34 . Our results clearly confirm that rogue wave generation is the result of the constructive interference (focusing) of elementary waves enhanced by bound nonlinearities in agreement with the theory of stochastic wave groups developed by Fedele and Tayfun (2009) 10   simulated rogue waves indicate that such waves were nearly incipient breaking 28,29,44 suggesting that larger rogue events are less likely to occur 21,44 . The saturation of the crest height is mainly due to the nonlinear dispersion and it is an energy limiter for rogue waves.
The four wave profiles are very similar suggesting a common generation mechanism of the rogue events. The manner waves are generated by Hurricane Joaquin or the northerly storm of the Draupner, Andrea and Killard sea states, all four waves and their statistics cannot differ in a fundamental way from each other as the    Further, we observe a set-up below the simulated El Faro rogue wave, most likely due to the multidirectionality of the sea state. A set-up is also observed for the actual Draupner rogue wave. Indeed, recent studies showed that Draupner occurred in a crossing sea consisting of swell waves propagating at approximately 80 degrees to the wind sea 46,47 . This would explain the set-up observed under the large wave 43 instead of the second-order set-down normally expected 48 .
Space-time statistics of the sea state encountered by El Faro before sinking. The largest crest height of a wave observed in time at a given point of the ocean represents a maximum observed at that point. Clearly, the maximum wave surface height observed over a given area during a time interval, i.e. space-time extreme, is much larger than that observed at a given point. Indeed, in relatively short-crested directional seas such as those generated by hurricanes, it is very unlikely that an observed large crest at a given point in time actually coincides with the largest crest of a group of waves propagating in space-time. In contrast, in accord with Boccotti's (2000) QD theory 35 , it is most likely that the sea surface was in fact much higher somewhere near the measurement point.
Space-time wave extremes can be modeled stochastically 3, 4 drawing on the theory of Euler Characteristics of random fields [49][50][51] and nonlinear wave statistics 14 . In the following, we present the Fedele's Space-Time (FST) stochastic model for the prediction of space-time extremes 3 that accounts for both second and third-order nonlinearities 5 . Fedele's work 3, 5 considers a 3-D non-Gaussian field η(x, y, t) in space-time over an area A for a time period of D (see Fig. 8). The area cannot be too large since the wave field may not be homogeneous. The duration should be short so that spectral changes occurring in time are not significant and the sea state can be assumed as stationary. Then, the third-order nonlinear probability ξ P A D ( ; , ) nl FST ( ) that the maximum surface elevation η max over the area A and during the time interval D exceeds the generic threshold ξH s is described by 5   The statistical interpretation of h q is as follows: the maximum surface height η max observed within the area A during D exceeds the threshold h q only in qN realizations of the above mentioned ensemble of N sea states.
Note that for large areas, i.e. L 0 , our FST model as any other similar models available in literature 47, 53-56 will overestimate the maximum surface height over an area and time interval because they all rely on Gaussianity. This implies that there are no physical limits on the values that the surface height can attain as the Gaussian model does not account for the saturation induced by the nonlinear dispersion 21 of ocean waves or wave breaking. Thus, This point is elaborated further and demonstrated explicitly by way of the results displayed in Fig. 9. Here, the theoretical (FST) ratio h h / FST T as a function of the area width  L / 0 is shown for the El Faro, Draupner and Andrea sea states respectively. The FST ratios for Draupner and Andrea are estimated using the European Reanalysis (ERA)-interim data 5 . For comparisons, the empirical ST ratio from the El Faro HOS simulations together with the experimental observations at the Acqua Alta tower 4 are also shown. Recall that h FST is the mean maximum surface height expected over the area  2 during a sea state of duration D = 1 hour and h T is the mean maximum surface height expected at a point. Clearly, the theoretical FST ratio for El Faro fairly agrees with the HOS simulations for small areas ( ≤  L 0 ), whereas it yields overestimation over larger areas. We argue that the saturation of the HOS FST ratio over larger areas is an effect of the nonlinear dispersion which is effective in limiting the wave growth as a precursor to breaking 21,44 .
Note that the FST ratios for all the three sea states are nearly the same for ≤  L 0 . These results are very encouraging as they suggest possible statistical similarities and universal laws for space-time extremes in wind sea states 5 . Moreover, for ∼  L 0 the mean wave surface maximum expected over the area is 1.35 times larger than that expected at a point in agreement with Acqua Alta sea observations 4 .

The occurrence frequency of a rogue wave by the El Faro vessel. The data suggests that the El Faro
vessel was drifting at an average speed of approximately 2.5 m/s prior to its sinking. This is considered in our analysis as follows. First, define the two events R = "El Faro encounters a rogue wave along its navigation route" and S = "El Faro sinks". We know that the event S happened. As a result, one should consider the conditional probability Here, Pr[S] is the unconditional probability of the event that El Faro sinks. This could be estimated from worldwide statistics of sunk vessels with characteristics similar to El Faro. Pr[S|R] is the conditional probability that El Faro sinks given that the vessel encountered a rogue wave. This probability can be estimated by Monte Carlo simulations of the nonlinear interaction of the vessel with the rogue wave field.
Our rogue wave analysis provides an estimate of the unconditional probability Pr[R] that El Faro encounters a rogue wave along its navigation or drifting route by means of the exceedance probability, or occurrence frequency P e (h). This is the probability that a vessel along its navigation path encounters a rogue wave whose crest height exceeds a given threshold h. The encounter of a rogue wave by a moving vessel is analogous to that of a big wave that a surfer is in search of. His likelihood to encounter a big wave increases if he moves around a large area instead of staying still. This is a space-time effect which is very important for ship navigation and must be accounted for 3, 57-60 . The exceedance probability P e (h) is formulated as follows. Consider a random wave field whose surface elevation at a given point (x, y) in a fixed frame at time t is η (x, y, t). Consider a vessel of area A that navigates through the wave field at a constant speed V along a straight path at an angle β with respect to the x axis. Define also (x e , y e ) as a cartesian frame moving with the ship. Then, the line trajectories of any point (x e , y e ) of the vessel in the fixed frame are given by where for simplicity we assume that at time t = 0 the center of gravity of the vessel is at the origin of the fixed frame.
The surface height η c (t) encountered by the moving vessel, or equivalently the surface fluctuations measured by a wave probe installed on the ship, is c e e e e If η is a Gaussian wave field homogeneous in space and stationary in time, then so is η c with respect to the moving frame (x e , y e , t). The associated space-time covariance is given by c e e c e e x y e and k is the wavenumber associated with the frequency f by way of the wave dispersion relation. As a result of the Doppler effect, the encountered, or apparent frequency is given by [57][58][59][60] e and S(f, θ) is the directional wave spectrum of the sea state. Note that when the vessel moves faster than waves coming from a direction θ, the apparent frequency f e < 0 and for an observer on the ship waves appear to move away from him/her. In this case, the direction of those waves should be reversed 57 , i.e. θ = θ + π, and f e set as positive.
The spectral moments m ijk e ( ) of the encountered random field readily follow from the coefficients of the Taylor series expansion of Ψ(X, Y, T) around (X = 0, Y = 0, T = 0). In particular, The nonlinear space-time statistics can then easily processed by using the encountered spectral moments m ijk e ( ) using the FST model 3,5 , which is based on Eq. (2) as described above. Note that for generic navigation routes the encountered wave field η c is a non-stationary random process of time. Thus, the associated spectral moments will vary in time. The space-time statistics can be still computed by first approximating the route by a polygonal made of piecewise straight segments along which the random process η c is assumed stationary. Figure 10 illustrates the HOS and theoretical predictions for the normalized nonlinear threshold h n /H s exceeded with probability 1/n, where n is the number of waves. In particular, consider an observer on the vessel moving along the straight path Γ spanned by El Faro drifting against the dominant sea direction over a time interval of 10 minutes. In space-time the observer spans the solid red line shown in Fig. 8. In this case, he has a probability P e ~ 3 · 10 −4 to encounter a wave whose crest height exceeds the threshold 1.6H s ≈ 14 m (blue lines), and the expected spatial shape is shown in Fig. 11. If we also account for the vessel size (base area A = 241 × 30 m 2 ), in space-time El Faro spans the volume of the slanted parallelepiped V a shown in Fig. 8. In this case, the exceedance probability P e (V a ) further increases to 1/400 (black lines in Fig. 10). Note that if the vessel would be anchored at a location for the same duration, in spacetime it would span instead the volume of the vertical parallelepiped V c shown in the same Figure. Note that the two parallelepipeds cover the same space-time volume A × D, with the base area A and height D = 10 min. For the case of the anchored vessel, the associated exceedance probability P e (V c ) is roughly the same as P e (V a ) since El Faro was drifting at a slow speed. Larger drift speeds yield larger P e (V a ) since the vessel encounters waves more frequently than if it was anchored, because of the Doppler effect 58,59 . Moreover, the drifting vessel covers the strip area (1500 × 30 m 2 ) in the 10-minute interval and the associated space-time volume is that of the parallelepiped V b shown in Fig. 8, which has a larger volume than that of V a . As a result, the occurrence frequency P e (V b ) of a rogue wave associated with V b is larger and it increases to ~1/100 (see red lines in Fig. 10). However, El Faro does not visit the entire volume V b , but it only spans the smaller volume V a . Thus, the conditional probability P e (V a |V b ) that the drifting El Faro encounters a rogue wave given that a rogue wave occurred over the larger spacetime volume V b is P e (V a )/P e (V b ) ~ 1/4. Furthermore, a fixed observer has a much lower probability P e ~ 10 −6 to pick randomly from a time series extracted at a point a wave whose crest height exceeds 1.6H s (see Fig. 6, TF model, black solid line). Finally, we observed that the exceedance probability P e (V a ) for the drifting El Faro does not scale linearly with time because of nonlinearities that reduce the natural dispersion of waves. Indeed, assuming that El Faro drifts over a time interval 5 times longer (50 minutes), P e (V a ) just increases roughly by 3 times, ~1/130.

Discussions
Our present studies open a new research direction on the prediction of rogue waves during hurricanes using the WW3 wave model combined with HOS simulations and the new stochastic FST model 5 for the prediction of space-time wave extremes 3,4 . Any wave growth associated with wind stresses, dissipation due to wave breaking and exact nonlinear resonant four-wave interactions is accounted for in the WW3 model. It is the wave growth associated with quasi-resonant and bound harmonics nonlinearities that is not modeled. Such nonlinear effects can locally increase the wave amplitude over the expected values of the WW3 simulations. In our analysis, quasi-resonant and bound nonlinearities are modeled by way of a HOS wave solver that simulated the sea state around the time and location of the accident. The HOS simulations provided an estimate of the most likely rogue wave amplitude over a given area and time interval indicating that bound nonlinearities are dominant, in agreement with recent rogue-wave studies 7 . Our analysis also suggests new studies on the possible effects of factors such as wind gustiness 27 and wave breaking 28, 29 on generating rogue waves and associated statistics.

Methods
Wave parameters. The significant wave height H s is defined as the mean value H 1/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 wave omnidirectional spectrum as H s ≈ 4σ, where σ = m 0 is the standard deviation of surface elevations, m j = ∫S o (f)f j df are spectral moments. Further, S(f, θ) is the directional wave spectrum with θ as the direction of waves at frequency f, and the cyclic frequency is ω = 2πf.
The dominant wave period T p = 2π/ω p refers to the cyclic frequency ω p of the spectral peak. The mean zero-crossing wave period T 0 is equal to 2π/ω 0 , with ω = m m / Here, overbars imply statistical averages and σ is the standard deviation of surface wave elevations. For second-order waves in deep water 10 Here, ν is the spectral bandwidth defined above and the characteristic wave steepness μ m = k m σ, where k m is the wavenumber corresponding to the mean spectral frequency ω m 8 . For narrowband (NB) waves, ν tends to zero and the associated skewness λ 3,NB = 3μ m [8][9][10] .
For third-order nonlinear random seas the excess kurtosis where λ 3,NB is the skewness of narrowband waves 8 . As for the dynamic component, Fedele 14 recently revisited Janssen's 62 weakly nonlinear formulation for λ d 40 . In deep water, this is given in terms of a six-fold integral that depends on the Benjamin-Feir index = BFI m v 2 / m and the parameter σ ν = θ R /2 2 2 , which is a dimensionless measure of the multidirectionality of dominant waves 11,15 . As waves become unidirectional (1D) waves R tends to zero and a random narrowband wave train becomes unstable if BFI > 1 65 .
The Tayfun-Fedele model. 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-Fedele model 9 , . In our studies Λ is approximated solely in terms of the excess kurtosis as Λ appr = 8λ 40 /3 by assuming the relations between cumulants 66 λ 22 = λ 40 /3 and λ 04 = λ 40 . These, to date, have been proven to hold for linear and second-order narrowband waves only 12 . For third-order nonlinear seas, our numerical studies indicate that Λ ≈ Λ appr within a 3% relative error in agreement with observations 67,68 .
For second-order seas, referred to as Tayfun sea states 34 R 2 Note that the Tayfun distribution represents an exact result for large second order wave crest heights and it depends solely on the steepness parameter defined as μ = λ 3 /3 10 26 is a numerical pseudo-spectral method, based on a perturbation expansion of the wave potential function up to a prescribed order of nonlinearities M in terms of a small parameter, the characteristic wave steepness. The method solves for nonlinear wave-wave interactions up to the specified order M of a number N of free waves (Fourier modes). The associated boundary value problem is solved by way of a pseudo-spectral technique, ensuring a computational cost which scales linearly with M 2 Nlog(N) 70,71 . As a result, high computational efficiency is guaranteed for simulations over large spatial domains. In our study we used the West formulation 26 , which accounts for all the nonlinear terms at a given order of the perturbation expansion. The details of the specific algorithm are given in Fucile 70 and Fedele et al. 2 . The wave field is resolved using 1024 × 1024 Fourier modes on a spatial area of 4000 m × 4000 m. Initial conditions for the wave potential and surface elevation are specified from the directional spectrum as an output of WAVEWATCH III 72 . Data Availability. All the publicly available data and information about the El Faro accident are posted on the National Transportation Safety Board (NTSB) website 1 .