An Operational Numerical System for Oil Stranding Risk Assessment in a High-Density Vessel Traffic Area

In the Mediterranean Sea unique environmental characteristics and sensitive assets coexist with intense maritime traffic that is represented by frequent daily passages of vessels along the main waterways. In order to assess the risk of oil stranding in case of at-sea emergencies and provide key products for environmental agencies or policymakers preparedness, a geographically relocatable, operational numerical system is implemented and tested. The system relies on the application of oceanographic and particle tracking models and is able to provide, on a high-resolution and unstructured computational grid, a 3-days forecast of those variables known as the main drivers of oil slicks at sea. The risk of potential oil stranding is computed through a combination of anthropogenic hazard and shoreline vulnerability. The sources of hazard vary on time and space in relation to local maritime vessel traffic. The shoreline vulnerability is based on the current knowledge of slope, main grain size, geology of rocks, and occurrence of manmade structures at coast. The operational system is enriched by a web graphical user interface and includes automatic and on-demand working modes. Its functionality is demonstrated in the Strait of Bonifacio (western Mediterranean Sea), area with a high potential risk of oil stranding due to an intense maritime traffic. Risk assessment is hence computed for a test year, the 2018. Critical values of risk are found in correspondence of long stretches of littoral while many of them are currently characterized by a low anthropogenic pressure. The results emphasize the geomorphological features of the shorelines as reducing or amplifying factors to any potential impact of oil stranding at coast.


INTRODUCTION
At the end of 2019 the European Maritime Safety Agency has published traffic density maps (European Maritime Safety Agency, 2019) which provide an astonishing description of ship movement patterns to both professionals and the public. These products raise citizen's awareness and encourage science and authorities to implement operational risk management systems. Their objective is to improve marine traffic planning to enhance ship safety and adequately protect the coastal environment from damaging consequences of potential pollution events.
The Mediterranean Sea is a major shipping routes with about 30% of the total world vessel traffic (Campana et al., 2017). Its space hosts important transit routes for international shipping activities. Thousands of cargos, oil tankers and all kind of vessels populate the sea routes and overcrowd sea straits like the busy Gibraltar and Bosporus straits, the Strait of Sicily and the straits of Messina and Bonifacio.
The Mediterranean Quality Status Report (2017) of the United Nations Environmental Programme cites a correlation between vessel traffic density and accidents causing pollution, being collisions, the first cause of accidents resulting in an oil spill. REMPEC (2008) estimated that 18%, or 421 million tons, of global seaborne crude oil shipments took place through the Mediterranean Sea in 2006. Thus, the routes of oil tankers are one of the main sources of danger that exacerbates the risk of an oil spill accident or a deliberate illegal discharge.
Moreover, maritime traffic may also affect the high biodiversity of the Mediterranean Sea being one of the major threats for marine environment in terms of collisions with cetaceans (Pennino et al., 2016;Ritter and Panigada, 2019), noise, alien species, and pollution (Abdulla and Linden, 2008;Coll et al., 2012). In particular, chemical pollution by oil or other toxic substances severely impact on water, fauna, flora and seabed sediments, and such an impact becomes more manifest when trails of pollution reach the littorals. There, except for relatively limited industrial areas, many social and economic activities (tourism, fisheries, and spin-off companies) often rely on a good environmental status of the ecosystem. In Italy, national parks like Cinque Terre, the Tuscany (northern Tyrrhenian Sea), and La Maddalena (Sardinia) archipelagos and the Sardinian Asinara Island embrace extended shorelines, representing, along with the marine protected areas, a high environmental, social and hence economic value for the country.
In response to international regulations and recommendations (see United Nations Convention on the Law of the Sea-UNCLOS) national and local authorities have currently available numerous operational risk management systems that can provide pollutants hazard assessment scenarios or a rapid hazard evaluation in case of unexpected at-sea accidents. These systems are mostly based on the application of the particle-tracking models (PTMs) to predict the trajectory of pollutants dispersion under specific meteorological and marine conditions. Over years, commercial companies or research institutions that are often involved to provide risk assessments in case of accident, have developed PTMs based on structured computational grid (e.g., Caballero et al., 2008;Sotillo et al., 2008;Coppini et al., 2011;Liubartseva et al., 2016;Zodiatis et al., 2016;Sorgente et al., 2020) or based on unstructured computational grid like those presented by Chen et al. (2007); Azevedo et al. (2009), Wang and Shen (2010), Olita et al. (2012), and . The former approaches offer successful results in the open sea, defining shape of the slick and its volume changes in relation to weathering processes, the latter find appropriate application into irregular sea domains with a complex bathymetry and provide, in addition, efficient and detailed indications on pollutants displacement.
The PTMs based on unstructured computational grids allow high spatial resolution for those parts of the domain where it is required, and they can provide hazard information at the same spatial scale of the littoral's vulnerability information. Adopting this type of tool, the authors aim to realize an operational numerical system to support the oil stranding risk preparedness and predict the risk that unexpected oil spills could impact the littorals of the high-density vessel traffic area of the Strait of Bonifacio and its surroundings (see www.seaforecast.cnr.it/sicomarplus/index.php/en). The preparedness can actually reduce the risk that many sensitive littorals could suffer of a potential oil stranding with consequent long natural recovery and it relies on the scientific knowledge of local environmental forcing and coastal vulnerability. The latter is defined on the base of the main geomorphological features of the littorals that will modulate the risk computation in the domain of investigation revealing unexpected critical states along pristine traits of coast.
The numerical system is indeed designed to combine the current knowledge of the local environment features with technological innovations, improving decision making. This is possible through the production of risk assessment scenarios that, for the first time in this region, include the evaluation of vulnerability in correspondence of the littoral and the adoption of a dynamical source of hazard, represented by the spatial and temporal variations of the vessel traffic information.
The numerical system consists of a coupled suite of oceanographic and PTM models, a vulnerability database and a risk computation engine. The essentials of the system have been applied in previous studies of waves, hydrodynamics and pollutant surface transport in shallow waters and coastal seas (e.g., Canu et al., 2015). Ad hoc designed web graphical user interface is also implemented to easily manage the system (Figure 1). Hazard is estimated by the PTM solutions that give information on pollutants that can potentially reach the coast within a specific time frame. The potential sources of hazard vary on space and time and they are based on the vessels traffic records 1 .
Vulnerability information are derived by the intrinsic properties of the land-sea interface like slope, geology of rocks, and grain size that can reduce or amplify the littoral vulnerability in case of pollutants stranding (Gundlach et al., 1978;Wang et al., 2005;Olita et al., 2012;Grottoli and Ciavola, 2019).
The risk is computed in purely quantitative terms, resulting from the linear interaction between anthropogenic hazard and the in-situ vulnerability (e.g., Christensen et al., 2003; United Nations/Inter-Agency Secretariat of the International Strategy for Disaster Reduction, 2004;Azevedo et al., 2009 for vocabulary;De Lange et al., 2010).
The functioning of the system is demonstrated through the computation of an oil stranding risk assessment scenario in the coastal domain of the Strait of Bonifacio (SoB domain, FIGURE 1 | Flow chart representing the functioning scheme of the numerical system for oil stranding risk assessment. hereafter), a busy strait to navigate between the north-western Tyrrhenian Sea at east and the Algero-Provençal basin (western Mediterranean Sea) at west. All products are computed for the test year 2018.
The main geophysical features of the SoB and its surrounding are described in section "Site Description." All components and the functioning of the operational numerical system along as the innovations in hazard, vulnerability and risk metric are described in section "Materials and Methods." The main outcomes of the system like the sea conditions operational products, vulnerability, hazard and risk assessment are explained in section "Results." In section "Discussion, " a synthesis of the adopted approach and the results are discussed in relation to the recent literature to draw the main conclusions (section "Conclusion").

