SHORT-TERM AND MEDIUM-TERM FORECAST OF OIL SPILL TRAJECTORIES: APPLICATION TO LOCAL AND REGIONAL SCALES

The number of accidental oil spills affecting the coastal environment in recent decades has led to a growing concern regarding oil spill preparedness and response, and has motivated the development and implementation of different tools to be used in these emergency situations. Oil spills can arise from a number of different sources ranging from offshore spills (e.g. drilling, collision or accidents of vessels carrying crude oil) to coastal events (e.g. oil-loading and –unloading operations in local ports or even deliberate and illegal spills to clean the ship’s tanks).

Operational tools based on numerical models that provide real-time forecasts of oceanographic variables and oil spill evolution have demonstrated their usefulness in recent oil spills (Castanedo et al., 2006, Liu et al., 2011). In order to improve emergency oil spill response, many national operational oceanographic systems have been setup around the world focusing on short-term forecasts (up to 5 days). It is worth mentioning that most of these systems have been developed at regional scales (O(km)), and there is a lack of systems suitable for local scales such as estuarine environments or ports (O(10m)). Offshore oil spills (as in the Prestige in Spain or the Deepwater Horizon in the Gulf of Mexico) have demonstrated the importance of having larger prediction horizons to model oil spills that remain longer in the marine environment. In such cases, response planning would depend largely on understanding its likely evolution over a longer timespan (15 days to 2 months). Therefore, tools for oil spill preparedness and response must take into account the management of the immediate response (short-term) and the planning for oil spills that may be more persistent (medium-term).

To develop this area, this paper presents two methodologies providing short and medium-term forecasts of oil spill trajectories at local and regional scales. For short-term predictions (within 48 hours), a high-resolution operational oil spill forecast system based on dynamic downscaling in Belfast Lough (Northern Ireland) is presented. For mid-term predictions (15–60 days), a stochastic methodology based on the combination of met-ocean databases, clustering methods, autoregressive logistic regression models and oil spill numerical modelling is proposed. The mid-term methodology is applied to the Bay of Biscay (Northern Spain) to simulate oil spills observed during the Prestige accident.

Belfast Lough is a shallow semi-enclosed marine bay situated on the east coast of Northern Ireland at the mouth of the river Lagan, with the city of Belfast at its head (Figure 1). The Port of Belfast is a busy shipping port, which manages over 80% of Northern Ireland’s petroleum and oil imports and handles over 5,000 vessels each year. In terms of ecological sensitivity, Belfast Lough has many conservation designations. There are currently two Natura 2000 (EC Habitats Directive) designated sites within Belfast Lough: the Belfast Lough Special Protection Area (SPA) and the Belfast Lough open water SPA, the latter of which covers most of the inner lough adjacent to the Port. In terms of socio-economics, subtidal aquaculture for bottom-grown blue mussel (Mytilus edulis) is a notable activity in the inner lough, with licensed beds occupying over 50% of the seabed in this area. In addition to leisure and recreational sailing, pot fishing, scallop dredging and bottom-trawling for Nephrops norvegicus in the outer lough are activities sensitive to oil pollution.

The purpose of the system is to provide two-day forecasts of oil spill fate and transport in Belfast Lough. To achieve this objective, the system comprises an operational oceanographic module coupled to an oil spill forecast module, integrated in a web application that allows the user to run simulations of hypothetical spills as well as model oil spill evolution in real-time during an emergency situation (Figure 2).

