MEDSLIK-II, a Lagrangian marine surface oil spill model for short-term forecasting – Part 2: Numerical simulations and validations

In this paper we use MEDSLIK-II, a Lagrangian marine surface oil spill model described in Part 1 (De Dominicis et al. , 2013), to simulate oil slick transport and transformation processes for realistic oceanic cases, where satellite or drifting buoys data are available for verification. The model is coupled with operational oceanographic currents, atmospheric analyses winds and remote sensing data for initialization. The sensitivity of the oil spill simulations to several model parameterizations is analyzed and the results are validated using surface drifters, SAR (synthetic aperture radar) and optical satellite images in different regions of the Mediterranean Sea. It is found that the forecast skill of Lagrangian trajectories largely depends on the accuracy of the Eulerian ocean currents: the operational models give useful estimates of currents, but high-frequency (hourly) and high-spatial resolution is required, and the Stokes drift velocity has to be added, especially in coastal areas. From a numerical point of view, it is found that a realistic oil concentration reconstruction is obtained using an oil tracer grid resolution of about 100 m, with at least 100 000 Lagrangian particles. Moreover, sensitivity experiments to uncertain model parameters show that the knowledge of oil type and slick thickness are, among all the others, key model parameters affecting the simulation results. Considering acceptable for the simulated trajectories a maximum spatial error of the order of three times the horizontal resolution of the Eulerian ocean currents, the predictability skill for particle trajectories is from 1 to 2.5 days depending on the specific current regime. This suggests that re-initialization of the simulations is required every day.


Introduction
MEDSLIK-II has been designed to provide timely information on oil spill advection-diffusion and weathering after a surface oil spill release.This model has the potential to become part of an operational detection-prediction system using observed oil slicks as initial conditions and prediction of their movement and transformation to guide oil spill response activities.
The MEDSLIK-II model described in Part 1 of this paper (De Dominicis et al., 2013) is capable of predicting physical changes of a surface oil spill and uses a Lagrangian particle representation for the transport and diffusion processes.MEDSLIK-II has been coupled to operational Ocean General Circulation Model (OGCM) outputs that provide analyses and forecasts for the deterministic components of the particle trajectory equations (Coppini et al., 2011;Zodiatis et al., 2012).Moreover, atmospheric forecast models provide surface winds for the transformation process, the surface current corrections and the computation of wind waves affecting the transport.Additionally the model can be initialized using the slick position and slick shape provided by satellite systems, Published by Copernicus Publications on behalf of the European Geosciences Union.
both SAR (synthetic aperture radar) and optical images.Such satellite images of surface oil slick have already been combined with Lagrangian trajectory models in rapid response to the Deepwater Horizon incident (Liu et al., 2011a, c).
Validation of oil spill models is usually carried out comparing surface buoy drifter trajectories with modeled trajectories.The papers of Reed et al. (1994) and Al-Rabeh et al. (2000) showed a qualitative comparison between drifting buoy trajectories and modeled trajectories.In more recent papers, quantitative metrics, based on the separation distance between modeled and observed trajectories, are presented (Price et al., 2006;Barron et al., 2007;Caballero et al., 2008;Sotillo et al., 2008;Huntley et al., 2011;Cucco et al., 2012).A skill score, based on separation distance normalized by the trajectory length, has been recently proposed by Liu and Weisberg (2011).This new metric has already been used by Röhrs et al. (2012) and Ivichev et al. (2012) to evaluate their model performances.Oil spill models' forecasting accuracy can be also evaluated by comparing the model results to remote sensing observations (Carracedo et al., 2006;Coppini et al., 2011;Berry et al., 2012;Mariano et al., 2011;Liu et al., 2011c), although it is difficult to have oil slick time series for long periods after the first observation, due to the long revisit time for satellites.Between those studies, the pioneering study of Reed et al. (1994) combined the drifters and remote sensing observations with chemical samplings.However, no study has been done up to now that systematically evaluates the predictability time of the oil spill evolution and the model sensitivity to many of the uncertain model parameters, such as the oil properties and the type of current information given for the transport of the oil.
In this paper, we illustrate three groups of experiments in order to understand the sensitivity of oil slick simulations to different model assumptions and validate the results with in situ and satellite data.First, we focus on the model skill in simulating single drifter trajectories as a function of the space and timescales of the Eulerian current field, the impact of local wind and the wave-induced velocity correction terms.Secondly, we show the sensitivity of the simulated oil slick, initialized from satellite observations, to uncertain oil input properties, such as oil type, slick thickness and age.Thirdly, sensitivity tests to the number of Lagrangian particles and tracer grid resolution is presented.
All these experiments are compared to observed data and the degree of predictability of the trajectories is evaluated in terms of: RMSE (root mean square error) between observed and simulated particle trajectories as a function of model parameters and skill score proposed by Liu and Weisberg (2011).This will set the limit of predictability of oil spill evolution as a function of the Eulerian input fields' horizontal resolution.
The manuscript is organized as follows: Sect. 2 overviews the model equations and parameters already presented in Part 1 of this paper (De Dominicis et al., 2013), the coupling with OGCM and atmospheric fields, the oil spill model pa-rameters and the initialization procedures; Sect. 3 presents the drifter data and the satellite images used to validate the model together with the methodology to quantitative assess the Lagrangian forecasts; Sect. 4 presents the results of the validation experiments; Sect. 5 offers the conclusions.

