Extreme environmental forcing on the container ship SS El Faro

ABSTRACT The sinking of the cargo ship SS El Faro is investigated by providing a comprehensive analysis of the wind, wave and ocean currents associated with Hurricane Joaquin. Using state-of-the-art reanalyses the event is assessed in high resolution and from a long-term climate perspective. The last known location of the SS El Faro was in the north-west eye-wall of Hurricane Joaquin when it was a category four major Hurricane. The maximum individual wave height in this region was over 10 m and the Benjamin-Feir index was 0.69 indicating a high likelihood of rogue waves. As the vessel tried to outrun the hurricane it was continually impacted by strong wind and waves on its port side. This was compounded with flooding that caused a starboard list which likely eventually caused the vessel to sink.


Introduction
The International Convention for the Safety of Life at Sea (SOLAS) ensures safety standards are set in the construction, equipment and operation of merchant ships. The risks posed to a vessel during a voyage include: shifting cargo; running aground; collision and encountering inclement weather. An example of how shifting cargo can affect ship stability occurred with the vessel M.V. Hui Long. Munro and Mohajerani (2016) describe how liquefaction of the wet granular cargo occurred which altered the stability of the vessel and ultimately caused it to sink. An example of a ship running aground is the Costa Concordia cruise ship which capsized after striking a submerged rock off the coast of Isola del Giglio, Italy in 2012. A recent collision event which got international attention was the United States Navy destroyer USS Fitzgerald and a Philippine-flagged container ship off the coast in Japan in Summer 2017.
It is estimated that about one third of major ship accidents are associated with inclement weather (Toffoli et al. 2005). There have been multiple cases of vessels being sunk or damaged by very extreme sea states as is given in Cardone et al. (2015); Heij and Knapp (2015). Extra-tropical cyclones in the North Atlantic often cause extreme waves (Breivik et al. 2014;Bell et al. 2017) because of a 'trapped fetch', where a weather system is moving at the same speed as the group velocity of the waves and therefore providing the waves with continual forcing for growth. This situation caused the sinking of the fishing vessel Andrea Gail which was made famous in the movie 'the perfect storm'. Extratropical cyclones pose a serious risk to vessels due to their size and localised strong winds (Gyakum 1983). Another factor which is a risk to container ships is parametric rolling in stormy seas (France et al. 2003). This is where large container vessels encounter head-on seas with a resonant wave period resulting in large roll angles that are coupled with significant pitch motions.
Hurricanes also pose an adverse risk to vessels as was the case with M. V. Derbyshire which was sunk by Typhoon Orchid in September 1980 (Faulkner 1998).
Hurricanes are known to generate extreme waves (Curcic et al. 2016) and the waves are often multi-directional and short-crested. Fedele et al. (2017) found that the possibility the SS El Faro encountered a rogue wave in hurricane Joaquin was high at 1/130 in a 50-minute window using Higher-Order Spectral simulations.
Ocean wave period plays a role in determining wave steepness and therefore the likelihood of a rogue wave occurring. This was the case in the North Sea in December 2012 where 19 rogue waves occurred during a storm (Gibson et al. 2014;Bell et al. 2017). Ocean currents also influence wave behaviour both on the large scale and for individuals waves (Ardhuin et al. 2017). Strong ocean currents in the Agulhas are known to modify the height and shape of ocean waves causing extreme sea states which pose a risk in a main shipping route (Quilfen et al. 2018). It is therefore important to take into consideration wind, wave and currents when understanding the forces acting on a ship (Bell and Kirtman 2018).
The mechanism by which a vessel sinks in rough weather conditions is determined by a combination of environmental forcing (mostly wind and waves) and ship hydrodynamics. Each ship, by design, responds differently to various wave groups (i.e. significant wave height, wave periods and directions). These are determined by response amplitude operators (RAO; Cattrell et al. 2018). Unfortunately, the RAOs are not available for the SS El Faro hence this manuscript focusses on the environmental conditions that the vessel experienced. The increased resistance to a ship's maneuverability from wind and waves also plays a part in the risk of vessel sinking during inclement weather. When a storm is approaching and the situation deteriorates the vessel can lose speed based on the environmental forcing, the hull, engine and propeller characteristics (Luo et al. 2016). This in turn limits the ability of the ship to outrun the storm.
The purpose of this study is to determine the evolution of the environmental conditions which led to the sinking of SS El Faro. There has only been one study pertaining to extreme waves during Hurricane Joaquin. Fedele et al. (2017) focused on the role of extreme waves in the sinking of the vessel but did not investigate the role of wind and currents. In our study, we use state-of-the-art reanalyses datasets to assess the wind, wave and currents associated with Hurricane Joaquin and the environmental forcing along the SS El Faro's final transit.
The layout of this paper is as follows. The SS El Faro event is explained in Section 2. The data and parameters are described in Section 3. Validation of the reanalyses and the environmental conditions of the event are given in Section 4. A discussion is provided in Section 5 and lastly the conclusions are given in Section 6.