SITE DESCRIPTION
The domain of the SoB is part of the Pelagos Sanctuary 2 and it hosts a French and an Italian marine national park. The strait is connecting the Tyrrhenian Sea and the western Mediterranean Sea (Figures 2A,B) with a narrow passage in the east, due to shallow waters around Lavezzi (Corsica) and La Maddalena (Sardinia) archipelagos and a wider passage in the west that is opening to the Gulf of Asinara. The average depth of the strait is approximately 70 m with eastern and western openings slowly increasing at about 100 m before approaching the continental slope, where strong bathymetric gradients occur .
The winds are blowing around 10 months to a year and the north-western Mistral and the south-western Libeccio, primary wind regimes, can often exceed the maximum speed of 16 m/s .
The river inflows are characterized by a torrential regime determining a weak variation of salinity that is ranging in between 37.7 and 38.8 PSU throughout the year. Sea surface waters temperatures approximatively vary from 15 • C 2 www.sanctuaire-pelagos.org in early spring, to 25 • C in late autumn (Artale et al., 1994;. The tidal forcing is weak ) and the sea current circulation is mainly wind driven with observed maxima near 1.46 m/s (Gerigny et al., 2011). The pressure gradients between the adjacent basins (Artale et al., 1994) and its interaction with local bathymetry determine the deep waters dynamics.
Wave events may be severe at the entrances of the SoB. The analysis of a 10 years (2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) dataset, derived by a wave model hindcast (Copernicus Marine Service, 2020) gives significant wave height maxima of 6.4 m at the east side (41.36 • N, 9.61 • E) and 7.54 m at the west (41.40 • N, 9.0 • E) side of the SoB. In the first case the mean wave direction is from the north, in the second case from the north-west (Life -European Project, 2020, Life17CCM/IT/000121).
In the narrow ship transit corridor of the SoB the navigation may hence be dangerous, representing a relevant threat for its sensitive space. During the period 2000-2009, the Italian national system of maritime traffic control has recorded 35,188 transits through the SoB with numerous transits of high bunker volumes vessels, which are associated with a medium and high potential risk index .
Although in the course of such a decade vessel transits started a sensibly reduction