The oceanographic module is based on Delft3D model (Roelvink and van Banning, 1994), which is a numerical hydrodynamic model that solves the Navier-Stokes equations for shallow water with hydrostatic pressure hypothesis and Boussinesq approximation. In order to obtain the high resolution required inside the lough, the model domain is discretized using a curvilinear grid composed by 222 × 61 elements from Belfast Harbor to the coast of the Galloway peninsula (south west Scotland). The grid presents a variable resolution increasing from 2 km in the North Channel, to 500 m in the outer lough, and reaching 20 m in the harbor area. The model domain is vertically discretized in 5 vertical σ-coordinate levels, with layer depths ranging between 5% (surface) and 35% (bottom) of water column. The hydrodynamic model is nested to the regional MyOcean-IBI (Iberic-Biscay-Irish) forecasting system (http://marine.copernicus.eu/) belonging to the European network COPERNICUS. The system interpolates to the open boundary grid nodes sea surface height data hourly, and 3D sea temperature and salinity data provided by the MyOcean-IBI system on a daily basis. Atmospheric forcing parameters are provided by the UK Meteorological Office (http://www.metoffice.gov.uk/). The forecasting system runs at 8 km spatial resolution in the study area and provides 48-hour atmospheric forecast with hourly resolution.

The operational oceanographic module provides a two-day forecast of currents, sea level, temperature and salinity with 15 minutes temporal resolution. These data and the meteorological forecast are used to force the oil spill forecast module.

The oil spill forecast module uses the TESEO oil spill model (Abascal et al., 2007), which consists of a transport and a weathering module to represent the evolution of oil spilled in the marine environment. The transport model derives from the two-dimensional Lagrangian transport model developed by the University of Cantabria as part of the operational forecasting system created to respond to the Prestige oil spill (Castanedo et al., 2006). In this model, the drift process of the spilled oil is described by tracking numerical particles equivalent to the oil slicks. The evolution of the horizontal position of each particle is defined by the combined effect of advective and diffusive velocities. In this system, the advective velocity is calculated as the linear combination of current and wind velocity. Given the relevance of tidal currents and winds to the hydrodynamics of the Lough, the Stoke’s drift is not integrated on the oil spill transport. The turbulent diffusive velocity is obtained using a Monte Carlo sampling in the range of velocities that are assumed proportional to the diffusion coefficients. A total of 1000 independent numerical particles were used in the simulations and a 60 s time step was used to calculate the time evolution of the particle positions. Based on the calibration of the model with drifting buoys, the wind drag coefficient was set to 2.5 % of the wind speed and the diffusive coefficient was set to 33 m²s⁻¹ (Abascal et al., 2016).

The weathering module includes the process of slick spreading under gravitational and viscous forces, evaporation of the lighter components of the spilled oil, the entrainment of water forming emulsions, and changes in physicochemical properties. The oil stranding process is also taken into account in the model; when the oil reaches the shoreline, it is considered to be fully entrapped and it remains on the coast. However, the oil can be re-incorporated into the sea to properly model the wet and dry areas are a feature of estuarine and coastal environments. Both the transport and the weathering components of the model have been validated in previous studies (Abascal et al., 2007, 2009; Castanedo et al., 2014).

The oil spill module provides two-day forecast of the transport and dispersion of the oil spill, the temporal evolution of the oil properties (density, viscosity and water content) and the mass balance of the spill (amount of evaporated oil, beached oil, oil remaining on water, emulsified product that reaches the coast and emulsified product that remains at sea).

The system was validated using surface drifter buoys released in Belfast Lough during three experiments performed between September and November 2014. The drifters (model MD03, http://www.albatrosmt.com/), consist of a of 25 cm by 10 cm diameter cylinder with foam protection (see Figure 3) that transmit their GPS position via the GSM network to a modem connected to a PC. As an example this paper presents the validation of the system with a drifting buoy deployed on November 4^th, 2014 between 10:30 and 16:40.

Surface currents and wind data provided by the operational oceanographic system were used to force the numerical simulations. The comparison between the actual and simulated trajectory is presented in Figure 3. In this figure, the observed buoy positions are shown as green circles, the simulated particles are shown as white points and as orange points for the last time step, and the center of mass of the simulation is shown as red points. As can be seen in the figure, the drifter trajectory lies inside the area simulated by the model. The distance between the final drifter position and the center of mass of the particle cloud at the last simulation step is about 540 m, which represents the 16 % of the distance covered by the drifter (3500 m). This normalization is convenient since a pure deviation measurement can indicate a better or worse performance of the transport model depending on how far from the initial position the distance is obtained. Thus it is shown the capability and the range of accuracy of the proposed model as an assessing tool during an oil spill accident in the Belfast Lough.

Further research is required including actual trajectory data coming from different types of drifters. This data would better emulate a real oil spill drift, which is a function of the weathering state and the physico-chemical characteristics of the spilled oil.

The Bay of Biscay is located in the northeast Atlantic Ocean (see Figure 4). It is bounded to the east by the west coast of France and to the south by the north coast of Spain. The general circulation in this area is connected with the general circulation of the North Atlantic (wind and sea level pressure) at the surface and especially on the platform (Ferrer et al., 2009). This area has been affected by several important oil spill accidents as the Erika (1999) or the Prestige accident (2002) which illustrates the nature of the risk. The Prestige disaster (Castanedo et al., 2006) demonstrated that it is possible for an oil spill occurring in the Atlantic Ocean to enter and impact the Bay of Biscay. During this event, more than 1000 km of coastline and a huge variety of habitats were affected, from the Galician coast (northwestern Spain) to the French coast.

The proposed methodology is broken down as follows (see Figure 5):

• Reanalysis databases: use of high resolution long-term historical hourly databases of: 1) sea level atmospheric pressure for the period 1957–2001 from global reanalysis-I dataset NCEP/NCAR (https://www.esrl.noaa.gov), 2) wind from the SeaWind-EraInterim reanalysis (Menéndez et al., 2014) available in Europe from 1989 up to 2009 and 3) currents from the GOS2.1 reanalysis (generated by wind and pressure gradients) (Cid et al., 2014) available in Southern Europe from 1989 up to 2009.

• Classification and selection of atmospheric patterns: The atmospheric-pressure and wind reanalysis databases are classified by application of principal component analysis (PCA) and classification techniques (K-means clustering) (Camus et al., 2014). The procedure established is as follows: (1) A PCA is applied to the monthly sea level pressure anomalies (MSLPA) (predictor); and (2) classification of the daily mean wind states (DMWS) (predictand), this gives K wind states representing the complete database.

• Autoregressive logistic regression (ALR) probabilistic forecasting of met-ocean conditions: in order to obtain the probabilistic prediction of the evolution of wind conditions and surface currents, statistical techniques based on autoregressive logistic regression processes are applied. This procedure is structured in the following steps: i) Model Fit. To adjust the logistic model, we take into account the variables of influence such as seasonality, long-term trends, covariates (MSLPA fields, climatic indices…) and autoregressive processes; ii) Logistic simulation. Once the fitting process is completed, the model allows the simulation of categorical synthetic time series of wind states. These have been generated using the Monte Carlo method; iii) Probabilistic current forecasting. Obtaining the surface currents associated with the wind fields calculated in the previous step. From the logistic model, N simulations of the evolution of the wind and currents conditions in the established prediction period (S days) are made, obtaining a set of N x S wind and currents simulations. We are currently working to develop an analogous methodology to predict ocean currents based on ALR models.