The SS El Faro
The SS El Faro is a container ship which ferried supplies between the United States and Puerto Rico. During a voyage in September 2015 SS El Faro encountered Hurricane Joaquin ( Figure 1) and communications were lost. On the 31st of October the United States Navy Ship Apache identified the location of the sunken vessel. It was not until August 2016 that the Voyage Data Recorder (VDR) was recovered which contained information of the vessels position, date and time; engine information and audio from the bridge. This information was used in the National Transportation Safety Board's (NTSB) investigation which concluded in December 2017. The NTSB criticised the captain's decision to advance into the oncoming storm and noted he had relied on outdated weather information. It found flooding occurred in a cargo hold from an undetected open scuttle which caused a starboard list. This in turn affected oil levels in the engine room and the vessel lost propulsion (NTSB 2017).

SS El Faro automatic identification system
The SS El Faro ship data is provided by the NTSB's investigation (https://dms.ntsb.gov/pubdms/search/document. cfm?docID=447557&docketID=58116&mkey=92109). The Automatic Identification System (AIS) provides vessel position logs as well the speed over ground (SOG: speed relative to the surface of the earth), the course over ground (COG: direction of progress) and heading (direction the bow is pointing). The data starts at 2015-09-28T20:34:044 UTC when the vessel was in the port of Jacksonville, Florida. The last data point is at 2015-10-01T11:57:07 UTC to the Northeast of Acklins and Crooked Island, Bahamas. Figure 2 provides the SOG and COG of the SS El Faro along its transit. The drop in SOG is notable at 2015-10-01T09:00:00 UTC associated with flooding onboard. The final drop in the vessels speed from 8 to 5 m s −1 occurred at 2015-10-01T11:49:13 UTC, shortly before contact was lost with the vessel. Associated with the drop in the vessel's speed was a change in the vessels heading. As the ship lost power the vessel was being turned towards the east and towards the Hurricane.

IBTrACS
Hurricane Joaquin data are obtained from the International Best Track Archive for Climate Stewardship (IBTrACS; Knapp et al. 2010). This provides positional data of the Hurricane and wind strength at six-hourly time steps. IBTrACS provides the best observations of hurricanes using all available satellite and in-situ observations (Landsea and Franklin 2013). The data for Joaquin in IBTrACS originates from HURDAT (Hurricane Databases).

Global drifter program
Surface current speed measurements are obtained from drifter observations and the methodology is fully described in Laurindo et al. (2017). The data is obtained from undrogued and 15-m drogue Global Drifter Program (GDP) drifters. This dataset contains more than 29 million six-hour position and velocity estimates from February 1979 to June 2015. The data is provided at 0.25°spatial resolution. Laurindo et al. (2017) note that the ocean surface current climatology is well represented in this dataset compare to other climatologies derived from in-situ data, as well as when compared to altimeter-derived geostrophic velocity fields.   Table 1.