The Core of the Operational Numerical System
The core of the system is an evolution of previous ones that have been developed in the context of the European and Italian research projects SOS-Piattaforme e Impatti Offshore (Ribotti et al., 2019;Sorgente et al., 2020) and the Mediterranean Decision Support System for Marine Safety (Zodiatis et al., 2016). It is based on SHYFEM (Umgiesser et al., 2004), a suite of high-resolution numerical finite elements models describing waves, sea current dynamics and the fate of advected pollutants at the sea surface Canu et al., 2015).
The numerical tools are implemented over a shared unstructured mesh (Figure 2A) with a spatial resolution ranging in between 20 km offshore and 100 m in correspondence of specific focuses in the SoB domain (Figures 2B-E). Indeed, size and shape of the elements is variable, and the element density is higher along the coasts and in correspondence of specific sites like small islands. The vertical dimension is discretized by 30 uneven z-levels and the depth of each layer is varying between 2 m, for the surface, and 200 m, for the bottom. The mesh is constituted by about 93,000 triangular elements with a total of about 50,000 computational nodes. The hydrodynamic 3D model resolves the shallow water equations integrated over each layer in their formulations with water levels and transports. The equation system reads as: where l indicates the vertical layer, (U l , V l ) the horizontal transport components in x and y directions for each layer, Adv x l and Adv y l the advective terms, f the Coriolis parameter, p a the atmospheric pressure, g the gravitational constant, ζ the sea level, ρ 0 the typical water density, ρ = ρ + ρ 0 the water density, τ the internal stress term at the top and bottom of each layer, h l the layer thickness, H l the depth of the bottom of the layer l, A H the horizontal eddy viscosity (Smagorinsky, 1993) and, F x l , F y l the gradients of wave-induced surface stresses in the x and y directions, respectively. The General Ocean Turbulence Model (Burchard and Petersen, 1999) is used for the computation of the vertical viscosity K l . Wind and bottom friction terms, corresponding to the boundary conditions of the stress terms (τ x , τ y ), are defined as: with c D as the wind drag coefficient, c B the bottom friction coefficient, ρ a the air density (w x , w y ) the wind velocity components and (u L , v L ) the bottom velocity components.
Two boundaries are open to the basin scale dynamics and they are located at the western and eastern side of the mesh in the Algero-Provencal basin and in the Tyrrhenian Sea, respectively. There, the large-scale sea conditions are defined via Dirichlet formulation for water levels, whereas a nudging procedure was adopted for the temperature and salinity variables . At the closed boundaries, only the normal velocity was set to zero, whereas the tangential velocity was a free parameter corresponding to a full slip condition (Umgiesser et al., 2004).
The wave module WWM, solving the random phase spectral action density balance equation for wavenumber-direction spectra (Roland et al., 2009), is used to describe wave patterns. The module describes growth, decay, advection, refraction and diffraction of wind waves due to the wind action, depth gradient and coastal geometry and it provides the gradients of waveinduced surface stress to the hydrodynamic module.
The coupled hydro-wave numerical model provides the oceanographic drivers to the PTM that solves the advection and diffusion equations (Eq. 2) in a Lagrangian framework of reference to predict the particles trajectories Quattrocchi et al., 2016). The PTM equations system reads as: where u adv , v adv are the advective velocities components and u dif , v dif are the diffusive velocities components in x and y directions, respectively.
The components u adv and v adv are solved following: where u 1 , v 1 are the horizontal velocities computed by the hydrodynamic model for the 1st vertical layer; w x , w y are the wind velocity components; C c and C w are the current transport and the wind transport factors. The horizontal components of the turbulent diffusive velocities, u dif and v dif , are computed using a random walk technique on the basis of the Fisher et al. (1979) study, with turbulent diffusion coefficients obtained using the Smagorinsky formula (Smagorinsky, 1993).
The transport of oil particles is mainly due to the direct action of the currents and winds (Caballero et al., 2008;Sotillo et al., 2008;Liubartseva et al., 2016), whereas the wave contribution is generally related to emulsification and weathering processes that were considered but not explicitly modeled in this version of the operational numerical system (see Supplementary Material).