MEDSLIK-II model setup
This section describes the main equations of MEDSLIK-II model and the oil spill parameter values chosen in our simulations, the description of the ancillary environmental fields needed as input to the oil spill model and the algorithms for initialization of MEDSLIK-II from observed satellite images.

MEDSLIK-II model equations
The MEDSLIK-II model equations, presented in Part 1, are overviewed in this section.The oil spill model state variables are reproduced in Table 1 from Part 1 of this paper (De Dominicis et al., 2013).Three kinds of state variables are defined in the model: the concentrations called structural state variables, the oil slick and particle state variables that are used to simulate weathering and transport-diffusion processes respectively and to reconstruct the concentrations.
MEDSLIK-II only allows for a simulation of the evolution of a surface oil volume release, indicated by V S .Using Mackay's approach (Mackay et al., 1979(Mackay et al., , 1980)), the oil slick is subdivided into thin (sheen) and thick parts, described by the oil slick state variables: the volumes of the thick and thin parts of the slick, V TK and V TN respectively, the thick and thin slick areas, A TK and A TN and the thick and thin thicknesses, T TK and T TN .The oil slick variables are then written as The thin and thick area initial values are taken from the known initial surface amount of oil released, V S (x C , t 0 ), using the F parameter, which is the area ratio of the two slick parts, A TK and A TN , and assuming the initial values for the thicknesses where t 0 is the initial time and x C is the slick's central geographical position.
where the suffixes indicate evaporation (E), dispersion (D) and spreading (S), and all the slick variables related to volume are defined at the slick centre.
The initial surface oil volume is broken into N constituent particles characterized by particle variables, which are the position vector and the particle volumes υ(n k , t), where n k is the particle identification number.Each particle is characterized by a status index (see Table 1 in Part 1) which indicates if the particle is at the surface, in the subsurface or on the coast.The variation of the oil particle volumes, υ(n k , t), are linked to the weathered oil slick volumes of Eqs. ( 6) and ( 7) using empirical relationships described in detail in Part 1.
The advection-diffusion processes are solved using the N Lagrangian particles and the prognostic equations for their displacements are where σ = 0, 1 is the particle index that describes if the particle is respectively at the surface or dispersed, U C is the current velocity term, U W is the local wind velocity correction term, U S is the wave-induced current term (Stokes drift velocity), K is the turbulent diffusion diagonal tensor and Z is a vector of independent random numbers used to model the Brownian random walk processes chosen for the parametrization of turbulent diffusion.The turbulent diffusion is considered to be horizontally isotropic and the three diagonal components of K are indicated by K h , K h , K v .The transformation of the particles from the surface to the subsurface status is only due to the dispersion processes, as described in Part 1. Once the particle is in the subsurface, at a particular depth z k , it is horizontally dispersed by the correspondent horizontal velocity field at that depth (Eq. 9 for σ = 1).If U C is the output of a baroclinic, wind-driven oceanographic model, the currents will contain a satisfactory representation of surface ageostrophic currents in the surface and deep layers of the water column.For surface currents in particular, the U W term can be neglected.The surface wind term in fact is necessary when U C is estimated from climatological data using the geostrophic assumption (Al-Rabeh et al., 2000) or when the oceanographic models do not resolve accurately the upper ocean dynamics.In these cases, U W can be considered as a correction term accounting for uncertainty and unresolved processes in U C at the surface.Furthermore, U S accounts for the presence of surface wave current drift: in MEDSLIK-II it is introduced using an analytical formulation that depends on wind amplitude, as explained in Appendix C of Part 1.In the future, swell and other wave processes should be considered using the Stokes drift coming from a numerical wave model.
Finally, the surface (C S ), dispersed (C D ), and on-coast oil concentrations (C C ) are reconstructed using an oil tracer 2-D coordinate system (x T , y T ) with an uniform horizontal resolution (δx T δy T ) as where ρ is the oil density, C S and C D are expressed in units of kg m −2 and C C (L i , t) as kg m −1 , I S and I D are the particles on the surface and dispersed and I C is the set of particles beached on the coastal segment L i that has a length δL i , discussed in details in Part 1.
The minimum/maximum number of particles used to represent the miminum/maximum concentrations (C S min and C S max ) for any given initial release V S can be calculated as where N S is the number of sub-spills in which the oil volume is subdivided for a continuous time spill (see Part 1).

Oil spill model parameters
As described in Part 1 of this paper, many empirical model parameters and parametrizations are considered in MEDSLIK-II and they have been listed in Table 2 of Part 1, together with their nominal values from published literature.
In this paper, we left all parameters equal to their nominal values except for the number of initial particles, N, the tracer grid cell size, (δx T , δy T ), the thickness of the thin slick (Eq.5) and the horizontal diffusivity coefficient, K h .In the simulation experiments of single drifter trajectories, see Sect.4.1, the horizontal diffusivity coefficient K h is set to zero, while simulating an oil slick from satellite, see Sect.4.2, K h has been set to 2 m 2 s −1 in the range 1-100 m 2 s −1 indicated by ASCE (1996) andDe Dominicis et al. (2012).

