A General Methodology for Beached Oil Spill Hazard Mapping

The current lack of a standardized approach to compute the coastal oil spill hazard due to maritime traffic accidental releases has hindered an accurate estimate of its global impact, which is paramount to manage and intercompare the associated risks. We propose here a hazard estimation approach that is based on ensemble simulations and the extraction of the relevant distributions. We demonstrate that both open ocean and beached oil concentration distributions fit a Weibull curve, a two-parameter fat-tail probability distribution function. The simulation experiments are carried out in three different areas of the northern Atlantic. An indicator that quantify the coastal oil spill hazard is proposed and applied to the study areas.


INTRODUCTION
Risk assessments have been widely employed by scientists and engineers as a tool to deal with harmful natural events (e.g., earthquakes) and anthropogenic events (e.g., oil spills) that we are unable to describe and/or predict with fully deterministic methods. A critical step within a risk assessment framework (e.g., ISO 31000), hazard quantification (or hazard mapping) is the probabilistic estimation on the frequency of occurrence and magnitude, based on historical observations of the events of interest. The seminal paper by Gutenberg and Richter (1944) clearly illustrates how the careful analysis of the probability distribution function (PDF) of seismic events, i.e., the probability of observing an earthquake of a certain magnitude in a limited area over time, led to an objective hazard quantification method and triggered comparisons within different tectonic areas around the globe (Shedlock et al., 2000). This well-established seismic hazard framework suggests that an appropriate statistical description of the variable of interest is of fundamental importance to classify hazard and then risk. As far as the authors are aware, there has been no study dedicated to the description of beached (or floating) oil concentration PDF, based on either simulated or observed data.
The oil spill hazard assessment framework lags behind these well-established hazard assessment fields. There is a severe lack of observational data, especially oil concentrations along coasts after spill events (in the paper also referred to as beached oil concentrations), eliminating the possibility of empirically describing the coastal oil hazard. Numerical ensemble oil spill modeling has been widely used to tackle the data gap, emulating reality, i.e., spill fate and generating information for hazard quantification (e.g., Price et al., 2003;Liubartseva et al., 2015;Sepp Neves et al., 2015;French-McCay et al., 2018). However, there is currently no standard way to describe the beached oil concentrations obtained by the ensemble experiments and, therefore, mapping the relative coastal hazard. Consequently, we are still unable to compare different hazard estimates, to evaluate objectively the completeness of the ensemble datasets and, most importantly, to be able to establish a global coastal oil spill hazard assessment framework.
The estimation of the probability distribution functions of atmospheric state variables and tracer/contaminant concentrations has attracted the attention of meteorologists for decades. It has long been recognized that wind speed PDFs fit a Weibull distribution (Justus et al., 1976;Pavia and O'Brien, 1986;Monahan, 2006), characterized by a large number of low magnitude events (e.g., low wind speed), large variance and, most interestingly, "anomalously" frequent high magnitude events, typical of fat-tailed distributions. Observational and modeling experiments (Tirabassi et al., 1991a,b) have shown that wind speed plays a major role in shaping tracer concentration distributions, which were also found to fit a Weibull PDF. In the ocean, recent papers have documented that surface coastal currents (Ashkenazy and Gildor, 2011), i.e., Gulf of Eilat/Aqaba, and global ocean currents (Chu, 2009) also fit a Weibull distribution, suggesting that ocean pollutants could also be similarly distributed as observed in the atmosphere.
The linkage between flow field characteristics and tracer concentration PDFs has been described by, among others, Hu and Pierrehumbert (2001) for stratospheric flows. They explain that every flow feature (e.g., eddies, filaments) has its own characteristics such as its shape and length, speed and direction, and how long it will last before disappearing/merging with other structures. Tracer concentration PDFs are primarily a result of different "parcels separation rates" taking place in the different parts of the concentration field, modulated by the flow field features they are exposed to. Areas with a relatively complex circulation are expected to present both features with lower (e.g., western boundary currents) and higher (e.g., filaments or eddies) "separation rates." In other words, as the flow field becomes complex and rich in features, the tracer concentration PDF increasingly depart from a Gaussian distribution. Therefore, we expect that if ocean feature-rich current patterns are Weibull distributed, the marine tracers concentrations, including beached oil, should be similarly distributed.
Recently Alves et al. (2016b) gave evidence that ocean current conditions determine the fate of the oil at the coasts. Currents are by themselves affected by bathymetry and they contain the information about the specific flow field that is consistent with the coastline constraint and the other forcings. In this paper we add to these findings by estimating the PDF of accidental oil releases, both at sea and along the coasts. Our definition is the one of "standard" accidental oil releases, corresponding to a volume of oil spilt equal to 10,000 tons, as pointed out but the analysis of Huijer (2005). Firstly, we will show that at the surface and in the open ocean, far from coastal boundaries, the oil concentration fits a Weibull PDF. Secondly we will simulate an ensemble of oil releases near the coasts and we will analyze the PDF of the beached oil. Three very different Atlantic coastline segments and ocean current regime areas will be considered for different current and wind conditions during the year 2013. A Weibull PDF function is estimated from the data and a coastal oil spill hazard index will be proposed in terms of the "Weibull tail distribution." After a description of the framework for oil spill hazard (section 2) and the ensemble experiments performed (section 3), the statistical distribution of oil concentrations in surface waters is analyzed (section 4). In section 5 we will then analyze the PDF for the beached oil on three coastlines around the Atlantic basin. Section 6 will offer the analysis of a new hazard indicator extracted from the statistical distributions of the beached oil in the three areas. Final overview of the results and future work is presented in section 7.

STANDARDIZED APPROACH TO EVALUATE OIL SPILL HAZARD
The framework for oil spill hazard and risk assessment has been defined by adapting the generic and recognized standard, i.e., ISO 31000, to oil spill model simulations so that the oil spill hazard assessment could be carried out in a scientific and standard way (Sepp Neves et al., 2015). The authors identify ensemble oil spill simulations as the method to produce the basic data to estimate the hazard.
Later, the technique of ensemble oil simulations was extensively analyzed for one coastal stretch, the Algarve coastal area in Portugal (Sepp Neves et al., 2016). The ensemble simulations were carried out using an extensive identification of modeling uncertainty sources: the meteo-oceanographic fields used for the oil spill simulations (different wind and ocean models were used), the inclusion of the Stokes drift, the oil spill release points, the oil spill model setup in terms of volume of oil spilled, time of spillage, and duration of the spill. Ten different model configurations were used to obtain the hazard of the beached oil at the Algarve coasts. The paper concludes that the beached oil distribution is non-Gaussian so that a hazard risk index which would use mean and standard deviations is not appropriate. In the same paper, Sepp Neves et al. (2016) show also a probability distribution of beached events that is analog of the pdf discussed and revealed in this paper (see Figure 10 of Sepp Neves et al., 2016). However, the specific PDF was not identified, and the hazard index not clarified.
In this paper we extend the simulations to 3 different coastal areas and to the open ocean, trying to show the generic character of the pdf resulting from the oil spill simulations as defined in Sepp Neves et al. (2015), Sepp Neves et al. (2016. In this paper we choose only to sample the uncertainties connected with the oil release points and the meteo-oceanographic conditions, together with the structure of the coastline. We consider operational volumes of oil spilt equal for all three areas. As shown in Sepp Neves et al. (2016) the inclusion of Stokes drift does not change the structure of the resulting beached oil PDF so we have not included the uncertainty connected to the Stokes drift.
The methodology of oil spill hazard mapping requires an extensive number of simulations due to the large number of uncertainties connected with oil spill conditions and model parameters. One principle of importance is the statistical independence of each release point simulation with respect to meteo oceanographic conditions that are responsible for the transformation and transport of the oil. In Sepp Neves et al. (2016) and in this paper we use the spatial and temporal decorrelation time of ocean mesoscale eddies as the threshold for successive simulations, but other criteria might be used if higher frequency ocean processes, such as tides, will be considered in the future. Furthermore, the length of the simulation should be a good compromise with ocean currents statistical independence and weathering processes time scale that has been estimated to be of the order of few days for similar amount of oil and duration of the spill (Coppini et al., 2011).
In this paper we would like to show that, independently from the coastal geometry, the beached oil PDF has the same structure, changing its parameters in dependence of the flow field characteristics and the wind induced transformations. The three areas considered in this study have very different mean current regimes: the Eastern Atlantic Archipelago is characterized by a relatively weak and broad southwestward flow, the Western Atlantic Island is characterized by the intense zonal Caribbean current flow and the Basilian coastal stretch is characterized by a western boundary current regime. Given these differences, the three areas differ greatly for the spatial scales and intermittency of the mesoscale eddy field, as well as the wind regimes.

ENSEMBLE EXPERIMENT SETUP
Ensemble oil spill simulations have been used in the past to assess hazard at the coasts and in the open sea (Price et al., 2003;Liubartseva et al., 2015;Sepp Neves et al., 2015;Olita et al., 2019) but surprisingly not much has been done with respect to the derived oil concentration frequency distributions. Here we use an ensemble simulation approach to define the PDFs for oil spill release events occurring near the coasts.
Our ensemble methodology consists of varying the oil release positions, changing the met-ocean conditions and collecting the concentration of oil either in the open ocean or along the coasts. This methodology allows to sample the uncertainty in hazard mapping due to the unknown location of the spill and different meteo-oceanographic conditions. Hazard mapping deriving from a sufficiently large number of met-ocean conditions and release points is then used to describe the PDF of oil pollution events at the coasts.
The oil spill transport and transformation model is MEDSLIK-II that solves an advection-diffusion equation for the oil concentration and its transformation processes (De Dominicis et al., 2013). The transformation of oil (or weathering) considers evaporation, emulsification, dispersion, and spreading processes connected to the slick. The advectiondiffusion part of the equation is solved by a series of lagrangian equations for the position increments, d x, of the parcels that compose the slick at the surface, e.g., where x = (x, y) is the horizontal position of the parcels, dt is the time increment, U( x, t) is the two dimensional ocean transport field, κ is a constant horizontal turbulent diffusivity value and z(t) are independent random numbers normally distributed, parameterizing lagrangian diffusion. The model considers only surface releases, the slick is evolved considering four type of parcels: surface, dispersed, sedimented and beached parcels. For the surface parcels, Equation (1) is used with a prescribed κ (taken to be equal to 2 m 2 s −1 and current velocity analysis fields, as described in Table 1. Earlier studies (Dominicis et al., 2014) of lagrangian diffusivities found κ ∼ 1 − 100 m 2 s −1 and here we have chosen the lowest value possible in order to diminish the importance of diffusion with respect to advection by the currents. The dispersed particles are generated by the oil transformation processes (De Dominicis et al., 2013), they evolve separately from the surface oil parcels and they compose a volume of submerged oil that we do not consider in this paper except if the dispersed oil becomes beached oil. In MEDSLIK-II oil transformation processes are associated only to density, expressed as American Petroleum Institute (API) gravity.
Oil parcels arriving at the coasts become beached oil parcels when their trajectories cross pre-defined coastal segments. The permanent beaching and the probability of wash-off will depend on the coastal type, as explained in De Dominicis et al. (2013), however the coastal oil holding capacity is the same for all types of coasts and it is currently set to 5, 000 bbl/km. In our analysis we use the total beached oil, the sum of the permanent and re-detachable oil. The simplistic beaching model is in accordance with the limited temporal and spatial resolution of the employed meteo-oceanographic models and in the past it has generated satisfactory results in real spill cases at similar scales (Coppini et al., 2011).
All experiments are carried out supposing several Release Points (RP) off the coastlines (Figure 1). The RP grid encompasses the main maritime routes simulating potential oil releases from the maritime corridors. Given the amount of oil released, 10 4 tons ( Table 1), our simulations encompass the release from potential accidents (Huijer, 2005). The current field U( x, t) used in (1) is derived from the surface analysis fields of the Copernicus Marine Environment Service (CMEMS) at 1/12 • spatial resolution (Lellouche et al., 2013) 1 . The atmospheric winds and air temperature used in the weathering processes are derived from the European Center for Medium range Weather Forecast (ECMWF) ERA-Interim at 1/8 • (Dee et al., 2011). The input winds are given every 6 h and linearly interpolated to the oil spill model time step which is 1 h.
For each release point, 240 h-long oil spill simulations were launched every 10 days throughout year 2013, using the experiment configuration described in Table 1 for three different areas in the Atlantic ocean. All different areas used the same wind and ocean current data set, the same spilled volume, type and duration of the spill (see Table 1). The coastline segments were derived from the global Wessel and Smith (1996) data set which represent the coastlines with polygons of a few hundred meters resolution for the world ocean. Coastal segments were considered at 1 km resolution for the three study areas because this is recommended if details of the coastline are not known from complementary in situ data.
The statistical independence of each release point simulation, i.e., the occurrence of one beached oil event not affecting the probability of occurrence of the other due to the same release point, was preserved in two ways. Firstly, the time interval between consecutive simulations from a single release point was longer than the typical Lagrangian decorrelation time in the Atlantic ocean (∼5 days, see Maximenko et al., 2012). Secondly, the distance between neighboring release points (i.e., 25 km) was chosen based on estimation of the mesoscale spatial decorrelation scale from satellite altimetry that is at its minimum 50 km (Le Traon et al., 1990). The 10 days long simulations were chosen as a good compromise between this statistical independence and the transformation processes occurring at the surface that will evaporate and disperse the oil in the first few days for the amount and modalities of oil released in the simulations (2 days release time of 10,000 tons), as shown for example in Coppini et al. (2011).

THE OPEN OCEAN SURFACE OIL CONCENTRATIONS PDF
The first aim of the paper is to show the PDF of surface oil in the open ocean, due to a single oil spill event repeated every 10 days for the entire year 2013. A single release point is chosen off the Portuguese coast, located in the open ocean flow field (Figure 2A) using the model configuration described in Table 1. The experiment considers a "virtual" coastline, sampling line, defined 50 km southward of this release point, where the oil concentration is monitored to study its statistical distribution due to the oil advective-diffusive dynamics with and without the oil transformation processes (weathering and no-weathering experiments). All non-zero surface oil concentration values observed along the sampling line grid cells are stored every 6 h throughout the duration of each simulation.
In total, weathering and no-weathering experiments generated 3,779 non-zero oil concentration values each along the sampling line grid boxes and their concentration histograms are presented in Figure 2. Both histograms show that most of the events accumulated in low concentration bins but that in general, the range of concentrations observed is large, spanning two orders of magnitude for the oil concentration. Very high concentration events were found to be present, depicting a fat-tailed histogram.
Based on the histogram characteristics we choose a Weibull distribution function, W(x; β, η) for a maximum likelihood estimation of the PDF: where x is the oil concentration, β and η are the shape and scale parameters respectively. In the Weibull function (2), the shape parameter is the one that controls how different from a Gaussian distribution is the PDF. For β less than 1 the distribution has its maximum at values close to zero. For β = 1 it is an exponential and for β greater than 1 the distribution has a single maximum for oil concentrations for relative large positive values and exhibits a marked asymmetry. The Weibull distribution function has a wide range of applicability and it was in fact invented by Weibull and Stockholm (1951) to describe the failure of mechanical systems under various types of stresses or the size distribution of fly ash.
Here we show for the first time that the surface oil in the ocean fits the same distribution as the currents (Chu, 2009;Ashkenazy and Gildor, 2011). It is interesting also to note that the Langevin-type equations, like the ones used in our oil spill modeling equations (1), produce probability distribution functions that are skewed and heavy-tailed (Sardeshmukh and Penland, 2015). Thus it    The error in the scale and shape parameter has been computed by a bootstrap algorithm subtracting ∼ 10 − 20% of the input data. In addition to the scale and shape parameters, we calculate also the mean and standard deviation of the estimated distribution, defined seems likely that an active ocean tracer, such as oil, would fit a statistical distribution with a heavy tail, i.e., have a distinct non-Gaussian behavior. The estimate of (2) to the sampling line events is shown in Figure 2 demonstrating that Weibull is a good distribution function for the oil arriving at the sampling line. The estimated shape parameter value is β = 0.399 ± 0.005, indicating a large peak at the low concentrations, for the weathering and noweathering simulations. The scale parameter is instead different between weathering and no-weathering simulations because it is connected to the distribution mean value and it is η = 0.50 ± 0.03 tons/km 2 and η = 0.93 ± 0.06 tons/km 2 for the weathering and no-weathering experiment respectively.
We conclude that the oil concentration distribution from a single point source release and along an offshore sampling line has the form of a Weibull distribution that is peaked around the lowest concentrations.

THE BEACHED OIL CONCENTRATION PDF
In this section we will show that, notwithstanding the complexity of beached oil processes, the beached oil distribution along the coastlines of three different areas in the Atlantic has a universal structure similar to the surface oil PDF found along the sampling line in the open ocean (section 4).
Following the experimental set up described in section 3, we collect the beached oil parcels on each segment of the coastline for each run. The number of RP points in the three areas are reported in Table 2 together with the number of simulations and beached oil events in each area. In Figure 3, we present the histograms of the beached oil concentration distributions for the Eastern Atlantic Archipelago, Western Atlantic Island, and state of Bahia (Brazil) coastlines.
The Weibull distribution function seems to estimate well the beached oil data, even if some discrepancies appear at concentration values between 0 and 50 tons/km. The scale, η, and shape, β, parameters are reported in Table 3. It is clear again that the shape parameter is less than 0.5 for all the three areas, showing that the maximum number of events occur at the lowest values of the distribution. It is puzzling that the shape value is the same for all the areas, within the calculated errors, and it is also similar to the one found for the open sea conditions in section 4. We argue that our results show that the beached oil distribution does not depend on the specific coastal geometry but only on the type of flow field conditions that impinge on the coasts. The coastal currents are themselves affected by the type of bathymetry and coastline, so the dependence is indirect. Furthermore the beached oil PDF depends also on the winds, both directly and indirectly: the former through the transformation processes and the latter because currents are themselves driven by winds.
The histograms shown in Figure 3 could be also displayed on a logarithmic scale that would enhance the fat tail structure of the three distributions (Figures 4A,C,E). To check visually if our data set comes from a population that would logically be fit by a 2-parameter Weibull distribution we use the Weibull plot (Nelson, 1982). It is clear from Figures 4B,D,F that the assumption of a Weibull distribution is reasonable and that there are no outliers.
In Table 3, we list also the mean and standard deviation of the distributions from the estimated Weibull function. The Eastern Atlantic Archipelago and Western Atlantic island have similar "mean" value of beached oil concentration, i.e., 23 tons/km, while Bahia ∼ 30% less, i.e., 17 tons/km. We conclude that the mean value of beached oil, due to surface current transport and wind induced transformations, is quite different between the two northern Atlantic areas and Bahia.
The mean of Weibull distributions with a β parameter smaller than 1 is shifted toward low concentrations. Thus, the distribution mean is not really sufficient to assess the hazard which should consider also the large events in the PDF tail. In order to do so we need to consider the cumulative distribution and define an appropriate indicator.

A GENERAL INDICATOR OF COASTAL OIL SPILL HAZARD
We will now show how beached oil concentration PDFs can be used to describe the hazard in a general way, paving the way for a quantitative method to answer a key question related to the management of the oil spill risk: what is the beached oil spill hazard for different areas and which area has the highest hazard?
Having now the PDF which describes the beached oil events from a set of RPs, simulating potential releases from accidents, we can define an index that contains the information about the large pollution events, concentrating on the tail of the Weibull distribution. It is indeed reasonable to assume that the hazard of large oil releases near the coasts should be evaluated by estimating the extreme events of the distribution.
To describe the extreme events it is then necessary to find how many events are contained in the tail with respect to the overall distribution events. To this end, we can use the Weibull tail distribution H, defined as: where F(x cut ) is the Weibull cumulative distribution function for the beached oil concentration x for which x x cut . When H is large, it means the distribution tails contain a relatively large number of high concentration events, thus the hazard is large while the contrary is true for small values of H.
The key parameter here is the threshold value x cut chosen for the estimation of extreme events. For illustrative purposes, a concentration of 25 tons/km was applied following the International Tanker Owners Pollution Federation (ITOPF) Technical Information Paper on "Recognition of oil on shorelines" (ITOPF, 2011). The threshold is proposed for classifying impacted coastlines as heavily oiled.
The tail distribution H was computed for each area (Table 4) for the chosen value of the threshold. The Western Atlantic island emerges as an area of relatively high hazard with respect to the Eastern Atlantic Archipelago and the Bahia areas. This means that large beached oil events are more likely to occur in the Western Atlantic island, about ∼ 12 − 14% more frequently than in the Eastern Atlantic Archipelago and the Bahia coasts,

CONCLUSIONS
This paper develops a straightforward and objective method to quantify the coastal oil spill hazard based on ensemble oil spill experiments which sample the uncertainties associated with oil spill accidental releases in the ocean areas next to the coasts. Ocean current variability at the mesoscales affects the flow field and in order to have a realistic sampling of the oil transport variability, ensemble simulations have been extended for an entire year, sampling the seasonal variability of the flow field at mesoscale permitting resolution. It has been found that both oil in the open ocean and beached oil concentrations are successfully described by the Weibull distribution. Such distribution is characterized for these cases by a shape parameter which is less than 1, thus having the maximum beached oil events at small concentrations. The large beach oil concentrations are contained in a "fat tail" which characterizes this distribution. Concentrations at the coasts span two orders of magnitude values, from 1 to 100 tons/km. Previous studies have indicated that also the winds, currents, and tracers in the atmosphere and oceans distribute as Weibull but this has never been shown for active ocean pollutants, such as oil.
Furthermore, the paper proposes a new hazard index for beached oil which allows to intercompare different world ocean areas and estimate the different hazards. The index is based upon the fit of a Weibull PDF to the simulation data and extracting the tail distribution for large events.
We would like to point out that the present study has potential limitations due to the relatively low number of uncertainties considered in our simulations. In the future, experiments should probably vary the simulation length, should check different oil types and beaching algorithms and use different meteooceanographic conditions including higher frequency winds and tides. This paper does not discuss how exactly the flow field modulates the oil spill hazard through advection-diffusion processes in the ocean. Nor our oil spill hazard estimates for the three regions examined are conclusive for future management decisions, especially due to the limited number of oceanographicweather conditions considered in the study. Such estimates can be obtained with the methodology described here, but will require more extensive investigation and analyses. It is important here that a common quantitative framework for intercomparison of different world ocean areas has been found. Furthermore, to evaluate risk, vulnerability data need to be acquired and used to modulate the hazard that we have computed here (Alves et al., 2016a).
Future work will consider an in-depth study of the ocean flow field parameters and how they modulate the coastal oil spill hazard. Additionally, the spatial resolution of the oceanographicweather fields used in our experiment offers reliable answers only for mesoscale dominated flows and coarsely resolved coastal morphology and bathymetry. The increase in resolution of the ocean and atmospheric models, and therefore their ability to reproduce smaller scale features and higher frequency variability, should impact the oil concentration PDF and, consequently, the oil spill hazard estimates. Furthermore, uncertainties connected to the amount of oil released should be also considered in the future, i.e., cases of larger oil spill releases.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
All authors contributed to conception and design of the study. AS generated the database and wrote the first draft of the manuscript. NP and AN wrote sections of the manuscript and fully revised the manuscript. FT made the additional figure on the revised manuscript and discussed the meaning. All authors contributed to manuscript revision, read and approved the submitted version.

FUNDING
This work has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement No. 633211, project AtlantOS.