Globwave
The 10-m surface wind speed and significant wave height data are obtained from the Globwave project (Gavrikov et al. 2016;Queffeulou and Croizé-Fillon 2016) as representative of our best estimate of observed conditions. This dataset consists of a 23-year time-period  and uses data from nine altimeters: ERS-1&2, TOPEX-Poseidon, GEOSAT Follow-ON (GFO), Jason-1, Jason-2, ENVISAT, Cryosat and SARAL. Data is restricted to 2000-2015 to match the same period as the ERA5 data. An empirical function is applied to the backscatter (e.g. Sepulveda et al. 2015) to calculate surface wind speed and significant wave height. The data is corrected using co-located buoy measurements (Skandrani et al. 2004). In general, the altimeter data are close to observations with some biases of overestimating low wind speed and significant wave heights, as well as underestimating high wind speed and significant wave heights. The satellite track files have been interpolated to match the spatial resolution of ERA5. It should be noted that altimeter winds are assimilated into ERA5.

HYCOM
The HYbrid Coordinate Ocean Model (HYCOM) is an ocean model which is isopycnal in the open ocean, terrain following in shallow coastal regions and z-level coordinates in unstratified seas. Data is used from the latest version of the Global Ocean Forecasting System (GOFS) 3.1 which is the HYCOM model combined with the Navy Coupled Ocean Data Assimilation (NCODA) system (Cummings 2005;Cummings and Smedstad 2013). GOFS 3.1 has improved physics and a well resolved surface layer. The horizontal resolution is 0.08°in the tropics and the temporal frequency is three-hourly. The draft of SS El Faro was 12 m, therefore surface current data is integrated over the top 12 m. Studies have shown HYCOM generally captures most energetic and persistent features such as western boundary currents but is limited with some mesoscale features (Savage et al. 2015;Carvalho et al. 2019).

ERA5
3.6.1. IFS The European Centre for Medium-Range Weather Forecasts (ECMWF) 5th Generation (ERA5) reanalysis is used to assess the wind and waves. ERA5 is based on the predictions of a weather forecast model (Integrated Forecast System (IFS) Cycle 41r2) constrained by assimilation of data. The horizontal resolution of the atmospheric component of ERA5 is 31 km (T639) and data are available hourly (C3S 2017). Eventually, this reanalysis will replace ERA-Interim (Dee et al. 2011) thus validation of ERA5 data is still ongoing. Studies point towards that the higher resolution provides an improved assessment of extreme weather events (Olauson 2018).