Ancillary ocean and atmospheric fields
MEDSLIK-II requires data on wind forcing, sea surface temperature and sea currents in order to compute the transport (Eq.9) and transformation processes (Eqs.6-7).Wind forcing, i.e., the wind velocity components at 10 m above the sea  Pinardi et al., 2003;Pinardi and Coppini, 2010), the Adriatic Forecasting System (AFS, Guarnieri et al., 2013) and the IRENOM relocatable model, explained below.
The MFS system (Tonani et al., 2008) is composed of an OGCM at 6.5 km horizontal resolution and 72 vertical levels (Oddo et al., 2009) and an assimilation scheme (Dobricic et al., 2008) which corrects the model's initial guess with all the available in situ and satellite observations, producing analyses that are initial conditions for ten days ocean current forecasts.In this paper, U C comes from daily and hourly mean analyses in order to eliminate the additional uncertainty connected with forecasts for both atmospheric and oceanographic input data.
The MFS basin-scale output provides initial and lateral boundary conditions for high-resolution models, thereby resolving the coastal dynamics better.AFS is one of the nested models with a horizontal grid resolution of 1/45 • (approximately 2.2 km) and 31 vertical sigma levels, and it also considers tidal motion (Guarnieri et al., 2013).AFS produces simulations and forecasts, which are provided as daily and hourly mean outputs.
The IRENOM relocatable model has been designed in order to provide high/very high time and space resolution forecasts starting from operational large-scale circulation models, such as MFS (Fabbroni, 2009).The hydrodynamics model core is based on the Harvard Ocean Prediction System (Robinson, 1999) and in this work IRENOM has been implemented with 3 km horizontal resolution, starting from approximately 6.5 km resolution MFS fields, and 40 vertical sigma layers.Initial and lateral boundary conditions are obtained from MFS.The atmospheric forcing is interactively computed using the ECMWF operational products.The model outputs are daily and hourly simulation fields.

Oil slick initialization from satellite images
The data required to define the oil slick initial condition are the total surface volume released, the geographic location, the time, the oil type, the area covered by the slick and its thickness, as well as the age of the oil slick from the initial release into the sea.
Most of this information can be estimated from satellite sensors.Synthetic aperture radar (SAR) and optical images can provide, as satellite image post-processing products, the area covered by the slick and the slick contour coordinates (Trivero et al., 2001;Nirchio et al., 2007Nirchio et al., , 2010)).The total oil slick area at the initial time t 0 is the sum of the thick and thin parts, A(t 0 ) = A TN (t 0 ) + A TK (t 0 ).Thus, by combining Eqs. ( 4)-( 5), the initial surface oil volume release can be calculated as The information on the area ratio F and thicknesses are normally unknown and have to be hypothesized.In our study, F and T TK are fixed and they are taken from the standard values listed in Table 2 of Part 1, while T TN will be varied between 1 and 10 µm.The N Lagrangian particles initial positions, x k (t 0 ) within the slick contour, are determined using the method described in the Appendix A.
A novel feature of MEDSLIK-II is its ability to initialize, within the satellite image slick area, the slick and particle state variables, such as the volume of the thick and thin slicks, V TK (t 0 ) and V TN (t 0 ), and oil particle volume υ(n k , t 0 ).In order to calculate these variables, the age of the slick has to be hypothesized.A simulation, with weathering processes only, is performed for a time period equal to the assumed slick age (see Fig. 1).During this phase the particles do not change their initial position, but the slick and particle state variables are evolved using Eqs.( 7)-( 8), starting at a time equal to the time at which the spill has been observed by satellites minus the assumed slick age.

Verification drifters and satellite data
Verification of oil spill forecasting is both a crucial issue and a difficult task to perform.The main reason for this is the lack of oil slick time series for long periods after the first observation, due to the long revisit time for satellites and the scarcity of in situ data.In this paper, we will use both drifters trajectories data and satellite imagery to validate MEDSLIK-II simulations.
The CODE drifters used in this paper were released in the Ligurian Sea in 2007 (Poulain et al., 2011) and will be used here to study the impact of U C horizontal resolution and depth.
The IESM-PTR buoys are independent floating ARGOS buoys and are parallelepipeds measuring 30 cm in height (30×10×10 cm) and they are designed as oil-spill-following surface drifters.The IESM-PTR drifters were deployed south of Nice in autumn 2007 (Brostrom et al., 2008) and were used to show the effects of wind corrections and Stokes drift, U W and U S respectively, in Eq. ( 9).
The newest drifters are the OSDs (oil spill drifters), which are 32 cm diameter cylinders with a low degree of submergence, designed to follow oil spills and surface pollution.Dev., 6, 1871-1888, 2013 www.geosci-model-dev.net/6/1871/2013/Inizializa'on of slick state variables WEATHERING PROCESSES calcula3on (evapora3on, dispersion, emulsifica3on) considering the wind and SST in the area where the spill is observed.

