A HIGH RESOLUTION OPERATIONAL OIL SPILL MODEL AT SANTANDER BAY (SPAIN): IMPLEMENTATION AND VALIDATION

Up to now, operational oceanography systems (OOS) which provide decision makers with oil spill trajectory forecasting, have demonstrated their usefulness in dealing with recent crisis involving dramatic environmental and socio economic impacts (Castanedo et al., 2006). Accordingly, numerous national and regional OOS have been set up all over the world in order to facilitate and support an adequate response and preparedness (e.g. Olabarrieta et al., 2007). However, occasionally, the available spatial data resolution of these systems is insufficient to be used for coastal applications, and therefore, this information must be transferred to shallow water, in order to increase the spatial resolution (namely downscaling).

With the aim of increasing our capacity to respond to oil spill pollution at the coast, a set of high resolution operational oceanography systems (HR-OOS) are being currently developed in the framework of the SPRES (OIL SPILL PREVENTION AND RESPONSE AT LOCAL SCALES) project, co-funded by the European Transnational Programme (Atlantic Area). This pilot study is being undertaken in four European sites: Santander's Bay, Aveiro's Lagoon, Belfast's Lough and the Port of Falmouth.

Certainly, one of the main challenges when approaching the HR-OOS, is to correctly model local scale physical processes whilst providing a response fast enough to be used in an emergency situation.

Bearing this in mind, the present paper describes the OOS that is being built in the Bay of Santander. Special attention has been paid to show the numerical models setup and their calibration.

Santander's Bay is one of the largest inlets (2270 ha) on the northern coast of Spain (Gulf of Biscay) (see Figure 1). The bay provides a natural shelter from the waves of the Gulf of Biscay which arrive from NNW and have an annual average significant wave height of 1 m with winter storms waves of 4 m. Tide is semidiurnal with a mean and spring tidal range of 3 and 5 meters respectively. Four rivers discharge into the Bay. Their average discharge does not reach 15 m³/s, although under flood conditions the discharge may be as large as 400 m³/s.

Seventy five per cent of the total extension of the estuary is intertidal area. Subtidal zones are dominated by shallow waters, with maximum depths of 10–12 m occurring along the navigation channel.

Regarding oil spill risk, the Prestige disaster (Castanedo et al., 2006; 2009) proved that it is possible for an oil spill occurring thousands of kilometres away to enter and impact the Bay. Moreover, the harbour of Santander and other companies located within the Bay manage oil products and are a real threat for this valuable environment which supports many socio economical activities such as shellfish extraction, fishing and recreational uses.