WAM
The ocean wave model, the Wave Assimilation Model (WAM) is coupled with the atmospheric model and they communicate through the Charnock parameter which determines the roughness of the sea surface. The ocean wave model has slightly coarser resolution at 40 km. It uses 24 directions and 30 frequencies. ERA5 outputs the two-dimensional spectrum at every grid point as well as integrated parameters. The derivation of the integral parameters is provided below.
3.6.2.1. Ocean wave integral parameters. The twodimensional wave spectrum describes how wave energy is distributed as a function of frequency f and direction θ at each grid point: (1) The one-dimensional wave spectrum describes how wave energy is distributed in f and is defined as: Bulk parameters can be obtained from the spectrum using spectral moments with the nth momentum given as: The significant wave height H s is defined as: Mean wave period T m−1 is given as: Mean zero-crossing wave period T m2 is defined as: Peak wave period T p is defined as: where f p is the peak frequency of the one-dimensional spectrum.
Mean wave direction u m is given in oceanographic form with direction the waves are travelling toward and 0°being true north. It is defined as: Mean directional spread s u has values between zero (uni-directional spectra) and 2 √ (uniform spectra) and is defined as: where u f is the mean direction u at frequency f : 3.6.2.2. Rogue wave parameters. Rogue waves are defined as induvial waves that are much larger than surrounding waves. One common simple metric is if an individual wave is twice as large as the significant wave height H = 2H s . A spectra wave model cannot compute individual wavestherefore rogue wavesexplicitly but it can indicate an increased probability of their existence. One approach in which this is possible is by relating the shape of the probability density function of the surface elevation to the mean sea-state as described by the two-dimensional frequency spectrum (Janssen 2003). In a normal sea-state the distribution of the surface elevation has a gaussian shape. However, under exceptional circumstances deviations from Normality may occur which indicates the waves are nonlinear and there is an increased probability of rogue waves (Tayfun 1980;Tayfun and Fedele 2007). The deviations from Normality are measured in terms of the kurtosis of the sea surface elevation: with h as the surface elevation.
According to the theory of wave-wave interactions the kurtosis C 4 is related to the frequency spectrum by: where g is acceleration due to gravity, P is the Cauchy principal value, v is the angular frequency, T 1,2,3,4 is a complicated function of the four wave numbers (k n ; Krasitskii 1990). v 4 is: However, for an operational wave model, defining C 4 using Equation (12), a six dimensional integral, is too computationally expensive (Cattrell et al. 2018). Rogue waves most likely occur when the spectrum is narrow banded. One can therefore make a narrow band approximation. Based on numerical simulations of the Nonlinear Schröndinger equation which follows from the narrow-band limit of the Zakharov equation (Zakharov 1968) a fit can be obtained giving the maximum of the kurtosis. The total kurtosis can be approximated as: with d u the relative width of the directional spectrum. BFI is the Benjamin-Feir Index given by: with ε the integral steepness parameter: and d v the relative width of the frequency spectrum. For deep water α = 6. The skewness C 3 of the probability density function of the surface elevation is useful to quantify the contribution of bound waves to the deviations from Normality. It is given as: k 2 3 = 5(k 2 30 + k 2 03 ) + 9(k 2 21 + k 2 21 ) + 6(k 30 k 12 + k 03 k 21 ) Where k n refers to the third-order cumulants of the of the joint probability density function of the surface elevation and its Hilbert transform (Janssen 2014): k 30 = 31; k 03 = 0.; k 12 = 0.; k 21 = (k 30 /3) = 1. These parameters to describe the deviation from Normality can be used to come up with an expression for the expected maximum wave height H max . For a large number of independent wave groups N and small C 4 , H max can be approximated as: whereẑ 0 is related to the number of independent wave groups N (Janssen 2014): T D is the duration of the time series and γ is Euler's constant (0.5772). The expected wave period associated with H max (T max ) can be computed using the joint probability of the normalised envelope R and normalised period T: where v is the width parameter (Longuet-Higgins 1983): For a given normalised envelope height, wave period follows the conditional distribution of wave periods p(T|R) = p(R, T)/p(R): T max becomes the normalised period multiplied by the expected value of the period: Maximum steepness is defined as:

Model validation
It is important that any biases or limitations of the models used to investigate the environment conditions during the SS El Faro's final voyage are clearly stated.
In this section, we compare the spatial biases of the models climatologies with the best estimate of observed conditions. Figure 3 shows a comparison of the derived annual surface current climatology from HYCOM and that of insitu drifters. Despite HYCOM having a higher spatial resolution than the drifter climatology it still simulates a weaker Florida current than observed. This has been noted in other studies. For example Duerr et al. (2012) found HYCOM under predicts mass transport at 27°N when compared to ADCP (Acoustic Doppler Current Profiler) measurements. The Antilles current, to the North of the Bahamas, is stronger in the drifter data compared to HYCOM by ∼0.25 m s −1 . In general, along the SS El Faro's final transit the HYCOM simulation is remarkedly similar to the drifter data, especially east of 75°W.

ERA5 10-m wind speed validation
The Globwave annual 10-m wind speed climatology (Figure 4(b)) shows the altimeter tracks in the interpolated dataset and is demonstrable of a limitation of

JOURNAL OF OPERATIONAL OCEANOGRAPHY
using satellite altimeter data. However, the altimeter data can distinguish large-scale features such as stronger wind speed around 15°N, 73°W and weaker wind speed to the west of Haiti. The 10-m wind speed simulated by ERA5 compares well to Globwave with most of the biases within −0.5-0.5 m s −1 . ERA5 wind speeds are slightly weaker in generally compared to Globwave although altimeters generally overestimate low wind speeds. The 10m wind speed is stronger in ERA5 around the Caribbean Islands, but it is hard to validate altimeter data in coastal regions. In the region of the vessel route from Jacksonville to Puerto Rico ERA5 10-m wind speed is in good agreement with Globwave.