Simula'on of the spill advec'on, diffusion and weathering
Full MEDSLIK--II dynamics considering WEATHERING and ADVECTION/DIFFUSION processes as explained in Part I.The comparison between observed and simulated drifter trajectories will be evaluated by the RMSE, calculated from the separation distance between the observed and the simulated trajectories as a function of the simulation time:

Observa(on (me + Simula(on length
where d i is the distance at the selected time t i , after a reference time t 0 , between the simulated drifter position, x s , and the observed positions, x o , and S is the total number of simulations using the same model parameters. The acceptable maximum separation between observed and modeled trajectory depends on the particular model application.An error between 7 and 19 km would allow the use of the model forecasts in situations of rapid response, such as oil spills and search and rescue operations.We should have in mind that the oil spill model results should be used to deploy booms, to place skimmers, to protect a particular piece of coast, to intervene with airplanes or vessels.Furthermore, it is common wisdom that in finite difference models eight grid points are required (Haidvogel and Beckmann, 1999) to resolve a structure.Thus, taking a quite conservative limit, we can consider acceptable a spatial error of the simulated trajectories of the order of three-four times the horizontal resolution of the Eulerian ocean currents.
The MEDSLIK-II model performances will be also compared with the state-of-the-art assessment of Lagrangian predictive skill.Price et al. (2006) found that separation distance between modeled and drifters trajectories is 78 km after 3 days, and 229 km after 20 days.RMSE over the integration time, estimated by Barron et al. (2007), ranges from 10 to 25 km after 1 day and from 50 to 150 km after 7 days.Caballero et al. (2008) showed that after 3 days observed and modeled trajectories separate by 23 km, while after 7 days the separation increases to 46 km.Sotillo et al. (2008) obtained a mean RMSE, among the 13 days of simulation, of 5 km.Liu and Weisberg (2011) found a separation distance between 13 and 34 km after 1 day, and between 58 and 177 km after 5 days.Huntley et al. ( 2011) indicated that the model trajec-tories separate from observations by roughly 15 km after the first day on average.Cucco et al. (2012) mean separation distance is 4 km after 2 days.
In addition to the RMSE calculation (see Eq. 13), the new skill score proposed by Liu and Weisberg (2011) has been used to further evaluate the performance of the modeled trajectories.A non-dimensional index is defined as an average of the separation distances weighted by the lengths of the observed trajectories: where S and d i have been already defined in Eq. ( 13) and l oi is the length of the observed trajectory at the corresponding time, t i , after a reference time t 0 .Such weighted average tends to reduce the evaluation errors that may rise using only the purely Lagrangian separation distance.The s index can be used to define a model skill score: where n is a tolerance threshold.In this work, as suggested by Liu and Weisberg (2011), we used n = 1, this corresponds to a criterion that cumulative separation distance should not be larger than the associated cumulative length of the drifter trajectory.The higher the ss value, the better the performance, with ss = 1 implying a perfect fit between observation and simulation and with ss = 0 indicating the model simulations have no skill.This skill score may have some limitations in case of very weak currents and hence small cumulative distances, that may imply a very large value of s and very low skill score ss.These limitations may be overcome by setting a proper tolerance threshold, n, as suggested by Liu and Weisberg (2011).Finally, the model has also been validated using remote sensing data from satellite images obtained using both synthetic aperture radar (SAR) (Trivero et al., 1998;Fiscella et al., 2000;Trivero et al., 2001;Nirchio et al., 2005Nirchio et al., , 2007Nirchio et al., , 2010) ) and MODIS optical sensors (Hu et al., 2003(Hu et al., , 2009)).The satellite data allowed for a study of the importance of shape initialization, the sensitivity to oil slick input properties (thickness, oil type and age) and to the number of constituent particles. of spatial and temporal current resolution in the U C term in Eq. ( 9), the local wind correction term U W and the Stokes drift U S , which are written as (see Part 1) where (W x , W y ) are the wind velocity components at 10 m, ϑ = arctg W x W y is the wind direction and D S is the Stokes drift velocity intensity in the direction of the wave propagation at the surface, defined as where ω is angular frequency, k is wave-number, and S(ω) is wave spectrum.The turbulent diffusion coefficient K h was set to zero in all the experiments that are described in this section.
The oceanographic fields (hourly and daily currents) are obtained from the operational MFS OGCM and the nested high-resolution IRENOM.The winds are from ECMWF analyses at 6 h time resolution.The total length of the simulation is 3 days, and all simulated and real drifters were launched at the same time on 14 May 2007 at 15:00 UTC.
Figure 2 shows the real drifter tracks (black lines) for three days and the simulated MEDSLIK-II trajectories for the five experiments of Table 1.The trajectories obtained using the daily MFS surface fields are not capable of reproducing the correct drifter direction.When high time frequency MFS fields (CURR-EXP2, Table 1) are used, the simulated drifters have the correct direction but are much too slow as compared to reality.When higher horizontal resolution IRENOM hourly fields are used (CURR-EXP3, Table 1), the trajectories are in better agreement with the observations.We therefore conclude that hourly and relatively high resolution currents are needed to reproduce the trajectories of observed drifters.
This is confirmed by the RMSE curves shown in Fig. 3a.Indeed, using daily currents (CURR-EXP1) the distance error is always higher than the one of CURR-EXP2 and CURR-EXP3.The best results are shown by CURR-EXP3: for the first 24 h of simulation the distance error calculated using Eq. ( 17) is of the order of the hydrodynamic model resolution (IRENOM, 3 km), after 48 h the error remains within two times the model resolution (6 km) and after 60 h the error is three times the model resolution (9 km).These separation distances are lower than the values generally obtained in previous works (see Sect. 3).The RMSE of CURR-EXP2 (MFS, 6.5 km) confirms the same behaviour observed for CURR-EXP3, although it is slightly worse at all times, making evident the fact that by increasing horizontal model resolution we can improve the predictability time for particle trajectories.Considering acceptable a spatial error of the simulated trajectories of the order of three times the horizontal resolution of the Eulerian ocean currents, the predictability time for this case is 2.5 days.A restart of the simulation should be required every day to maintain the distance error of the same order of the model resolution (not shown).The skill scores in Fig. 3b confirm that the best results are obtained by CURR-EXP3.The CURR-EXP3 skill score reaches 0.86 after 1 day of simulation, it remains constant during the second day of simulation and it starts to slowly decrease to 0.81 at the end of the simulation.While the maximum skill score obtained in CURR-EXP2 is 0.79 at the beginning of the simulation and then it decreases up to 0.59 at the end of the simulation.
In the CURR-EXP4 and CURR-EXP5 simulations (see Table 1), we test the impact of using the surface currents provided by the MFS OGCM versus the 30 m currents, assumed to be the geostrophic components, with the addition of a 3 % wind velocity Ekman current correction estimate (see Eq. 16), using a wind angle equal to 0 • and 25 • (the wind angle range indicated by Al-Rabeh, 1994).This is to correct for OGCM inaccuracies in the simulation of the Ekman dynamics.In Fig. 2 we can observe that this correction and composition of the surface currents does not give as accurate a representation as the direct MFS surface fields, as confirmed by the RMSE trends shown in Fig. 3a and by the skill scores in Fig. 3b.A similar result was found by Al-Rabeh et al. (2000) but it is difficult to generalize since we argue that this depends on the specific Ekman process occurring at the surface and the vertical resolution of the OGCM.Ekman corrections should be carefully tested for coastal currents where the major forcings are local bathymetries and coastlines.
Other model sensitivity experiments were carried out for the IESM-PTR drifters, taking the currents from MFS hourly analyses and winds from ECMWF 6-hourly analyses.The simulations were carried out applying different wind and Stokes drift corrections as described in Tables 2 and 3.In Fig. 4 the observed drifters were released on 10 October 2007, while the numerical numerical drifters were launched on 14 October 2007 at 1 a.m. and followed up to 22 October 2007.We want to show first this case because we have an interesting positive impact of the wind correction here even if for a particular case.Figure 4 shows that the observed drifters move parallel to the coasts between 5 and 7 • E and between 4 and 5 • E they translate offshore, probably under the influence of winds.We note that using the wind correction (WIND-EXP4) we reproduce the observed drifter movement offshore and southward, which is not reproducible using only the MFS currents.As shown in Fig. 6b, the skill score trend of WIND-EXP4 is always above the skill scores of the other experiments and it reaches the maximum value of 0.85.The distance error (see Fig. 6a) is of the order of three times the model resolution (MFS, 6.5 km) after 24 h.Although higher than in the previous experiments (CURR-EXPs), these separation distances are still lower or   of the same order of the values generally obtained in previous works (see Sect. 3).Thus, we argue that the predictability skill for the particle trajectories in this current regime is 1 day.In Fig. 5 the simulation is then re-initialized every day, showing the capability of the model to reproduce the entire drifters trajectories (8 days) maintaining the error within three times the model resolution (in WIND-EXP3, WIND-EXP4, WIND-EXP5, see Fig. 6c).The separation distance are significantly lower by using the re-initialization and this is confirmed by the higher skill scores shown in Fig. 6d.In WIND-EXP4 and WIND-EXP5, the skill scores reach the value of 0.7 after 24 h and grow up to 0.9 after 8 days.

In order to understand what the wind correction means in the experiments of
experiments, SD-EXP1 and SD-EXP2, listed in Table 3.We note that, from Eqs. ( 13) and ( 16), the wind correction with an angle of 0 • is analogous to the Stokes drift correction parameterization except for the fact that the correction amplitude is determined by a fixed parametrization for the Stokes drift while it is arbitrary in the wind case.In Fig. 4 we show that the Stokes drift correction (SD-EXP2) is less effective than the wind correction to reproduce the observed trajectory.
We therefore argue that in this case the wind correction has parameterized the direct effect of wind drag on the IESM-PTR buoy rather than accounting for missing wave induced surface drift.The effect of Stokes drift correction was also studied using the OSD drifter in the coastal area near Cesenatico (northern Adriatic Sea).The drifter was launched on 21 July 2009 at 09:40 UTC. and was at sea for nearly a week.The simulations were carried out using the hourly current fields provided by the AFS model and the ECMWF 6-hourly wind fields.The different experiments are described in Table 3 and the results are shown in Fig. 7.
The simulated drifters were deployed daily and simulations lasted 24 h, starting from a simulation on 21 July 2009 at 09:40 UTC and lasting 15 h.As shown in Fig. 7 the model once again appears to underestimate the current intensity in the northward direction, with the result that the inertial oscillation loops are tighter than they are in the observations.Using the re-initialization, the error is maintained with five times the AFS model resolution (2.2 km) for 6 days of simulation.However, if we consider again as maximum acceptable error three times the model resolution, the predictability skill is now only 16-18 h after each re-initialization, as shown in Fig. 8a.From Figs. 8a and 8b we argue that the simulated trajectories obtained by adding 1 % of the wind intensity of the current velocities, or considering the Stokes drift, are in better agreement with the observations than without the corrections.In this case, adding 1 % of the wind intensity and considering the Stokes drift gives almost identical results,   (Brostrom et al., 2008).Green lines are the trajectories simulated without any correction (WIND-EXP1/SD-EXP1); the red lines are the trajectories simulated using the MFS surface currents and a wind correction of 1 % (WIND-EXP2); the grey lines are the trajectories simulated using the MFS surface currents and a wind correction of 2 % (WIND-EXP3); the light blue lines are the trajectories simulated using the MFS surface currents and a wind correction of 3 % (WIND-EXP4); the blues lines are the trajectories simulated using the MFS currents at 30 m depth and a wind correction of 3 % (WIND-EXP5) and the pink lines are the trajectories simulated using the MFS surface currents and considering Stokes drift velocity (SD-EXP2).Note: in panel (b) the WIND-EXP1 trajectory is not visible because the simulated drifter arrived onto the coast after few hours of simulation.indicating that the wind correction can be interpreted as a parameterization of the wind/wave-induced current effects.

Geosci
The skill scores for SD-EXP4 and SD-EXP5 (see Fig. 8b) reach 0.8 after 24 h and grow up to 0.9 after 6 days of simulation, in agreement with the values find in the literature (Liu and Weisberg, 2011;Röhrs et al., 2012;Ivichev et al., 2012).
In order to validate the Stokes drift formulation described in Part 1 of this paper and the significant wave height calculations using the JONSWAP wave spectrum parameterization, the wave simulated by MEDSLIK-II has been compared with the data measured by a wave buoy during the period 21-27 July 2009.The buoy is located about 5.5 km off Cesenatico, over a depth of 10 m.Assuming that wave conditions offshore from Ravenna are comparable to those measured offshore from Cesenatico by the wave buoy, the comparison between measured and simulated waves by MEDSLIK-II is presented in Fig. 9.The waves simulated compare quite satisfactorily with observations, supporting the simplified calculation of the Stokes drift described in Part 1 of this paper.

Sensitivity of oil concentration to uncertain input parameters, number of particles and oil tracer grid resolution
In this section we validate the MEDSLIK-II simulation with SAR and optical satellite images.In the Mediterranean Sea it is very difficult to have subsequent satellite images over the same area in subsequent days, and even more difficult to have satellite images from the same research group or same sensor (SAR or optical images).This is why we are obliged to use different sources of satellite images, despite the inherent technological limitations.In Fig. 10 two slicks are shown: the first is observed by ASAR sensor (Trivero et al., 1998;Fiscella et al., 2000;Trivero et al., 2001;Nirchio et al., 2005Nirchio et al., , 2007Nirchio et al., , 2010) ) for 6 August 2008 and the other is observed by the optical sensor MODIS (Hu et al., 2003(Hu et al., , 2009) ) 25 h later.
We consider that the two images represent the evolution of the same oil slick, so we have both an initialization image and a verification one for the successive 25 h.The time of observation, the slick shape and area from the ASAR image are taken as initial surface slick variables for the simulation.
In Table 4 the parameters of the central simulation experiment are listed.Here no wind or Stokes drift corrections are used and two sets of sensitivity experiments were conducted: (1) to uncertain initial oil slick state variables such as oil age, oil type and thickness; (2) to number of constituent particles and tracer grid resolution.
The first set of experiments is described in Table 5.The oil slick age is taken to be varying between 0 and 24 h.We hypothesized an oil with an API of 22, which corresponds to an oil density of 0.92 tons m −3 and of 45, which corresponds to a lighter oil (density 0.804 tons m −3 ).The thin oil slick thickness, T TN , was changed between 1 µm and 10 µm.We assume an area factor, F , equal to 1000 and we consider the thickness of the thick part of the slick, T TK , equal to 0.1 mm (see nominal values in Part 1).We did not perform sensitivity experiments to T TK and F .Using Eq. ( 12), we obtained an initial surface oil volume, V S (t 0 ), equal to 764 m 3 when T TN = 10 µm and to 83.2 m 3 when T TN = 1 µm (A(t 0 ) is listed in Table 4).
Figure 11 shows the simulated oil slick location and concentration 25 hours after the initial detection of the oil.The modified shape of the slick is well captured by the model but the movement toward the north is probably too slow.No sensitivity to the age parametrization was observed in this case and in the following we will discuss only the experiments with age equal to 24 h.In Fig. 11 we compare the thinner slick and lighter oil simulation (ALGERIA-EXP3, Fig. 11a), with the thicker slick and heavier oil simulation (ALGERIA-EXP8, Fig. 11b).We can observe that after 25 h of simulation time, the oil concentration is almost zero for API 45 and T TN = 1 µm, whereas for API 22 and T TN = 10 µm the oil concentration is still high.Since the satellite optical image confirms the presence of the oil slick, we argue that ALGERIA-EXP8 is more realistic than  resolution to 150 m and the minimum detectable concentration limit to 0.1 tons km −2 and 30 tons km −2 , obtaining a maximum number of Lagrangian particles to be 300 000 and 1000 respectively.Figure 12a shows that using a coarse oil tracer grid, the concentration gradients are not correctly represented and the slick area is too large.Using a grid resolution of 50 m (Fig. 12b), we obtain a realistic estimate of the slick shape and area comparable to ALGERIA-EXP8 of Fig. 11b.However, the oil seems to be too uniformly distributed in the slick area.A smaller number of particles for the 150 m grid (Fig. 12c) generates a slick appearing as a large number of isolated and equal concentration oil slick sub-areas, while using a larger number of particles again a reasonable concentration is obtained (Fig. 12d).In conclusion we argue that an oil tracer grid of about 100 m and a number of particles around 100 000 gives the best results in terms of smoothness and consistency of the simulation with the area of a satellite detected oil slick.

Conclusions
In this paper we have shown an extensive calibration and validation of the MEDSLIK-II Lagrangian marine model for oil slicks described in detail in Part 1.The aim is to show the sensitivity of the oil slick simulations to choices of ancillary environmental conditions, advecting velocity parametrizations, oil slick parameters and number of Lagrangian particle and tracer grid resolution.In addition, the aim is to find for the first time the limit of predictability of simulated drifter trajectories compared to different observations, in different current regimes.
In the sensitivity experiments we found that Lagrangian trajectories forecast skill largely depends on the accuracy of the input ocean currents: an hourly time frequency and an open-ocean horizontal resolution of only a few km are necessary for recovering drifter trajectories.The present MEDSLIK-II model is then accurate in reproducing drifter trajectories for 1-2.5 days depending on the current conditions.Those results are consistent with the experience of the rapid response to Deepwater Horizon oil spill (Liu et al., 2011b, c).
In the past (Al-Rabeh, 1994;Reed et al., 1994), the drift velocity of the surface oil was considered to be the sum of a fraction of the wind velocity and an estimate of the current fields from OGCM.The wind correction was necessary in order to reproduce the surface Ekman currents, i.e., the local wind effects that were not properly resolved by lowresolution, climatological models.Nowadays, with the advent of accurate operational oceanographic circulation models, a correct representation of the ageostrophic surface current velocity field is provided by the operational OGCM.
Comparing the MEDSLIK-II simulations with drifter trajectories, we therefore prove that there is no need to add a wind correction to reconstruct a correct Ekman current for state-of-the-art operational models such as MFS and AFS, which have upper ocean vertical numerical resolutions of the order of a few metres.Where models have a lower resolution, then corrections allowed by MEDSLIK-II may still be necessary, and each model may develop a calibration matrix for the correction factors.
The use of the wind corrections can still be justified to account for wind drag directly on the drifter, as we argue it is necessary for the IESM-PTR drifter, but for oil slicks it seems unlikely that this correction would be needed unless the quantity of oil is so large that it could modify the airsea interaction physics (Hoult, 1972).In this case, we have yet to obtain a proper representation of the processes, and further investigation is required, especially when there are strong winds.Finally, further investigations are needed to obtain the correct representation of the physical processes in the first mm of the water column, since the thin, interfacial viscous layer could be important in the surface oil spill dynamics and this is not included in any of the present OGCM.In general, wind and wave effects are lumped together and represented by a wind correction coefficient, but the specific role of waves in the slick's drift is important, especially in nearshore areas.Transport by waves (Stokes drift) has been introduced in MEDSLIK-II using an analytical formulation that depends on wind amplitude (using the JONSWAP wave spectrum).We found that adding 1 % of the wind intensity is almost equivalent to considering the Stokes drift velocity.This offers evidence that the wind-correction factor may be used to account for missing wave physics at the air-sea interface.In the future, however, swell and other wave processes should be considered, and MEDSLIK-II should due coupled with a fully-resolved surface wind waves model.
One of the experiments was conducted with an oil slick detected by satellite imagery.We have shown that by changing some uncertain input parameters, such as oil type and slick thickness, the oil concentration simulations are different and the comparison with the satellite imagery can indicate approximately the most likely API value.Moreover, realistic oil concentration distributions are obtained by an optimal oil grid tracer resolution of the order of 100 m and number of particles of the order of a hundred thousand.
Last but not least, the predictability time for oil spill forecasting is of the order of few days maintaining the spatial errors for trajectories within three times the OGCM numerical grid resolution.This implies a frequent re-initialization of the simulation approximately every day along the drifter trajectory positions.The same predictability time window and the need of frequent re-initialization have been also found by Liu et al. (2011b, c).
We believe in the future it will be promising to start an ensemble approach to combine the different model output simulations with uncertain oil spill model parameters, as it has been done during the operational oil spill trajectory modeling effort in the Deepwater Horizon oil spill response (Liu et al., 2011b, c).Among them the most important seem to be the time and space resolution of the advecting current field, the volume of the oil, its thickness and the API value.

Fig. 1 .
Fig. 1.Initialization and forecast of oil spill evolution phases.During initialization the thin and thick areas and thicknesses of the slick state variables are changed.

4
Oil spill simulation and validation experiments 4.1 Sensitivity to the current horizontal resolution, local wind correction and wave correction terms In this first part of the validation study, MEDSLIK-II is used to simulate CODE and IESM-PTR drifters trajectories.CODE drifters were released in the Ligurian Sea (northwestern Mediterranean Sea) in order to understand the importance www.geosci-model-dev.net/6

Fig. 2 .Fig. 3 .
Fig. 2. Observed drifter trajectories (black lines) and the MEDSLIK-II trajectories from 14 May 2007 at 15:00 UTC to 17 May 2007 at 15:00 UTC.(a) The light blue lines are the trajectories obtained using the surface daily MFS currents (CURR-EXP1), the green lines are the trajectories obtained using the surface hourly MFS currents (CURR-EXP2) and the pink lines are the trajectories obtained using the surface hourly currents produced by the IRENOM (CURR-EXP3).(b) The dark blue lines are the trajectories obtained using the 30 m hourly currents produced by MFS and adding a 3 % wind correction with a wind angle of 0 • (CURR-EXP4) and the red lines are the trajectories obtained using the 30 m hourly currents produced by MFS and adding a 3 % wind correction with a wind angle of 25 • (CURR-EXP5).14M. De Dominicis: MEDSLIK-II -Part 2: Numerical simulations a

Fig. 3 .
Fig. 3. RMSE (a) and skill score (b) between the observed and simulated trajectories of Fig. 2 as a function of the prediction time.

Fig. 4 .
Fig. 4. Observed drifter trajectories from 10 to 22 October 2007 (black lines) and the MEDSLIK-II trajectories from 14 to 22 October 2007: (a) drifter 75661, (b) drifter 75662, (c) drifter 75663, (d) drifter 75664, (e) drifter 60212, (f) drifter 60213(Brostrom et al., 2008).Green lines are the trajectories simulated without any correction (WIND-EXP1/SD-EXP1); the red lines are the trajectories simulated using the MFS surface currents and a wind correction of 1 % (WIND-EXP2); the grey lines are the trajectories simulated using the MFS surface currents and a wind correction of 2 % (WIND-EXP3); the light blue lines are the trajectories simulated using the MFS surface currents and a wind correction of 3 % (WIND-EXP4); the blues lines are the trajectories simulated using the MFS currents at 30 m depth and a wind correction of 3 % (WIND-EXP5) and the pink lines are the trajectories simulated using the MFS surface currents and considering Stokes drift velocity (SD-EXP2).Note: in panel (b) the WIND-EXP1 trajectory is not visible because the simulated drifter arrived onto the coast after few hours of simulation.

Fig. 5 .
Fig. 5.As in Fig. 4, but the simulated trajectories last for 24 h and are re-initialized every day, from 10 to 22 October 2007.

Fig. 7 .
Fig. 7.Observed drifter trajectory (black lines) and the MEDSLIK-II trajectories, obtained using the surface hourly AFS currents, the first simulation starts on 21 July 2009 at 09:40 UTC and lasts 15 h, the next simulations start every day at 01:00 UTC and last 24 h: green lines are the trajectories simulated without any correction (SD-EXP3); red lines are the trajectories simulated using the AFS surface current and a wind correction of 1 % (SD-EXP4) and the pink lines are the trajectories using the AFS surface currents and considering Stokes drift velocity (SD-EXP5).

Fig. 10 .Fig. 11 .
Fig. 10.The slick observed by SAR (red) on 6 August 2008 at 09:51 UTC (post-processed data from the related ASAR image, wide swath mode, 400 km, with a 150 m spatial resolution) and the slick observed by the optical sensor (black) on 7 August at 10:50 UTC (post-processed data from the MODIS image).

( a )Fig. 12 .
Fig. 12. Results of experiments of the sensitivity to oil tracer grid resolution and number of particles: MEDSLIK-II-predicted position and concentration corresponding to 7 August 2008 at 10:50 UTC compared with the slick observed by MODIS.

Table 1 .
Table of sensitivity experiments to horizontal current resolution, time frequency and depth of currents.
Table of experiments designed to study the model trajectory's sensitivity to current depth and to local wind correction.
Table of experiments designed to study the model's sensitivity to Stokes drift velocity.

Table 4 .
Oil slick input data provided from satellite image analysis and wind/current fields used.

Table 5 .
Table of the experiments designed to study the model's sensitivity to oil type, slick thickness and slick age.

Table 6 .
Table of experiments designed to study the model's sensitivity to the horizontal resolution of the oil tracer grid and to the number of particles.