The main goal of the system presented here is to provide short-term (within 48 hours) oil spill trajectory forecasting at local scale (Santander's Bay, Spain).

The system is based on nesting high resolution models and is made up by the following components: (1) Daily boundary conditions (sea level, ocean currents, salinity and temperature) and meteorological forcing are obtained from the European network MYOCEAN (http://www.myocean.eu/) and from the Spanish met office, AEMET, respectively; (2) COAWST modelling system which allows to exchange data fields between the ocean model ROMS, the atmosphere model WRF and the wave model SWAN, is the engine of the OOS (at this stage of the project only ROMS is implemented); (3) an oil spill transport and fate model, TESEO (Abascal et al., 2007); (4) a web service that manages the operational system and allows the user to run hypothetical as well as real oil spill trajectories using the daily forecast of wind and high resolution (O(20 m)) ocean variables carried out by COAWST. Figure 2 shows a scheme of the structure of the system.

As can be observed, the HR-OOS is composed of two main modules: the HR operational oceanographic module (HR-OOM) and the HR oil spill forecast module (HR-OSM).

The HR-OOM is based on COAWST modelling system (A Coupled-Ocean-Atmosphere-Wave-Sediment Transport Modelling System) (Warner et al., 2010). COAWST is comprised of the Model Coupling Toolkit to exchange data fields between the ocean model ROMS, the atmosphere model WRF, the wave model SWAN, and the sediment capabilities of the Community Sediment Transport Model. In this first development phase only the ocean model ROMS (Regional Ocean Modelling System) is running operationally.

ROMS (Shchepetkin and McWilliams, 2005) is a three-dimensional, free-surface, terrain-flowing numerical model that solve the Reynolds-averaged Navier-Stokes equations using the hydrostatic and Boussinesq assumptions. It has been specially designed for accurate simulations of regional oceanic systems for a diverse range of applications (Warner et al., 2005; Wilkin et al., 2005). The algorithms that comprise ROMS computational nonlinear kernel are described in detail in Shchepetkin and McWilliams (2003, 2005), and the tangent linear and adjoint kernels and platforms are described in Moore et al. (2004). ROMS includes accurate and efficient physical and numerical algorithms and several coupled models for biogeochemical, bio-optical, sediment, and sea ice applications. It also includes several vertical mixing schemes (Warner et al., 2005), multiple levels of nesting and composed grids and wetting and drying capabilities.

The HR-OSM is based on the oil spill model TESEO developed by the Environmental Hydraulics Institute of the University of Cantabria (IH Cantabria) (Abascal et al., 2007). TESEO is composed of a transport module and a weathering module that simulate the processes governing the evolution and behaviour of oil spilled on the marine water. Oil spill transport is described by tracking numerical particles equivalent to the oil slicks by means of the two-dimensional Lagrangian approximation. At each time step, the new position of the particles is computed by the superposition of the transports induced by the mean flow, tides, wind/waves and turbulent dispersion. The processes encoded in TESEO to describe physical and chemical weathering processes of spilled oil include: spreading, evaporation, emulsification and the change of rheological properties. The model also takes into account oil beaching.

To get the desirable resolution (46 m) two embedded grids with two-way nesting are used. The HR bathymetry used to build these grids is the result of the assemblage of well-known regional bathymetric databases: Etopo 1' and Admiralty charts, local bathymetric databases provided by the Santander Port Authority and the Topography Department of the University of Cantabria and two bathymetric surveys carried out in the framework of SPRES project. Figure 3 shows the set of grids used and Table 1 summarizes their main characteristics.

In relation to the vertical distribution, 5 vertical levels in σ coordinates have been defined. ROMS has a generalized vertical, terrain-following, σ coordinate system. In this application the transformation equation chosen was the first one described in Shchepetkin and McWilliams (2005) and the stretching function is the one developed by Song and Haidvogel (1994).

Once the model domain was defined, the downscaling procedure consisted of designing a pre-processing scheme to prepare the boundary and initial conditions (from MyOcean-IBI, http://www.myocean.eu.) and the atmospheric forcing (from AEMET, Spanish MetOffice, http://www.aemet.es.) to feed the COAWST model (see Figure 2) .

Initial and boundary conditions for the HR-OOM are provided by MyOcean-IBI (Iberian-Biscay-Irish) Physical Forecasting System. This system is based on a NEMO model application run at 1/36º horizontal resolution. The System is daily run by Puertos del Estado and Mercator-Ocean and provides 5-day hydrodynamic forecast (+ 1 day hindcast as best estimate) including high frequency processes to characterize regional scale marine processes (i.e. tidal forcing, surges, high frequency atmospheric forcing, fresh water river discharge, etc.). It provides 3D daily means field of water temperature, salinity, sea surface height, zonal velocity and meridional velocity, as well as hourly means of surface fields (sea height, water temperature and currents). Figure 4 shows the IBI domain.

The initial conditions consist of a three dimensional structure of salinity and water temperature, depth averaged currents, and free surface elevation for both grids. These initial data can come from: (1) the HR-OOS previous (day-before) simulation; (2) in case of any failure, from IBI-MyOcean interpolated data to the two nested grid.

To obtain the boundary conditions, the HR-OOS interpolates the following IBI-MyOcean variables to the intermediate grid nodes:

• Vertical profiles of water temperature and salinity (daily resolution)

• Depth averaged currents (hourly resolution)

• Baroclinic currents (daily resolution)

• Free surface elevation (hourly resolution)

COAWST model has different options to include these boundary conditions. After several calibration tests, the following options were activated:

• Temperature and salinity: Flather condition

• Depth averaged currents: Chapman condition

• Baroclinic currents: Radiation condition

• Free surface elevation: Radiation and nudging condition

The atmospheric forcings are provided by AEMET, the Spanish MetOffice (http://www.aemet.es). AEMET, through the forecasting system HIRLAM (High Resolution Limited Area Model) which runs at 1/7º in Atlantic area, provides every 6 hours, 72-hours atmospheric forecast (hourly data). The atmospheric variables used by HR-OOS are northward and eastward wind at 10 m over the sea surface, sea level pressure, SLP, and heat and freshwater fluxes, all of them with hourly resolution. Some pre-processing is needed to obtain the final atmospheric input file (e.g. extraction of the land values and extrapolation of the sea values to land).

ROMS model has different options to include heat and freshwater fluxes. The one selected in the HR-OOS consists of activating the bulk fluxes module (Fairall et at., 1996a, 1996b), which, although is more time-consuming, gives the best results.

The HR-OOS is run daily by IH-Cantabria and provides 2-day hydrodynamic forecast (sea level, currents, temperature and salinity) and 1-day of hindcast as best estimate. Furthermore, a web service has been developed that manages the operational system and allows the user to run hypothetical as well as real oil spill trajectories using the daily forecast of wind and the high resolution ocean variables as input for the TESEO model. Figures 5 and 6 show and example of the HR-OOS interface.

To improve the quality of the simulations two field campaigns (winter/summer) to calibrate/validate the hydrodynamic model COAWST have been performed. Moreover, in order to improve the weathering formulations included in the oil spill model, TESEO, laboratory test are being conducted by the Centre of Documentation, Research and Experimentation on Accidental Water Pollution (Cedre).

The HR-OOM has been calibrated and validated using data from a permanent tidal gauge located inside the Bay (managed by Puertos del Estado) and two field campaigns carried out in summer of 2012 and winter of 2013 (see Figure 7). The summer campaign lasted 19 days, from 18th July to 6th August, and the winter campaign lasted 16 days, from 7th March to 22nd March.

The hydrodynamic model calibration consists of quantifying the optimal hydrodynamic model parameters, comparing modelled results with real observed values by using statistical descriptors. In this case, the unknowns were the hydrodynamic model parameters (bottom drag coefficient C_D and the horizontal eddy viscosity, k) and the turbulence closure model (TCM). We have used several statistical descriptors such as the root mean square error (RMSE), the correlation index (^ρ) and the skill index (^skill), which was finally selected to evaluate the model performance:

where ^x_obs are the instrumental data, ^x_model are the model results and ⁿ are the length of both datasets. A perfect agreement between model results and observations involves ^skill parameter values equal to one (^skill=1), whereas a complete disagreement is associated with zero (^skill=0). Therefore, the objective is to find the best combination of (CD, k, TCM) that maximizes the skill index:

where ^nc is the number of campaigns (studied periods), ^nm_c is the number of measured magnitudes, and ^ns_c,m is the number of sensor of each campaign that measured the studied magnitudes.

Since the hydrodynamic models are very time consuming, this ideal calibration methodology only can be applied in a discrete approximation: the decision variables, ^C_D and ^k are only allowed taking several values, 3 values each; and the decision variable TCM has only 2 options (constant, or ^k−ε model). According to the global descriptors obtained by solving Eq. (2), the best result were achieved with the following parameters: ^C_D=0.00392, ^k=0.05 and ^k−ε as TCM.

It is noteworthy that, after doing a detailed study of the statistical descriptor values per each instrument (i.e. per each location), the results are very dependent on the location of the sensor. For example, in some inner parts of the Santander Bay (the narrowest areas) the horizontal resolution is insufficient and because of that their statistical descriptors are the lowest. These less accurate descriptors mask some excellent ones in other locations which are more interesting for the target of the project.

Figure 8 shows an example of the validation using the best parameter estimates for water level and horizontal currents. Water level results obtained with the calibrated HR-OOM are well suited to the real conditions. With regard to the current velocity results, the model gives satisfactory results obtaining ^RMSE less than 0.12 m/s in all the cases.

Although the HR-OOS results reproduce adequately the Bay hydrodynamics, future improvements, in the case of having more computational resources, could be to increase the vertical and horizontal resolution in order to achieve better results mainly in the narrowest channels of the Bay. It could be possible by nesting 3 computational grids instead of 2 grids, or by using curvilinear grids with more than 10 vertical layers.

As was mentioned, the HR-OSM is made up by the oil spill model TESEO (Abascal et al., 2007). The transport module of TESEO was calibrated and validated in previous works using data from drifting buoys (Abascal et al., 2007, 2009a, 2009b). Specifically, buoys deployed in the Bay of Biscay during the Prestige accident were used to calibrate the model and to obtain the optimal coefficients for the study area (Abascal et al., 2009a).

The weathering module included in TESEO simulates spreading, evaporation and emulsification using well know formulations (e.g. Fay (1969, 1971); Lehr et al. (1984); Stiver and MacKay (1984); MacKay (1980)). In order to adjust the weathering formulations parameters, pilot scale experiments were conducted by Cedre in the framework of the SPRES project. Two different refined products were chosen in this study. A light product, a diesel oil, and a heavier oil, an IFO 220 (Intermediate Fuel Oil characterized by a viscosity of 220 cSt at 50°C). Tests were performed considering various environmental conditions: 3 temperatures (10, 15 and 20°C), 3 salinities (5, 20 and 35 g/L), and 3 different states of energy (low, average and high, the highest level corresponding to about sea state 3). The tests lasted 3 days (diesel oil) or 7 days (heavy fuel oil) in Cedre's hydraulic canal (the Polludrome), which was used with the water being continuously circulated to simulate open sea conditions (see Figure 9). During the tests, surface oil samples were taken periodically to determine the oil characteristics (viscosity, density, water content and evaporation rate) and the concentration of PAHs dissolved in the water column (only in the case of the IFO 220). Conditions of weathering experiments were in agreement with those described previously (Guyomarch et al., 2012).

The model parameters to be calibrated are the empirical coefficients for each weathering process considered in the model: evaporation, emulsification and the change in density and viscosity. The parameters used to calibrate the weathering processes are: d (evaporation); Kw (emulsification); CDE (change in density); CT, CE, CV and CM (change in viscosity). The statistical descriptors used to compare model results with experimental values are the root mean square error (RMSE), the correlation index (r) and the residual scatter index (RSI).

The shuffled complex evolution method (Duan et al., 1992) developed by the University of Arizona (SCE-UA) has been applied to estimate the optimal coefficients of the oil spill model. The SCE-UA method is based on a synthesis of four concepts: (1) a combination of deterministic and probabilistic approaches; (2) a systematic evolution of a “complex” of points spanning the parameter space in the direction of global improvement, (3) competitive evolution; (4) complex shuffling. The synthesis of these elements makes the SCE-UA method effective and robust, and also flexible and efficient.

Following this methodology, the oil spill model calibration was formulated as an optimization problem where an objective function J has to be minimized:

where FT is the number of experimental tests, WP is the number of weathering processes considered, and θ=(d, Kw, CDE, CT, CE, CV, CM) is the vector of unknown parameters.

Eq. (3) represents the mean of the values of the dimensionless standardized statistical index (residual scatter index (RSI)) obtained by comparing the predicted weathering processes of the model (evaporation, emulsification and changes of physicochemical properties), with the experimental results.

According to the global descriptors obtained by applying Eq. (3), the coefficients minimizing the objective function were: d= 2.034E-3, Kw= 1.085E-06, CDE= 0.18006, CT= 1.00003, CE= 0.6504, CV= 0.23812, CM=15.44. The value of the global scatter index obtained was 0.62402.

After a preliminary analysis of the calibration and validation results, we obtained the best results for the IFO 220. And specifically, for the evaporation with RMSE= 0 and RSI=0. We found that the poor results obtained with the Diesel oil are due to the overestimation of the initial area of the simulated oil slick. Ongoing analysis of these results is being conducted. As an example of the results, Figure 10 shows the comparison between experimental data and model results for one of the test carried out with the heavier oil, IFO 220. We can observe that although the model reproduces quite well the time evolution of the evaporation, emulsification and viscosity, density curves does not show a good agreement. These poor results for the density are consistent for all the IFO 220 tests. More research is required to explain this behavior. Overestimation of the initial area of the simulated oil slick and procedures for measuring density at laboratory are factors that will be investigated in future studies.

This work has been co-funded by project SPRES in the framework of the EU-Atlantic Area Programme and by project PLVMA (TRA2011-28900-Spanish Ministry of Economy and Competitiveness). Authors B.P. and M.C. would like to thank the Spanish Ministry for its support within the FPI Program. The authors would like to thank AEMET, AZTI-Tecnalia, Cedre, MyOcean and Puertos del Estado for the data and collaboration provided in this project.