ERA5 significant wave height validation
ERA5 H s is in good agreement with Globwave ( Figure 5 (c)).  Figure 6 (a). The red cross denotes the centre of Hurricane Joaquin and the purple triangle is the position of SS El Faro. There is a Hurricane wind field signature in the ocean surface current response showing cyclonic motion. Figure 6(b) shows the same time slice, 2015-10-01T12:00:00 UTC, but as an anomaly compared to the climatological mean. Current speeds were twice as fast compared to climatology during Hurricane Joaquin especially on the south side of the storm. The ocean current speed and direction the vessel encountered was mostly consistent during the transit (Figure 7). The current speed was generally low between 0.2 and 0.5 m s −1 and the ocean current direction was unfavourable for the journey as the current was moving towards the west. Figure 8 shows SS El Faro was located in the front right quadrant of Hurricane Joaquin where the highest winds speed are often found (Kossin et al. 2007). The vessel is also located in the eye wall which is the region with the most intense wind speeds. in ERA5, the strongest 10-m wind speed on 2015-10-01T12:00:00 UTC was 31 m s −1 . The observed 10-m wind speed at this time was 59 m s −1 as given in IBTrACS. All reanalyses and to some extent operational global forecasts struggle to capture the strongest winds in hurricanes (Murakami 2014). It should be noted that the position of the hurricane given in ERA5 is in good agreement with observations. The time series of 10-m wind speed and wind direction at the location of SS El Faro is given in Figure 9. The SS El Faro experienced the largest winds speeds at approximately 2015-10-01T10:30:00 UTC, 90 min before contact was lost. The wind direction was constant at approximately north-easterly for about 18 h from 2015-09-30T18:00:00 UTC until 2015-10-01T10:30:00 UTC. This would have ensured the sea-state was well developed and the waves had a constant directional source of energy. The decrease in wind speed between 2015-09-30T11:00:00 UTC and 2015-09-30T12:00:00 UTC is associated with SS El Faro entering the eye of hurricane Joaquin. Figure 10(a) presents the maximum individual wave height H max on 2015-10-01T12:00:00 UTC and Figure  10(b) shows how anomalous this period is compared to climatology. SS El Faro is located close the region of maximum H max (the north-east quadrant) and is in a region where H max is greater than 10 m. The maximum H max is located to the north-east of the storm centre and the waves are smaller to the south-west as the islands dissipate the waves energy.

Ocean waves
The mean wave direction u m the SS El Faro experienced shifted southward along the voyage as the Hurricane passed (Figure 11(a)). The wave period associated with the maximum individual wave height T max got shorter as the Hurricane approached SS El Faro from 9.3 to 7.2 s (Figure 11(b)) and as a result the wave steepness increased.

Wave spectra and rogue waves
Wave energy increased as the SS El Faro continued its journey south-eastward. This can be seen in the onedimensional wave spectra (Figure 12). The blue line gives the wave spectra integrated over all directions at the location of SS El Faro at 2015-10-01T03 UTC. At this time the wave field was almost fully developed with the peak wave energy at 0.1 s radian −1 . Throughout the next nine hours the wave energy increased and the spectra became more narrow banded. The peak of wave energy slightly shifted to shorter waves as the vessel The two-dimensional wave spectra for the same location and time are displayed in Figure 13. The contours represent wave energy on a logarithmic scale. At 2015-10-01T03:00:00 UTC (Figure 13(b)) the wave energy is mostly travelling to the west and south-west. As the Hurricane neared, the wave energy backed and moved toward a southward direction. At the last known location, the waves became multi-directional and a secondary partition is apparent with some waves travelling north-westward. These can be seen in the directional spreading (s u ) which increased from 0.56     Table 1. to 0.91 between 2015-10-01T09:00:00 and 2015-10-01T12:00:00 (Table 1).
Along with the wave energy increasing there was also an increased likelihood that rogue waves occurred. The wave spectral kurtosis (C 4 ), Benjamin-Feir Index (BFI) and wave spectral skewness (C 3 ) all increase during this period. The BFI, which represents the ratio between wave steepness and spectral bandwidth, is large at 0.69. As a comparison the BFI was 0.24 during the Andrea rogue wave event (Fedele et al. 2016;Donelan and Magnusson 2017). The C 4 and C 3 values increase but only slightly over the period of nine hours. It can also be seen that H max /H s remains constant during this period at ∼1.9.