Model Calibration and Simulation Setup
The simulation setup of the SHYFEM model for calibration and further runs of the operational numerical system was designed for the case study application of the SoB domain (Figures 2A-E) and the necessary inputs are described in the following.
A calibration procedure was carried out to improves the accuracy of the SHYFEM models in predicting the sea surface transport when different wind regimes take place in the SoB domain. A Trial and Error method was used to reduce the discrepancies between Lagrangian currents measurements, collected during October 2018 and September 2019, and simulated trajectories.
According to statistical procedures found in , 225 numerical particles were released at each observed drifter position and their trajectories computed. At each run of the model calibration procedure, the values of C D , in Eq. 2, C w and C c , in Eq. 4, were changed within their typical ranges. For each run of the model, the obtained results were compared with observations and a trajectories relative error (TRE, see Cucco et al., 2016bCucco et al., , 2019. where, t is the prediction time lag, x and y the components of the drifter position at the time t, o and m indicating observations and numerical particles and D o indicating the length of the path followed by the drifters at the time t. The average of the TRE computed within the first 6 h from the release (TRE6) was selected as reference metric for estimating the model accuracy.
The calibration procedure brought to an optimal combination of model parameters consisting in to C D = 2.8 × 10 −3 , C c = 0.95, and C w = 4 × 10 −3 .
In Figure 3, all panels, the observed (black lines) and computed trajectories (blue and gray lines) are reported for the whole set of Lagrangian measurements. The model was able to reproduce the observations with a TRE24 ranging between 0.12 and 0.34. The discrepancies between numerical particles and drifter's location over the 24 h predictions varied on average between 1.3 km up to around 7.7 km with a best match for shortterm experiments A, B, C and E. Long-term experiments, like G and H showed a less accurate but satisfactory agreement between the modeled and observed trajectories, reproducing multiple days trajectories and, in the case G, the drifter landing area.
Atmospheric and ocean analysis and forecast are provided at the top and side open boundaries of the models with a 3-and 24-h frequency, respectively. Atmospheric fluxes are obtained by the European Centre for the Medium-Range Weather Forecasts (ECMWF), the state of art weather forecast system in Europe. Geomorphologic data, needed to setup the vulnerability index, were obtained by the Regione Sardegna Geodatabase (2020) 4 , the European Environmental Agency (2005) and with the adoption of the Google Earth software. In order to discretize the available information, squared cells were centered along the shoreline. The choice of a 100 m facet provides an effective description of the main geomorphological features of the coastlines in the SoB domain. Every cell contains relevant information about the enclosed tract of shoreline: the slope of the beaches, the main grain size of their sediments (gravel, sand, or mud), the geology of the rocks (hard rock, eroding rock) and the existence of manmade structures (bank, ports, jetties).
Sources of anthropogenic hazard were derived, for the year 2018, by monthly values of the vessel density, a measure of the shipping traffic in relation to the time and the unit of surface (see text footnote 1). Vessel density is released by means of a georeferenced grid composed by squared cells of 1 km per side and indicating the total time, expressed in hours, per surface unit, expressed in km 2 , during which the cells are occupied by ships.
For the computation of the risk assessment the sources of hazard and their space-time variation are represented by thousand numerical particles that are released where vessel density is defined for tanker, cargo, and passenger vessels. This choice is justified by the significant threat that such type of vessels represents for the SoB domain, considering their dimensions, dangerous load, and number in Mediterranean waters.
The top picture in Panel A of Figure 4 displays the monthly extent of those areas that, within the SoB domain, are crossed by such type of vessels, highlighting the macroscopic difference between tanker vessels and other categories. The bottom picture in the same panel displays vessel density, per month and for each vessel category showing a weak variability for tanker and cargo vessels all year long and an increase for passenger category during the best weather season.
In order to provide a geographical representation of the areas crossed by each category of vessels, Figure 4, panels B, C and D, gives annual mean values with the related vessel density. Tanker vessels typically follow a wide route at east, cross the SoB and, at west, they turn around the Asinara Island to reach or leave the industrial harbor of Porto Torres in Sardinia and Corsica, where they probably shelter in the bay. Cargo vessels traffic covers a wide area with high vessel density in correspondence of the SoB and toward the main commercial ports. Passenger vessels, including Ro-Ro, route toward the ports of the main towns where touristic interests are also found.

Hazard, Vulnerability and Risk: The Adopted Metric
According to Cutter (1993) the authors consider hazard as the probability of occurrence that a natural or anthropogenic event will affect a specific area within a given time period. The solutions of the hydro-wave model and PTM are used to compute the hazard statistic in correspondence of each discretized coastal segment.
On the basis of the monthly information about the extent of the areas and related vessel density values, numerical particles, representing the dangerous sources are released and they are advected by the PTM in relation to the hydrodynamic and wave conditions. Indeed, for each category of vessel, with daily frequency, a total of 600,000 particles are released per month within a variable area that is representing the density of the shipping routes.
The number of particles is consistent with the maximum spatial resolution of the hydrodynamic model to guarantee that numerical particles are released from all elements of the model mesh. Each particle location was updated computing the transport induced by the currents and the wind with a 300 s time step for a duration of 10 days after the release or until its stranding at coast. The release frequency is the same of the open ocean boundaries update while estimates of weathering processes time scale suggest that evaporation, dispersion, spreading, emulsification, and dissolution acting on the oil volume with immediate effects and within a few days (Coppini et al., 2011;Fernandes, 2018).
During the advection, numerical particles retain memory of the release time and the vessel density at the release location. Therefore, if one or a set of particles strand at coast, hazard is computed for each coastal segment by combining the following indicators: (i) impact density [m −1 ], the number of stranded particles per coastal segment length expressed as: where L is the length of the j th coastal segment and n is the total number of particles that reach the j th coastal segment.
(ii) averaged vessel density [h km −2 ], the average of the vessel density related to the release point of each stranded numerical particle where VD, the vessel density, is the number of hours per month that any vessel spends over a squared area of 1-km side and n is the total number of particles that reach the j th coastal segment.
(iii) residual volume of light (VLO) and heavy oil (VHO) associated to the stranded set of particles where for the j th coastal segment the residual volume over time (t) was obtained by evaluating the average age of such set of particles through specific theoretical volume reduction functions g and f.
Indeed, the values of the residual volume of oil associated to a set of stranding particles is computed on the basis of their average age (unit of measure: hours), the average of their ages since release The functions g and f were defined by fitting the solutions of a process study experiment that was conducted with MEDSLIK_II (Zodiatis et al., 2008;De Dominicis et al., 2013b, Part 1 and Part 2), a community model that is designed to predict the transport and weathering of an oil spill. In fact, MEDSLIK_II was here used to evaluate the effects of the weathering processes on a unitary oil volume under synthetic sea state conditions. The reader can refer to Supplementary Materials for a detailed description of the model simulation setup and related interesting solutions. It is here essential to know that function g describes how for a medium density Arabian crude oil, a type of oil carried by tankers, the unitary volume is reduced by the 55% after 10 days and 45% after 2 days since its release at sea, while function f describes how, for a highly viscous oil typically used as fuel for vessels, the unitary volume is reduced by 72% after 10 days and 52% after 2 days since its release at sea. Hazard is hence computed for the j th coastal segment and for each type of vessel according to the following equation where VO s j = k s * VLO s j + m s * VHO s j with the upper case s = {T, C, P} indicating T for tanker, C for cargo and P for passenger vessels, k and m used as weights to consider the type of oil transported by the vessels: k = m = 0.5 if one considers the volume of oil carried by tanker (crude oil plus fuel), k = 1 and m = 0 if we consider the volume of oil carried by cargo or passenger vessels (fuel only). Vulnerability can be defined as a measure of the weakness of the environmental system that can suffer one or multiple threats (Azevedo et al., 2017 and its reference list). Vulnerability of complex contexts, involving relationship between people and their ecological space, should primarily be related to the intrinsic properties of the environment that can represent weakness or strong points to potential hazard sources, that is, resistance and resilience of the exposed environment modulate the measure of its vulnerability (Salter, 1997;Kleissen et al., 2007;Lahr and Kooistra, 2010).
In the specific case of the intertidal environment, intrinsic properties include shoreline geomorphology that should be considered to evaluate the vulnerability of littorals, where human activities mainly insist and the exposure to pollution events is very high. For these reasons and in the context of the present case study an environmental sensitivity index (ESI rank = {1, 9}) was defined on the basis of the local geomorphologic features (i.e., slope, main grain size, geology, and manmade structures of the littoral), that are described per each 100 m squared cell, centered at the shoreline. The higher is the ESI rank the higher is the vulnerability of the shoreline.
Indeed, in order to assess the vulnerability of specific shorelines to oil stranding, local geomorphologic features were evaluated in relation to the shoreline classification that is found in ESI guidelines, ad hoc established by the National Oceanic and Atmospheric Administration to draw maps of vulnerability (NOAA, 2002).
Similarly as done by NOAA, the authors assigns ESI rank = {1, 2} to high-energy shorelines, regularly exposed to large waves or strong tidal currents during all seasons, ESI rank = {3, 4, 5, 6} to Medium-energy shorelines that often have seasonal patterns in storm frequency and wave size, ESI rank = {7, 8} to low-energy shorelines that are sheltered from wave and tidal energy, except during unusual or infrequent events, ESI rank = {9} to anthropized shorelines. A rapid (from days to weeks) natural removal of stranded oil is associated to those littorals with low values of ESI rank, a slow natural removal (years) is associated to those littorals with high ESI rank values, days to months is the time-scale of oil removal associated to intermediate ESI rank values (NOAA, 2002).
Embracing the expression reported by Azevedo et al. (2017), the risk concept can be defined as the expected losses due to the occurrence of the adverse event for a given time period and specified area, that result from the interaction between the natural or anthropogenic hazard and the vulnerable assets (United Nations/Inter-Agency Secretariat of the International Strategy for Disaster Reduction, 2004;De Lange et al., 2010). The numerical computation of the risk (RI) was hence implemented, along the shoreline and according to the following equation, as the product of the hazard and the ESI indices: where the risk is computed for j th coastal segment and for each type of vessels s = {T, C, P}.
The approximation in Eq. 13 excludes the values of other exposed elements, such as human lives, sensitive infrastructure and cultural heritage but it provides a straight and reasonable assessment of the risk (Blong, 1996;Sayers et al., 2002), leaving decision makers free to adopt proper management solutions in relation to a wider range of vulnerable elements.
The computed risk can change on the base of the adopted definition of hazard sources and/or vulnerability. Thus, for the sake of a clear representation and further comparisons with other risk assessment, the risk values are normalized by their maximum value, providing a scale that is ranging in between 0 and 1.

RESULTS
In the following the authors explain the main operational products of the prediction system and the results obtained by the modeling scenarios of risk assessment carried out for the high-density vessel traffic area of the SoB domain, during the test year 2018.

Sea Conditions Operational Products
The operational hydrodynamic and wave models provide, with daily frequency, sea conditions information for a 3-days forecast horizon.   Figure 5, where snapshots of the sea state condition products are reported for the entire domain of investigation (Panel A) and for the focus sub-areas of the Bonifacio Strait (Panel B) and Asinara Island (Panel C), that are characterized by a high naturalistic relevance. These pictures represent the sea state conditions that typically arises from the most common meteorological forcing in correspondence of the SoB domain , a Mistral wind that blows with high intensity, up to 12 m/s, for around 2 days.
The Panel A shows the pattern of the sea surface current. At west, far from the entrance of the SoB, the sea current dynamics is dominated by basin scale general circulation structures that start to interact with the morphology of the local coastal boundaries. In the coastal area, near the SoB, the current come from northwest, it is strong and enters the SoB following the same direction of the wind. Entering the SoB, current maxima of 50 cm/s result from local bathymetric and coastal morphology constraints. Once to the east the current follows the eastern coast of Sardinia meandering around the numerous islands and join the largescale circulation that, at the surface, also displays the effects of the strong northwesterly wind. Similar oceanographic conditions were observed by Gerigny et al. (2011) that, through the analysis of meteorological and oceanographic measurements in the SoB, explained how in many cases the wind and current are strong and orientated in the same direction. Also, their correlation is strong on the whole water column in case of strong winds (>7.0 m s −1 ). Panel B shows an example of significant wave height and direction of the wave field which highlight how complex coastal geometry can protect specific traits of coast, for instance, by reducing the strength of the wave and the related current field. Panel C shows an example of the sea surface temperature displaying a relatively homogeneous field characterized by a high temperature located at north of Asinara Island and a coastal high and low temperature dipole, within the Gulf of Asinara, that is modulated by local upwelling and downwelling due to a combined wind-coast orientation effect.
As part of the output, the Stranding Time is also reported in the Panel D of Figure 5. It consists on the time needed by numerical particles that are released at the sea surface to reach the littoral. The spatial distribution of the stranding time is computed each 6 h for the 1 st day of prediction using the hydrodynamic model and accounting for the transport processes induced by wind and currents of the next 48 h. Specifically, every day four simulations are carried out initializing the model domain with a homogeneous distribution of numerical particles that are advected by wind and currents for 48 h. The ages of the particles at the time of stranding is therefore computed and such vales referred to the particles release positions, obtaining, for each single run, the spatial distribution of the so-called Stranding Time.
This quantity reverses the paradigm of the risk from pollutant at sea identifying potentially endangered waters instead of coastal zones. The stranding time distribution, in fact, allows the detection of the areas where a hypothetical oil-spill could reach the coast quickly, for stranding time lower then few hours, or after longer period. In the first case, those areas can be defined as critical and more attention should be given to the control the vessel passages.

Vulnerability Assessment
The sensitivity of the littorals to potential oil spill stranding was evaluated in the SoB domain. As stated above, the coastline slope, main grain size, geology and anthropogenic infrastructures, contribute to the definition of the littoral's sensitivity that is showed in the map of Figure 6 by means of the assigned ESI rank.
There, high-energy littorals with high slope gradients, typically rocky cliffs, have low sensitivity with ESI rank = {1, 2}. This is due to their intrinsic features that could promote a relatively fast natural recovery of the site after potential oil stranding. In the SoB domain, the 75% of the littorals, including the archipelagos, have such vulnerability.
Mixed-energy littorals with medium slope gradients, such as gravel and boulder beaches, as well as sandy beaches, have medium sensitivity with ESI rank = {3, 6} and they represent the 22% of the SoB domain. Some examples are given by open beaches in the surrounding of Porto Torres and by pocket beaches of the Corsica and east Sardinia, where complex coastal geometries can protect from severe sea conditions but, at the same time, preclude from any fast, natural recovery of the site after a potential oil stranding.
Low-energy and low slope littorals, like open and/or sandy beaches, and littorals with anthropogenic pressure, have high vulnerability with ESI rank = {7, 9}. They are located along the littorals of Porto Torres, Olbia, La Maddalena (Sardinia), Bonifacio and Porto Vecchio (Corsica), representing a percentage of the 3% in the SoB domain.
The vulnerability that is defined by ESI rank will be used in linear combination with hazard values to assess the risk in 2018, according to Eq. 13.

Hazard Assessment
The assessment 2018 of computed hazard along the littorals of the SoB domain is obtained according to Eq. 11 where the indices of impact density, vessel density and residual volume of oil are modulated by sea state conditions. Before displaying, the hazard values are normalized by the maximum values found in the area and its non-gaussian distribution, with positive skewness, gives a median value near 0.02. Thus, it is plausible to define a value of 0.03 as a threshold for hazard values indicating a high one.
In Figure 7 the hazard maps are shown for each vessel category that defines a different dangerous source for the SoB domain that is tanker (Panel A), cargo (Panel B), and passenger vessels (Panel C).
In the Panel A high values of hazard are found in correspondence of the industrial harbor of Porto Torres in Sardinia where noticeable traffic of tanker vessels is also found (Figure 4). Other littorals in the SoB domain display low hazard values in relation to tanker vessel traffic that, with some exception, mostly occurs offshore.
In the Panel B hazard values and their geographical distribution are comparable with those of Panel A, indeed high values are found in correspondence and near Porto Torres. In addition, medium-high values are found along the limits of the Gulf of Olbia, where the homonym touristic and relatively big city fosters cargo vessels traffic, and in correspondence of the SoB that is often crossed by cargo vessels, generating hazardous conditions also for the Lavezzi (Corsica) and La Maddalena (Sardinia) archipelagos.
In the Panel C, medium-high values of hazard are obtained for long traits of coast in correspondence of the Gulf of Olbia, the La Maddalena Archipelago and Bonifacio, where the traffic of passenger vessel is dense for roughly 6 months per year. The values are geographically scattered due to the complex coastal geometry that protect many littorals from prevalent meteorological and sea conditions.
Hazard values are used to compute the risk assessment 2018 for the SoB domain in a linear combination with ESI rank values, identifying those parts of the littoral that deserve a more attentive monitoring to reduce the effects of potential accidents at sea.

Risk Assessment
As done for hazard, the risk assessment in 2018 is analyzed on the basis of different sources of danger that, in this application, are represented by the passages of three distinct vessels categories.  In Figure 8, the risk assessment gives critical values along the littoral of the Gulf of Asinara with maxima in correspondence of Porto Torres and on the beaches in its surrounding. The littorals in the SoB and the Gulf of Olbia have low values of risk. Considering the threshold of 0.03, the percentage of littoral that could be impacted by oil stranding does not vary a lot during the year (Panel D), displaying two maxima near 15%, in March and July.
In Figure 9, critical values of risk are found along the littorals of the gulfs of Asinara and Olbia while they are scattered in correspondence of the SoB. Considering the threshold of 0.03, the percentage of potentially impacted littoral widely varies during the year (Panel D), displaying a maximum close to 30%, in February and a minimum close to 5%, in July.
In Figure 10, the risk assessment gives critical values along some littorals of the La Maddalena Archipelago and the northern part of Sardinia including some sites near Olbia. The percentage of the impacted littorals, for threshold of 0.03, increases during summer (from May to August), with a peak of 5% in August, and in autumn (from September to November), potentially related to changing sea conditions.

DISCUSSION
In this paper a novel operational numerical system for oil stranding risk assessment, with a dedicated web interface, is implemented. The aim of the work is to provide a user friendly and innovative tool to assess the risk at coast and support the preparedness to the management at-sea emergencies that could arise from oil spillage or floating pollutants in the Strait of Bonifacio. Its functioning core is based on a suite of coupled numerical models that can predict the sea state conditions and the potential risk of oil stranding along littorals.
Similar systems already exist and are widely used for comparable purposes (e.g., Azevedo et al., 2009;Wang and Shen, 2010;Olita et al., 2012;Canu et al., 2015). This demonstrates how the scientific community has a common view about the adopted methodology. However, local geophysical constraints, sources of hazard and vulnerable assets, must be considered to enrich the methodology and adequately implement the operational numerical system.
For these reasons, the presented operational numerical system includes a number of innovations for coastal applications. First, the adoption of a high-resolution grid to accurately reproduce the complex coastal geometry of the investigation domain and provide a detailed risk evaluation, at the same resolution of the geo-morphological information. Second, the adoption of timevarying maritime traffic density information as potential sources of hazard that, in the SoB domain, primarily originated offshore, where wide extents are yearly crossed by thousands of vessels, many of which with dangerous loads. Third, the web interface that was especially designed to define the offshore origin of the hazard, allowing the user to map the shape of the slick of pollutant, the time of the accident and, finally, displaying the risk assessment at coast. Fourth, an ad hoc formulation of hazard index and vulnerable assets.
The hazard index computation changes between authors and specific applications with the final aim of properly define the origin and the features of local hazardous sources. For example, three different approaches have been followed by Lan et al. (2015), Azevedo et al. (2017), and Al Shami et al. (2017). They introduced predisposing factors (e.g., favorable or not weather conditions) and control states (e.g., emergency input and emergency plan of sea area), historical records on past oil spill events or indicators of the time of exposure of the oil and susceptibility of shorelines to oiling.
In the present formulation, the hazard index accounts for: (i) the number of numerical particles that, representing floating oil droplets, hit the coast (i.e., impact density); (ii) the theoretical evaluation of the oil volume reduction at the sea surface that is due to weathering exposure and is computed as a function of particle age; (iii) the vessel density at the release point, that virtually assigns a high probability of accident where the maritime traffic is denser.
For what concerns the vulnerable assets, the adoption of social and economic sensitivities requires long data collection in the investigated area. Thus, in terms of a primary need for environmental protection, the adoption of geomorphological information gives, here, a reasonable measure of the littoral sensitivity to the oil exposure. Indeed, geomorphological features may promote or reduce the effect of oil stranding.
Littoral sensitivity is often assessed (Gugliermetti et al., 2007;Olita et al., 2012;Al Shami et al., 2017) with the adoption of the environmental sensitivity index (ESI) that is defined to classify shoreline sensitivity in relation to oil stranding (NOAA, 2002).
The combination of hazard index and ESI rank is a formulation for risk assessment (Eq. 13). The maps in Figures 8-10 uncover the geographic locations that are exposed to critical values of risk on the basis of the intrinsic sensitivity of the littoral and the potential threat derived by seasonal variability of the maritime traffic density.
The authors expected, for instance, to find critical values of risk in the surrounding of Porto Torres, if tanker routes are considered a source of hazard (Figure 8). However, it was not obvious to find critical values near Olbia and in correspondence of the La Maddalena Archipelago when cargo routes are considered a source of hazard (Figure 9). Also, passenger vessel traffic is very dense in correspondence of La Maddalena Archipelago with the result that long stretches of coast display critical values of risk (Figure 10). The south of Corsica displays non-critical values of risk along the littoral, with the exception of some of the numerous embayed sandy beaches that typically inhibit any natural recovery after unexpected accident of oil spillage at sea. Figure 11 synthesizes the results obtained by the risk assessment 2018 and highlights a seasonal variation of risk that should be considered for the implementation of any risk management plan in the area. By considering risk values above the critical threshold of 0.03, the Panel A displays the cumulated percentage of the impacted littoral due to each hazard source, tanker, cargo and passenger vessels routes. A seasonal modulation of such percentage is shown by each vessel category, presumably related to the variability of the vessel density values and the area crossed by vessels, but it is clear that navigation can really be a threat for an averaged impacted littoral of 24%, with a minimum of 16% on December and a maximum of 45% on February. By considering the average of those risk values above 0.03, the Panel B displays a risk that does not significantly vary during the year, except for tanker vessels, showing that the dangerous loads, moving offshore, give high risk values at coast but affecting a lower percentage of littoral if compared with cargo vessels (Panel A).

CONCLUSION
The oil spillage at sea remains a global threat affecting world coastal areas. The scientific community is engaged in the study of the whole set of chemical, physical, biological, social and economic phenomena that can result from marine pollution events and provides continual upgrades of the last innovations in terms of software, infrastructure and indicators. In fact, the appropriate integration of the current knowledge and scientific advances in relation to local environmental assets is the key of successful solutions for the management of marine pollution events. In this framework the authors proposed the inclusion of specific technical advances and research-based innovations into an operational numerical system for oil stranding risk assessment that is based on a core of hydrodynamic, wave and PTMs. Its functioning is demonstrated in the marine area of the SoB and its surroundings where the vessel traffic is a real threat for local assets and where the environmental agencies and policy makers ask for innovative responses from the scientific community. Technical advances include the adoption of a high-resolution unstructured computational grid for domain discretization and the implementation of a web interface to easily meet the needs of the final users of the system. Researchbased innovations includes the identification of the spatial and temporal variation of the vessel traffic and routes as source of hazard and the ad-hoc formulations for hazard indicators and vulnerability that for the first time in the SoB domain was included in the risk assessment computation. Vulnerability is computed on the basis of the main geomorphological features of the littoral that can favor or not the negative consequences of a potential oil stranding. This approach provided a clear identification of critical risk values also in correspondence of very low-anthropized littorals and high environmental and indirect economic value.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the online repositories. The names of the repository/repositories and accession number(s) can be found below: https://drive.google. com/drive/folders/1R6oE-IN4_6-XqHr496jlqrh2GOKdQ5M0? usp=sharing.

AUTHOR CONTRIBUTIONS
GQ and AC conceived the research approach and with SS provided pictures of all figures. GQ wrote the manuscript and processed the input and output data of the operational numerical system. AC applied the hydrodynamic, wave, and trajectory models in the SoB domain. SS also provided ESI rank information. AP and AC designed and implemented the Web GUI. RS applied the oil spill model MEDSLIK-II in a synthetic configuration, providing necessary indication to fit a theoretical curve of oil volume reduction due to weathering processes. GQ, AC, and AR designed the at-sea Lagrangian survey. AR provided the post-processing of the Lagrangian survey. All authors provided valuable contribution to the text contents.

FUNDING
This study was financed by the Projects SICOMARplus-SIstema transfrontaliero per la sicurezza in mare COntro i rischi della navigazione e per la salvaguardia dell'ambiente MARino and GEREMIA "GEstione dei REflui per il MIglioramento delle Acque portuali both funded by the EU program Interreg France-Italy Maritime 2014-2020."

ACKNOWLEDGMENTS
We thank the Centro Operativo per la Meteorologia of the Italian Aeronautiaca Militare to have furnished meteorological forcing, and the Office de l'Environnement de la Corse for friendly support.