• Probabilistic forecasting of the oil slick trajectory: working from the stage 3 results, N simulations are made of oil-slick trajectories during a period of S days. The simulations are carried out using the TESEO oil-spill transport model (Abascal et al., 2009). This gives us S x N simulations of equiprobable trajectories, which are then used to calculate the probability of spill pollution.

The methodology has been applied to simulate the evolution of the trajectory followed by oil slicks observed during the Prestige accident, and the results are presented below.

• i) Predictor characterization: The first step of the methodology is to apply a PCA to the MSLPA fields to obtain the dominant spatial variability patterns (EOFs) and their corresponding temporal coefficients (PCs), to reduce the dimensionality of the temporal MSLPA fields whilst preserving the maximum variance of the data. The first 10 modes (see Figure 6) from the PCA analysis (explaining the 97% of the variance) are selected as the predictor of the wind data in the North Atlantic. In the step 2 of the methodology, the 10 PCs are the covariates of the ALR model.

• ii) Predictand characterization: Before the data clustering, a PCA is applied to the DMWS in order to reduce the temporal dimensionality of the data and to improve clustering. From this analysis, 23 PCs were obtained explaining the 99% of the variance. The K-means algorithm is then used to classify wind states. In this case, 20 years of daily mean wind data (n= 7638 days) are classified into K= 49 clusters (see Figure 7). Once the classification is obtained, each daily record of wind belongs to one of the 49 clusters named “Wind States” (WSs). A time series of WSs is generated which becomes an input to the ALR model. The event sequence time series represents the evolution in time of the WSs (Antolinez et al. 2016).

This step of the methodology consists of modelling the WS sequence using ALR that allows the simulation of synthetic sequences of WSs while taking into account different covariables such as: seasonality, PCs of MSLPA, autoregressive terms of WSs and long-term trends. These kinds of models are useful when working with categorical variables, such as the results of clustering data (Guanche et al., 2013).

Let Y_t;,t=1,…,N be the observations of WS at time t, with the following possible outcomes Y_t ∈{1,…,n_s} related to each WS. Considering X_t, t=1,…,N to be a time-dependent vector of covariates with dimensions (n_c x 1) (i.e., seasonal cycle, principal components PC_k(t) of synoptic circulation, long-term trend, etc.), the ALR model is stated as follows:

where t is given in years, corresponds to annual mean values, and are the harmonic amplitudes. X_i are the vectors of the n_c covariates, and is the corresponding parameter vector. represents the long-term trends effects. Y_j−d is the WS of the previous d-states, γ_j−d is the parameter associated with the previous d-states and D corresponds to the number of previous states that influence the actual WS.