Discussion and conclusions
The evolution of the environmental forcing on the container ship SS El Faro during the final nine hours is presented in Figure 14. The transit from Jacksonville to San Juan has a COG of approximately 145°. The vessel experienced strong persistent winds on her port side as the captain made the decision to continue the journey and attempt to outrun the hurricane. This was compounded by the fact the hurricane had an unusually south-westerly track so the situation deteriorated over time. The red arrow is the 10-m wind speed direction which is almost perpendicular to the vessel throughout the time-period shown. This explains the starboard list of the vessel. The blue arrow is the mean wave direction (u m ) and indicates that it mostly follows the wind except for 2015-10-01T09:00:00 where u m is southward and the 10-m wind direction is easterly. When the vessel lost power shortly before 2015-10-01T12:00:00 the      Figure 12. The grey arrow indicates the course over ground; the green arrow denotes the ocean surface current direction. The red arrow shows the direction of the 10-m wind and the blue arrow denotes the mean wave direction. environmental conditions caused the vessel to turn to the port side thus the angle of the wind to the COG was ∼90°. The 10-m wind speed, u m and current direction were all impacting the vessel's port side. This is an extremely dangerous situation as the vessel's roll is largest when wind and waves are perpendicular to the vessel's heading. Once the vessel lost power it was unable to correct its position relative to the wind and waves. This meant there was a greater chance that the vessel was struck by a rouge wave. This hypothesis is strengthened when looking at parameters that represent the likelihood of rouge waves. The individual maximum wave height and the Benjamin-Feir Index were largest at the point in time and space of the last known location of the vessel.
As ERA5 is a spectra wave model it cannot provide information of individual wave heights as were provided in Fedele et al. (2017). Thus, it cannot be used to explicitly quantify the likelihood that a rogue wave sunk the SS El Faro. However, the study of Fedele et al. (2017) can be expanded upon using initial wave field conditions of ERA5 for Higher-Order pseudo-Spectral (HOS) simulations.
The models used in this study provide good agreement with observed datasets to represent the wind, wave and ocean current climatology of the Caribbean. GOFS 3.1 shows similar results with observed drifter data although HYCOM simulates the Florida current and Antilles current weaker than observed. This indicates that while HYCOM is a relative high-resolution ocean model with 1/12°grid spacing, higher resolution or improved physical parametrizations may be needed to accurately simulate the strength of Ocean boundary currents. The meanstate 10-m wind speed and significant wave height in ERA5 compared to altimeter measurements are mostly in good agreement. The difference is largest around the small Caribbean islands. This is not surprising as ERA5 is a global model and uses sub-grid parametrization to represent wave energy dissipations by islands smaller that a grid cell. It is an improvement over ERA-Interim but more research is needed in this space.
ERA5 simulated the maximum individual wave height H max to be greater than 10 m at the last recorded location of SS El Faro which was in the north-west eye wall of Hurricane Joaquin. The wave period associated with the maximum individual wave height T max decreased along the vessel's voyage and the steepness increased. The coefficients of skewness and kurtosis increased and peaked after the vessel transmitted its last location. In addition, the Benjamin-Feir index increased to 0.69 indicating there was a high likelihood that rogue waves occurred.
A comprehensive analysis of the ocean current direction, 10-m wind direction and two-dimensional ocean wave spectrum respective to SS El Faro's course over ground is portrayed. This provides qualitative information of the dynamic behaviour of the ship without the existence of the vessels RAOs (Ibrahim and Grace 2010). The strong persistent wind and waves on the vessels port side amplified the starboard list. This coupled with a high probability that rogue waves were present are what were likely to have caused the vessel to sink. An improved knowledge of ocean conditions caused by extreme weather can improve safety during vessel routing (Edwing 2018). This study demonstrates that HYCOM and ERA5 can be used further to understand the relationship between extreme environmental forcing and vessel safety.