• i) Model Fitting: To perform the estimation of the parameters for the model, a maximum likelihood estimator is used which requires the definition of the likelihood function (Hosmer and Lemeshow, 2000). In this case, the criteria to choose the final model are based on the statistical significance, and specifically, the likelihood ratio (LR) statistic (see details in Guanche et al., 2013). The goodness-of-fit is obtained calculating the deviance ratio (ΔD), which measures the change in the quality of fit between two different parameterizations, and the chi-square distribution (χ²) with Δdf = Δn_p × (K − 1) degrees of freedom, where Δn_p is the difference in terms of numbers of parameters for both parameterizations, and K is the number of clusters of WSs selected. To identify the best model, we have evaluated the goodness-of-fit of different nested ARL models related with different predictor variables. The results (not shown here) determine that the final model takes into account the following predictor variables: seasonality, covariables of MSLPA and two autoregressive terms. The model is formed by a total of n_p = 109 parameters.

• ii) Model simulation and validation: Once the model has been adjusted and the coefficient matrix obtained, the probabilistic forecast of the evolution of the wind states can be obtained. To do this, the synthetic time series of the WSs are generated using the Monte Carlo method. Due to the stochastic nature of the process, a total of 1000 simulations of the evolution of the wind conditions for each of the simulation scenarios were run. To validate the model, we compare the probability transition matrix between WSs for historical and simulated data. The validation is carried out over 20 years of daily wind data (n = 7638). The probability transition matrix represents the probability of changing from group i to group j between consecutive days. Thus, in the case of having 49 WSs, the transition matrix (TP) has dimensions 49×49 (N= 2401 elements), and each cell TP_ij is the probability of changing from wind state i to wind state j. Figure 8 shows the scatter plot related to the TP matrix. The diagonal matrix represents the probability of staying in the same group. As can be observed, the model is able to reproduce correctly the transitions between wind states within the sequence. The correlation coefficient (CORR) is greater than 0.99, with a confidence interval CI < ±0.05 and RMS error lower than 0.3%.

  Once the proposed model is validated as suitable for obtaining the synthetic time series of wind types, the next step develops the simulated wind series, which will provide the probabilistic predictions of the wind fields for each of the N simulations of the periods selected for the study.

• iii)Probabilistic current forecasting: Once the evolution of the wind is known, the synthetic series of daily mean surface currents are obtained. To do this, we associate the daily mean currents obtained from GOS2.1 database corresponding to the same date as the WSs.

This methodology has been previously applied to simulate the main black tide (MBT) of the Prestige accident observed on 19 November 2002 when the ship split in two and sank at 130 nautical miles from the Galician coast (see Figure 9). Taking into account the previous results, 1000 simulations of the oil slick trajectory are carried out using the TESEO oil-spill transport model (Abascal et al., 2009) for the period comprised from 19/11/2002 to 01/17/2003 (60 days of simulation). This gives us 60 × 1000 simulations of equiprobable trajectories, which are then used to calculate the probability maps of spill pollution. Figure 9 shows the cumulative probability maps of pollution compared to observations of spill stains (MBT) during the incident. The cumulative probability maps have been obtained taking into account a grid with 0.5° of cell size, and in addition in each simulation we have run the model with 1000 particles. This shows the evolution of the probability maps at three (figure 9 (a)), five (figure 9 (b)), eight (figure 9 (c)) and twelve (figure (9 (d)) days after of sinking. Each map show the cumulative probability of the predicted oil versus the observed oil for the prediction date and the observations of the previous days to better see the evolution of the time progression of the prediction. It is observed as the highest probability of contamination coincides with the path followed by observations of the spots, whose origin is the area where the ship sank, and their evolution towards to northeast impacting on the southern coast of Galicia, on 30 November 2002 (see figure 9 (d)). In general, the probability of contamination agrees with the areas where the spots were observed, being important the pollution in the Galician coast (Montero et al., 2003). The results obtained show both the capacity of the developed methodology as its usefulness to provide necessary information in the planning and response before an oil spill in the marine environment.

This work has been partially funded by the European Transnational Programme (Atlantic Area) under the SPRES project and by the Ministry of Economy and Competitiveness (MINECO) under the research project TRA2014-59570-R (PLVMA3D). B.P. and M.C. would like to thank MINECO its support within the FPI Program. J.A.A.A. would like to thank MECD (Ministerio de Educación, Cultura y Deporte, Spain) its support the FPU Program (BOE-A-2013-12235). The authors would like to thank to Ángel David Gutiérrez (Sands Solutions and Services, SANDS CORP. S.L.) for his contribution in the implementation of the operational system in Belfast Lough.