MEDSLIK-II , a Lagrangian marine surface oil spill model for short-term forecasting – Part 1 : Theory

The processes of transport, diffusion and transformation of surface oil in seawater can be simulated using a Lagrangian model formalism coupled with Eulerian circulation models. This paper describes the formalism and the conceptual assumptions of a Lagrangian marine surface oil slick numerical model and rewrites the constitutive equations in a modern mathematical framework. The Lagrangian numerical representation of the oil slick requires three different state variables: the slick, the particle and the structural state variables. Transformation processes (evaporation, spreading, dispersion and coastal adhesion) act on the slick state variables, while particle variables are used to model the transport and diffusion processes. The slick and particle variables are recombined together to compute the oil concentration in water, a structural state variable. The mathematical and numerical formulation of oil transport, diffusion and transformation processes described in this paper, together with the many simplifying hypothesis and parameterizations, form the basis of a new, open source Lagrangian surface oil spill model, the so-called MEDSLIK-II, based on its precursor MEDSLIK (Lardner et al. , 1998, 2006; Zodiatis et al. , 2008a). Part 2 of this paper describes the applications of the model to oil spill simulations that allow the validation of the model results and the study of the sensitivity of the simulated oil slick to different model numerical parameterizations.


Introduction
Representing the transport and fate of an oil slick at the sea surface is a formidable task.Many factors affect the motion and transformation of the slick.The most relevant of these are the meteorological and marine conditions at the air-sea interface (wind, waves and water temperature); the chemical characteristics of the oil; its initial volume and release rates; and, finally, the marine currents at different space scales and timescales.All these factors are interrelated and must be considered together to arrive at an accurate numerical representation of oil evolution and movement in seawater.
Oil spill numerical modelling started in the early eighties and, according to state-of-the-art reviews (ASCE, 1996;Reed et al., 1999), a large number of numerical Lagrangian surface oil spill models now exist that are capable of simulating three-dimensional oil transport and fate processes at the surface.However, the analytical and discrete formalism to represent all processes of transport, diffusion and transformation for a Lagrangian surface oil spill model are not adequately described in the literature.An overall framework for the Lagrangian numerical representation of oil slicks at sea is lacking and this paper tries to fill this gap.
Over the years, Lagrangian numerical models have developed complex representations of the relevant processes: starting from two-dimensional point source particle-tracking models such as TESEO-PICHI (Castanedo et al., 2006;Sotillo et al., 2008), we arrive at complex oil slick polygon representations and three-dimensional advection-diffusion models (Wang et al., 2008;Wang and Shen, 2010).At the time being, state-of-the-art published Lagrangian oil spill models do not include the possibility to model threedimensional physical-chemical transformation processes.
The novelty of this paper with respect to the state-of-theart works is the comprehensive explanation on (1) how to reconstruct an oil concentration field from the oil particles advection-diffusion and transformation processes, which has never been described in present-day literature for oil spill models; (2) the description of the different oil spill state variables, i.e. oil slick, oil particles and structural variables; and (3) all the possible corrections to be applied to the ocean current field, when using recently available data sets from numerical oceanographic models.
Our work writes for the first time the conceptual framework for Lagrangian oil spill modelling starting from the Eulerian advection-diffusion and transformation equations.Particular attention is given to the numerical grid where oil concentration is reconstructed, the so-called tracer grid, and in Part 2 sensitivity of the oil concentration field to this grid resolution is clarified.To obtain oil concentrations, here called structural state variables, we need to define particle state variables for the Lagrangian representation of advection-diffusion processes and oil slick variables for the transformation processes.In other words, our Lagrangian formalism does not consider transformation applied to single particles but to bulk oil slick volume state variables.This formalism has been used in an established Lagrangian oil spill model, MEDSLIK (Lardner et al., 1998(Lardner et al., , 2006;;Zodiatis et al., 2008a), but it has never been described in a mathematical and numerical complete form.This has hampered the possibility to study the sensitivity of the numerical simulations to different numerical schemes and parameter assumptions.A new numerical code, based upon the formalism explained in this paper, has been then developed, the so-called MEDSLIK-II, for the first time made available to the research and operational community as an open source code at http://gnoo.bo.ingv.it/MEDSLIKII/(for the technical specifications, see Appendix D).In Part 2 of this paper MEDSLIK-II is validated by comparing the model results with observations and the importance of some of the model assumptions is tested.
MEDSLIK-II includes an innovative treatment of the surface velocity currents used in the Lagrangian advectiondiffusion equations.In this paper, we discuss and formally develop the surface current components to be used from modern state-of-the-art Eulerian operational oceanographic models, now available (Coppini et al., 2011;Zodiatis et al., 2012), considering high-frequency operational model currents, wave-induced Stokes drift and corrections due to winds, to account for uncertainties in the Ekman currents at the surface.
The paper is structured as follows: Sect. 2 gives an overview of the theoretical approach used to connect the transport and fate equations for the oil concentration to a Lagrangian numerical framework; Sect. 3 describes the numerical model solution methods; Sects.4 and 5 present the equations describing the weathering processes; Sect.6 illustrates the Lagrangian equations describing the oil transport processes; Sect.7 discusses the numerical schemes; and Sect.8 offers the conclusions.

Model equations and state variables
The movement of oil in the marine environment is usually attributed to advection by the large-scale flow field, with dispersion caused by turbulent flow components.While the oil moves, its concentration changes due to several physical and chemical processes known as weathering processes.The general equation for a tracer concentration, C(x, y, z, t), with units of mass over volume, mixed in the marine environment, is where ∂ ∂t is the local time-rate-of-change operator, U is the sea current mean field with components (U, V , W ); K is the diffusivity tensor which parameterizes the turbulent effects, and r j (C) are the M transformation rates that modify the tracer concentration by means of physical and chemical transformation processes.
Solving Eq. ( 1) numerically in an Eulerian framework is a well-known problem in oceanographic (Noye, 1987), meteorological and atmospheric chemistry (Gurney et al., 2002(Gurney et al., , 2004) ) and in ecosystem modelling (Sibert et al., 1999).A number of well-documented approximations and implementations have been used over the past 30 yr for both passive and active tracers (Haidvogel and Beckmann, 1999).Other methods use a Lagrangian particle numerical formalism for pollution transport in the atmosphere (Lorimer, 1986;Schreurs et al., 1987;Stohl, 1998).While the Lagrangian modelling approach has been described for atmospheric chemistry models, nothing systematic has been done to justify the Lagrangian formalism for the specific oil slick transport, diffusion and transformation problem and to clarify the connection between the Lagrangian particle approach and the oil concentration reconstruction.
The oil concentration evolution within a Lagrangian formalism is based on some fundamental assumptions.One of the most important of these is the consideration that the constituent particles do not influence water hydrodynamics and processes.This assumption has limitations at the surface of the ocean because floating oil locally modifies air-sea interactions and surface wind drag.Furthermore, the constituent particles move through infinitesimal displacements without inertia (like water parcels) and without interacting amongst themselves.After such infinitesimal displacements, the volume associated with each particle is modified due to the physical and chemical processes acting on the entire slick rather than on the single particles properties.This is a fundamental assumption that differentiates oil slick Lagrangian models from marine biochemical tracer Lagrangian models, where single particles undergo biochemical transformations (Woods, 2002).
If we apply these assumptions to Eq. ( 1), we effectively split the active tracer equation into two component equations: and where C 1 is the oil concentration solution solely due to the weathering processes, while the final time rate of change of C is given by the advection-diffusion acting on C 1 .The model solves Eq. ( 2) by considering the transformation processes acting on the total oil slick volume, and oil slick state variables are defined.The Lagrangian particle formalism is then applied to solve Eq. (3), discretizing the oil slick in particles with associated particle state variables, some of them deduced from the oil slick state variables.The oil concentration is then computed by assembling the particles together with their associated properties.While solving Eq. ( 3) with Lagrangian particles is well known (Griffa, 1996), the connection between Eqs.
(2) and (3), explained in this paper, is completely new.MEDSLIK-II subdivides the concentration C as being composed by the oil concentration at the surface, C S , in the subsurface, C D , adsorbed on the coasts, C C , and sedimented at the bottom, C B (see Fig. 1a).These oil concentration fields are called structural state variables, and they are listed in Table 1.
At the surface, the oil slick is assumed to be represented by a continuous layer of material, and its surface concentration, C S , is defined as with units of kg m −2 , where m is the oil weight and A is the unit area.Considering now volume and density, we write In the subsurface, oil is formed by droplets of various sizes that can coalesce again with the surface oil slick or sediment at the bottom.The subsurface or dispersed oil concentration, C D , can then be written for all droplets composing the dispersed oil volume V D as The weathering processes in Eq. ( 2) are now applied to C S and C D and in particular to oil volumes: dV S dt and ( 7) The surface and dispersed oil volumes, V S and V D , are the basic oil slick state variables of our problem (see Table 1).Equations ( 7) and ( 8) are the MEDSLIK-II equations for the concentration C 1 in Eq. ( 2), being split simply into V S and V D that are changed by weathering processes calculated using the Mackay et al. (1980) fate algorithms that will be reviewed in Sect. 4.
When the surface oil arrives close to the coasts, defined by a reference segment L C , it can be adsorbed and the concentration of oil at the coasts, C C , is defined as where V C is the adsorbed oil volume.The latter is calculated from the oil particle state variables, to be described below, and there is no prognostic equation explicitly written for V C .
The oil sedimented at the bottom is considered to be simply a sink of oil dispersed in the water column, and again it is computed from the oil particles dispersed in the subsurface.
In the present version of the model, the oil concentration on the bottom, C B , is not computed, and it is simply represented by a number of oil particles that reach the bottom.
In order to solve Eqs. ( 7) and ( 8) we need now to subdivide the surface volume into a thin part, V TN , and a thick part, V TK .This is an assumption done in order to use the Mackay et al. (1980) transformation process algorithms.Despite their simplification, Mackay's algorithms have been widely tested, and they were shown to be flexible and robust in operational applications.The surface oil volume is then written as where and where A TK and A TN are the areas occupied by the thick and thin surface slick volume and T TK and T TN are the thicknesses of the thick and thin surface slicks.V TN , V TK , A TN , A TK , T TN and T TK are then oil slick state variables (Table 1) and are used to solve for concentration changes due to weathering processes as explained in Sect. 4. In order to solve the advection-diffusion processes in Eq. ( 3) and compute C S , C D and C C , we define now the particle state variables.The surface volume V S is broken into N constituent particles that are characterized by a particle volume, υ(n k , t), by a particle status index, σ (n k , t), and by a particle position vector: where n k is the particle identification number.The particle position vector x k (n k , t) time evolution is given by the Langevin equation described in Sect.6.Following Mackay's conceptual model, the particle volume state variables are ulteriorly subdivided into the "evaporative" υ E (n k , t) and "non-evaporative" υ NE (n k , t) particle volume attributes: The particle volumes υ(n k , t) are updated using empirical formulas that relate them to the time rate of change of oil slick volume state variables, see Sect. 5.
The particle status index, σ (n k , t), identifies the four particle classes correspondent to the four structural state variables: for particles at the surface, σ (n k , t) = 0; for subsurface or dispersed particles, σ (n k , t) = 1; for sedimented particles, σ (n k , t) = 2; and for particles on the coasts, σ (n k , t) = −L i , where L i is a coastline segment index, to be specified later.
To solve the complete advection-diffusion and transformation problem of Eq. ( 1), we need to specify a numerical grid where we can count particles and compute the concentration.There is no analytical relationship between the oil slick and the particle state variables, and we will then proceed to define the spatial numerical grid and the solution methodology.

MEDSLIK-II tracer grid and solution methodology
In order to connect now Eqs. ( 2) and (3), we need to define a discrete oil tracer grid system, x T = (x T , y T ), with a uniform but different grid spacing in the zonal and meridional directions, (δx T , δy T ) (see Fig. 1b).The unit area A defined in Eqs. ( 5) and ( 6) is then A = δx T δy T , and the spatially discretized time evolution equations for the structural and oil slick state variables are dC S dt (x T , y T , t) = ρ δx T δy T dV S dt (x T , y T , t) and ( 15) The coastline is represented by a polygonal chain identified by a sequence of points connecting segments of length δL i , identified by the coastline segment index, L i (see Fig. 1c).The coast is digitised to a resolution appropriate for each segment, which varies from a few metres to a hundred metres for an almost straight coastal segment.The discrete form of Eq. ( 9) is then When the particle state variables are referenced to the oil tracer grid, we can write the relationship between structural and particle state variables, i.e. we can solve for evolution of the oil concentration at the surface, in the subsurface, and at the coasts.The countable ensembles, I S , I D , of surface and subsurface particles contained in an oil tracer grid cell are defined as The discrete surface, C S , and dispersed, C D , oil concentrations are then reconstructed as The oil concentration for particles on the coasts, C C (L i , t), is calculated using I C (L i , t), which is the set of particles "beached" on the coastal segment L i : (20) The concentration of oil on each coastal segment is calculated by In order to solve coherently for the different concentrations using the oil slick and particle state variable equations, a sequential solution method is developed, which is represented schematically in Fig. 2. First, MEDSLIK-II sets the initial conditions for particle variables and slick variables at the surface (see Sect. 3.1).Then, the transformation processes (evaporation, dispersion, spreading) are solved as described in Sect. 4 and in Appendices B1, B2 and B4.The weathering processes are empirical relationships between the oil slick volume, the 10 m wind, W , and the sea surface temperature, T .Next, the particle volumes, υ NE (n k , t) and υ E (n k , t), are updated (see Sect. 5).Then, the change of particle positions is calculated as described in Sect.6, together with the update of the particle status index.Finally, MEDSLIK-II calculates the oil concentration as described by Eqs. ( 19) and ( 21).
The most significant approximation in MEDSLIK-II is that the oil slick state variables depend only on the slick's central geographical position, which is updated after each advection-diffusion time step.The oil spill centre position, x C = (x C (t), y C (t)), defined by is then used for all the slick state variables of MEDSLIK-II (see Table 1).To evaluate the error connected with this assumption, we estimated the spatial variability of sea surface temperature and compare with a typical linear length scale of an operational oil slick, considered to be of the order of 10-50 km.In the Mediterranean, the root mean square of sea surface temperature is about 0.2 • C for distances of 10 km and 0.5 • C for distances of 50 km.Naturally, across large ocean frontal systems, like the Gulf Stream or the Kuroshio, these differences can be larger, of the order to several • C in 10 km.The calculation of the oil weathering processes, considering the wind and sea surface temperature non-uniformity for the oil slick state variables, will be part of a future improvement of the model.

Initial conditions
The surface oil release can be instantaneous or continuous.
In the case of an oil spill for which leakage may last for several hours or even months (Liu et al., 2011a), it may happen that the earlier volumes of oil spilled will have been transported away from the initial release site by the time the later volumes are released.In order to model the oil weathering in the case of a continuous release, the model divides the total spill into a number of sub-spills, N S , consisting of a given part of the oil released during a time interval, T C .As each sub-spill is moved away from the source, the total spill becomes a chain of sub-spills.In the case of an instantaneous release, the surface oil release at the beginning of the simulation is equal to the total oil released V S (x C , t 0 ).For a continuous oil spill release, every T C a sub-spill is defined with the following oil volume: where R C is the oil spill rate in m 3 s −1 and T C is the time interval between each sub-spill release.The number of subspills released is equal to where D C (s) is the release duration.
During an instantaneous release, N particles are released at the beginning of the simulation, while for a continuous release N C particles are released every T C : Each initial particle volume, υ(n k , t 0 ), is defined as where in the case of an instantaneous release N S is equal to 1.
The initial evaporative and non-evaporative oil volume components, for both instantaneous and continuous release, are defined as )υ(n k , t 0 ) and ( 27) where ϕ NE is the percentage of the non-evaporative component of the oil that depends on the oil type.The initialization of the thin and thick area values is taken from the initial surface amount of oil released using the relative thicknesses and F , which is the area ratio of the two slick parts, A TK and A TN .Using Eqs. ( 10), ( 11) and ( 12), we therefore write ) and ( 29) The same formula is valid for both instantaneous or continuous release.The initial values T TK (x C , t 0 ), T TN (x C , t 0 ) and F have to be defined as input.F can be in a range between 1 and 1000, standard T TK (x C , t 0 ) are between 1 × 10 −4 − 0.02 m, while T TN (x C , t 0 ) lies between 1 × 10 −6 and 1 × 10 −5 m (standard values are summarized in Table 2).For a pointwise oil spill source higher values of T TK (x C , t 0 ) and T TN (x C , t 0 ) and lower values of F are recommended.For initially extended oil slicks at the surface (i.e.slicks observed by satellite or aircraft), lower thicknesses and higher values of F can be used.In the latter case, the initial slick area, A = A TN + A TK , can be provided by satellite images and the thicknesses extracted from other information.

Time rate of change of slick state variables
Using Eq. ( 10), the time rate of change of oil volume is written as The changes of the surface oil volume are attributable to three main processes, known collectively as weathering, which are represented schematically in Fig. 3. Since the initial volume is at the surface, the first process is evaporation.In general, the lighter fractions of oil will disappear, while the remaining fractions can be dispersed below the water surface.In addition, for the first several hours, a given spill spreads mechanically over the water surface under the action of gravitational forces.In the case of a continuous release, the weathering processes are considered independently for each sub-spill.
The weathering processes are considered separately for the thick slick and thin slick (or sheen) and the prognostic equations are written as and where the suffixes indicate evaporation (E), dispersion (D) and spreading (S), and all the slick state variables are defined only at the slick centre.The slick state variables' time rate of change is given in terms of modified Mackay fate algorithms for evaporation, dispersion and spreading (Mackay et al., 1979(Mackay et al., , 1980)).In Appendices B1, B2 and B4, each term in Eqs. ( 32) and ( 33) is described in detail.The model can also simulate the mixing of the water with the oil, and this process known as emulsification is described in Appendix B3.Following Mackay's assumptions, T TN does not change and T TN (x C , t) = T TN (x C , t 0 ).Thus, A TN is calculated as where V TN is updated using Eq. ( 33).
For the thick slick, on the other hand The area of the thick slick, A TK , only changes due to spreading, thus where the time rate of change of the thick area due to spreading is given by Eq. (B20).V TK is updated using Eq. ( 32) and the thickness changes are calculated diagnostically by 5 Time rate of change of particle oil volume state variables The particle oil volumes, defined by Eq. ( 14), are changed after the transformation processes have acted on the oil slick variables.For all particle status index σ (n k , t), the evaporative oil particle volume changes following the empirical relationship where f (E) is the fraction of oil evaporated defined as and V TK (x C , t) (E) and V TN (x C , t) (E) are the volumes of oil evaporated from the thick and thin slicks, respectively, calculated using Eqs.(B1) and (B5).
For both "surface" and "dispersed" particles (σ (n k , t) = 0 and σ (n k , t) = 1), the non-evaporative oil component, υ NE (n k , t), does not change, while a certain fraction of the non-evaporative oil component of a beached particle can be modified due to adsorption processes occurring on a particular coastal segment, seeping into the sand or forming a tar layer on a rocky shore.For the "beached" particles, the particle non-evaporative oil component is then reduced to where t * 0 is the instant at which the particle passes from surface to beached status and vice versa, T S (L i ) is a half-life for seepage or any other mode of permanent attachment to the coasts.Half-life is a parameter which describes the "absorbency" of the shoreline by describing the rate of entrainment of oil after it has landed at a given shoreline (Shen et al., 1987).The half-life depends on the coastal type, for example sand beach or rocky coastline.Example values are given in Table 2.

Time rate of change of particle positions
The time rate of change of particle positions in the oil tracer grid is given by n k uncoupled Langevin equations: where the tensor A(x k , t) represents what is known as the deterministic part of the flow field, corresponding to the mean field U in Eq. ( 1), while the second term is a stochastic term, representing the diffusion term in Eq. ( 1).The stochastic term is composed of the tensor B(x k , t), which characterizes random motion, and ξ(t), which is a random factor.If we define the Wiener process W (t) = t 0 ξ(s)ds and apply the Ito assumption (Tompson and Gelhar, 1990), Eq. ( 41) becomes equivalent to the Ito stochastic differential equation: where dt is the Lagrangian time step and dW (t) is a random increment.The Wiener process describes the path of a particle due to Brownian motion modelled by independent random increments dW (t) sampled from a normal distribution with zero mean, dW (t) = 0 and second order moment with dW • dW = dt.Thus, we can replace dW (t) in Eq. ( 42) with a vector Z of independent random numbers, normally distributed, i.e.Z ∈ N (0, 1), and multiplied by The unknown tensors A(x k , t) and B(x k , t) in Eq. ( 43) are most commonly written as (Risken, 1989): where A was assumed to be diagonal and equal to the Eulerian field velocity components, B is again diagonal and equal to K x , K y , K z turbulent diffusivity coefficients in the three directions, and Z 1 , Z 2 , Z 3 are random vector amplitudes.For particles at the surface and dispersed, Eq. ( 45) takes the following form: where for simplicity we have indicated with dx k (t), dy k (t), dz k (t) the turbulent transport terms written in Eq. ( 45).For particles at the surface, the vertical position does not change: z k = 0 and dz k (t) = 0.The z k can only change when the particles become dispersed and the horizontal velocity at the vertical position of the particle is used to displace the dispersed particles.
The deterministic transport terms in Eq. ( 45) are now expanded in different components: where U C , is the Eulerian current velocity term due to a combination of non-local wind and buoyancy forcings, mainly coming from operational oceanographic numerical model forecasts or analyses; U W , called hereafter the local wind velocity term, is a velocity correction term due mainly to errors in simulating the wind-driven mean surface currents (Ekman currents); and U S , called hereafter the wave current term, is the velocity due to wave-induced currents or Stokes drift.In the following two subsections we will describe the different velocity components introduced in Eq. ( 46).

Current and local wind velocity terms
Ocean currents near the ocean surface are attributable to the effects of atmospheric forcing, which can be subdivided into two main categories, buoyancy fluxes and wind stresses.Wind stress forcing is by far the more important in terms of kinetic energy of the induced motion, accounting for 70 % or more of current amplitude over the oceans (Wunsch, 1998).One part of wind-induced currents is attributable to non-local winds, and is dominated by geostrophic or quasi-geostrophic dynamic balances (Pedlosky, 1986).By definition, geostrophic and quasi-geostrophic motion has a timescale of several days and characterizes oceanic mesoscale motion, a very important component of the largescale flow field included in U .It is customary to indicate that geostrophic or quasi-geostrophic currents dominate below the mixed layer, even though they can sometimes emerge and be dominant in the upper layer.The mixed layer dynamics are typically considered to be ageostrophic, and the dominant time-dependent, wind-induced currents in the surface layer are the Ekman currents due to local winds (Price et al., 1987;Lenn and Chereskin, 2009).All these components should be adequately considered in the U C field of Eq. ( 46).In the past, oil spill modellers computed U C (x k , t) from climatological data using the geostrophic assumption (Al-Rabeh et al., 2000).The ageostrophic Ekman current components were thus added by the term U W (x k , t).It is well known that Ekman currents at the surface U W = (U W , V W ) can be parameterized as a function of wind intensity and angle between winds and currents, i.e.
where W x and W y are the wind zonal and meridional components at 10 m, respectively, and α and β are two parameters referred to as drift factor and drift angle.There has been considerable dispute among modellers on the choice of the best values of the drift factor and angle, with most models using a value of around 3 % for the former and between 0 • and 25 • for the latter (Al-Rabeh et al., 2000).
With the advent of operational oceanography and accurate operational models of circulation (Pinardi and Coppini, 2010;Pinardi et al., 2003;Zodiatis et al., 2008b), current velocity fields can be provided by analyses and forecasts, available hourly or daily, produced by high-resolution ocean general circulation models (OGCMs).The wind drift term as reported in Eq. ( 47) may be optional when using surface currents coming from an oceanographic model that resolves the upper ocean layer dynamics, as also found by Liu et al. (2011b) and Huntley et al. (2011).In such cases, adding U W (x k , t) could worsen the results, as shown in Fig. 2 of Part 2. When the wind drift term is used with a 0 • deviation angle, this term should not be considered as an Ekman current correction, but a term that could account for other nearsurface processes that drive the movement of the oil slick, as shown in one case study of Part 2 (Fig. 4).This theme will be revisited in Part 2 of this paper, where the sensitivity of Lagrangian trajectories to the different corrections applied to the ocean current field will be assessed.

Wave current term
Waves give rise to transport of pollutants by wave-induced velocities that are known as Stokes drift velocity, U S (x k , t) (see Appendix C).This current component should certainly be added to the current velocity field from OGCMs (Sobey and Barker, 1997;Pugliese Carratelli et al., 2011;Röhrs et al., 2012), as normally most ocean models are not coupled with wave models.Stokes drift is the net displacement of a particle in a fluid due to wave motion, resulting essentially from the fact that the particle moves faster forward when the particle is at the top of the wave circular orbit than it does backward when it is at the bottom of its orbit.Stokes drift has been introduced into MEDSLIK-II using an analytical formulation that depends on wind amplitude.In the future, Stokes drift should come from complex wave models, run in parallel with MEDSLIK-II.
Considering the surface, the Stokes drift velocity intensity in the direction of the wave propagation is (see Appendix C) where ω is angular frequency, k is wave-number, and S(ω) is wave spectrum.Equation ( 48) has been implemented in MEDSLIK-II by considering the direction of wave propagation to be equal to the wind direction.The Stokes drift velocity components, U S , are www.geosci-model-dev.net/6/1851/2013/where ϑ = arctg W x W y is the wind direction, and W x and W y are the 10 m height wind zonal and meridional components.

Turbulent diffusivity terms
It is preferable to parameterize the normally distributed random vector Z in Eq. ( 42) with a random number generator uniformly distributed between 0 and 1.We assume that the particle moving through the fluid receives a random impulse at each time step, due to the action of incoherent turbulent motions, and that it has no memory of its previous turbulent displacement.This can be written as where d is the particle mean path and r is a random real number taking values between 0 and 1 with a uniform distribution.The mean square displacement of Eq. ( 50) is while the mean square displacement of the turbulent terms in Eq. ( 45) is simply dx k (t) 2 = 2Kdt.Equating the mean square displacements, we have Finally, the stochastic transport terms in MEDSLIK-II are then written as where K h and K v are prescribed turbulent horizontal and vertical diffusivities.As for modern high resolution Eulerian models, horizontal diffusivity is considered to be isotropic and the values used are in the range 1-100 m 2 s −1 , consistent with the estimation of Lagrangian diffusivity carried out by De Dominicis et al. ( 2012) and indicated by ASCE (1996).
Regarding the vertical diffusion, the vertical diffusivity in the mixed layer, assumed to be 30 m deep, is set to 0.01 m 2 s −1 , while below it is 0.0001 m 2 s −1 (see Table 2).This values is intermediate between the molecular viscosity value for water, i.e. 10 −6 m 2 s −1 , usually reached below 1000 m, and the mixed layer values.

Numerical considerations
Numerical considerations for MEDSLIK-II are connected to the interpolation method between input fields and the oil tracer grid, to the numerical scheme used to solve Eqs. ( 32), ( 33) and ( 45), to the model time step and to the oil tracer grid selection.

Interpolation method
The environmental variables of interest are the atmospheric wind, the ocean currents and the sea surface temperature.They are normally supplied on a different numerical grid than the oil slick centre or particle locations.For the advection calculation, interpolation is thus required to compute the currents and winds at the particle locations.While for the transformation processes calculation, sea surface temperature and winds are interpolated at the slick centre.
Let us indicate with (x E , y E , z E ) the numerical grid on which the environmental variables, collectively indicated by q, are provided by the Eulerian meteorological/oceanographic models.
First, a preprocessing procedure is needed to reconstruct the currents in the zone between the last water grid node of the oceanographic model and the real coastline.MEDSLIK-II employs a procedure to "extrapolate" the currents over land points and thus to add a velocity field value on land.If (x E (i),y E (i)) is considered to be a land grid node by the model, the current velocities component, q x E (i), y E (i) , at the coastal grid point (x E (i), y E (i)), is set equal to the average of the nearby values, when there are at least two neighbouring points (N WP >= 2); that means The result of this extrapolation is shown in Fig. 4. If the current velocities components are given on a staggered grid, a further initial interpolation is also needed to bring both components on the same grid point before the extrapolation is done.
Then, the winds and currents are computed at the particle position (x k , y k ), for a fixed depth z E , with the following interpolation algorithm: where (x k , y k ) is the particle position referenced to the oil tracer grid, (x E (i), y E (i)), (x E (i + 1), y E (i)), (x E (i + 1), y E (i+1)), and (x E (i), y E (i+1)) are the four external field grid points nearest the particle position and x E , y E are the horizontal grid spacings of the Eulerian model (oceanographic or meteorological).Using the same algorithm, the wind and sea surface temperature are interpolated to the oil slick centre, (x C (t), y C (t)), defined by Eq. ( 22).
A vertical interpolation of the currents at the particle positions is also needed and it is computed as follows: where z E (i) and z E (i + 1) are the two Eulerian model levels nearest the particle depth.

Numerical time integration scheme
The Lagrangian horizontal particle motion Eq. ( 40) are solved using a Euler forward scheme.The particle position at time step t + t is calculated as follows: where x k (t) represents the particle position at the current time step, t is the Lagrangian time step, normally taken to be 1800 s, U (x k , t) is the Eulerian ocean current velocity for the current time step at the particle position, and x k (t) is the particle displacement due to turbulent motion.To obtain the Eulerian velocity field at the current time step, another linear interpolation in time between successive input velocity field is carried out.Equations ( 32) and ( 33) are solved again with a Euler forward time stepping scheme but with a different time step, so-called weathering time step, indicated by δt, i.e.
V TK (t + δt) = V TK (t) + dV TK dt δt and (58) where dV TK dt and dV TN dt are given by Eqs. ( 30) and ( 31).The model contains both fast processes (transformation processes) and slower processes (advection-diffusion processes).This generally creates problems for most numerical methods of solving ordinary differential equations.The transformation equations are stiff and to integrate them, the time step should be a fraction of the Lagrangian time step, as done in other active tracer modelling (Butenschön et al., 2012).That is why in MEDSLIK-II the weathering time step has been imposed to be smaller than the Lagrangian time step, typically δt = t 30 .

Particle status updates
The particle oil volumes and the particle status are updated after the particles have moved for a Lagrangian time step ( t).After this movement, the surface particle can become a dispersed particle if the probability function becomes greater than a random number, r, defined to be between 0 and 1.In other words, Here, f (D) (x C , t) is defined as where V TK (x C , t) (D) and V TN (x C , t) (D) is the volume of oil dispersed beneath the thick and thin slicks, respectively calculated using Eqs.(B8) and (B12).The change of oil particle status due to adhesion onto the coast is done by checking whether the parcel intersects any of the line segments, L i , that are used to approximate the coastline.If the particle crosses the coastline, it is moved to the intersecting position.The particle status thus changes from "on surface" to "beached": The beaching of a particle may not be permanent and it is assumed that at subsequent time steps there is a probability that the parcel may be washed back into the water (Shen et al., 1987;Al-Rabeh et al., 2000).The probability of washback is given by  2. At each time step, for each "beached" particle a random number generator, r, is called up and the parcel is released back into the water (its status returns to "surface") if When a particle is washed back, it is of course depleted by the oil that has become permanently attached (Eq.40) and its new position is calculated using Eq. ( 46).The model does not follow any further the oil fraction that is permanently deposited on the coast.It is important to realize that whole particles are not lost as permanently beached, but only a fraction of them.The actual number of particles remains constant.
The deposition of a particle on the bottom is the only case in which a particle is lost from the model.However, no proper parameterization of sedimentation is now included into the model.The particles are considered to be lost from the water column when they are less than 20 cm from the bottom.Thus, the particle status changes from "dispersed" to "sedimented" can be written as where H B (x k , y k ) is the bottom depth below the particle position.

Oil tracer grid and number of particles
The oil tracer grid resolution, δx T , and the total number of particles, N , used to discretize the oil concentration for advection and diffusion processes are important numerical considerations for ensuring the correct reproduction of oil distribution in space and time.
Regarding the oil tracer grid resolution, the scale analysis of the stochastic Eq. ( 42) gives us two limiting spatial scales: where L A is the advective scale (considering U = 0.1 m s −1 and a model time step t = 1800 s), while L T is the diffusion scale (considering a diffusivity K = 2 m 2 s −1 ).The oil tracer grid spatial resolution, δx T , and the model time step must be chosen in order to have A method is needed for estimating the required number of particles and the minimum oil tracer grid spatial resolution.The oil concentration on the water surface (see Eq. 19) at the initial time can be written in the limit of one particle in the tracer grid cell and assuming no evaporation and beaching, using Eq. ( 26), as Deciding which minimum/maximum concentration is possible for any given problem, we can use Eq. ( 69) to find the maximum/minimum number of particles, given a (δx T , δy T ).Thus where N S is equal to 1 in the case of an instantaneous release.Equation ( 70) can be used to provide an estimate of the number of particles for a given spill scenario and oil tracer grid discretization, knowing the lower concentration level of interest.In Part 2 of this paper (De Dominicis et al., 2013), several sensitivity experiments will be carried out to show the impact of different choices regarding number of particles and tracer grid spatial resolution.

Conclusions
This paper presents a formal description of a Lagrangian marine surface oil spill model with surface weathering processes included.An accurate description of a state-of-the-art oil spill model is lacking in the scientific literature.Hand in hand with the release of MEDSLIK-II as an open source model, we want to make available the accurate description of the theoretical framework behind an oil spill model, so as to facilitate understanding of the many modelling assumptions and enable the model to be further improved in the future.
In particular, this paper focuses on the description of the Lagrangian formalism for the specific oil slick transport, diffusion and transformation problem, with particular attention on the clarification of the connection between the Lagrangian particle approach and the oil concentration reconstruction.In order to solve the advection-diffusion-transformation equation for the oil spill concentration, MEDSLIK-II defines three kinds of model state variables: slick, particle and structural variables.Oil slick state variables are used to solve the transformation processes, that act on the entire surface oil slick, and they give information on the volume, area and thickness of the oil slick.The advection-diffusion processes are solved using a Lagrangian particle formalism, meaning that the oil slick is broken into a number of constituent particles characterized by particle state variables.The model reconstructs the oil concentration by considering three concentration classes: at the surface, dispersed in the water column and on the coast.Those concentration fields are structural state variables that are computed by an appropriate merging of information for oil slick and particle state variables.
The transformation processes considered in MEDSLIK-II are valid for a surface oil release: the oil at the surface can be changed by evaporation and spreading, submerged by dispersion processes or adsorbed on the coast for a certain amount of time.Once the oil is dispersed in the water column, it is affected only by the diffusion and advection processes.In this paper, the oil transformation processes are written in terms of empirical analytical functions that have been generalized from existing finite-difference equations, given originally by Mackay et al. (1980).In the near future we will update the model formulation to consider, first, the improvement of the interpolation and resolution of the environmental conditions for the transformation processes, the space dependent thick : thin ratio, and we will develop a complete threedimensional model maintaining the present formulation at the surface.Moreover, in this paper we have presented in detail the deterministic and stochastic components of the particle trajectory equations and discussed the corrections needed to account for missing or imperfectly resolved transport processes with reference to the operational oceanographic analyses and forecasts available.The model now includes a proper representation of high frequency currents, wave-induced currents and wind field corrections used in the advective components.In Part 2 of this paper, MEDSLIK-II is applied to realistic case studies and the importance of model assumptions and corrections is tested.
At this time, MEDSLIK-II does not include the modelling of three-dimensional physical-chemical transformation processes.A complete three-dimensional oil spill model needs to be developed and we argue that MEDSLIK-II offers a good platform for this.While surface processes could remain practically unchanged, new state variables should be defined for the subsurface transformation processes and again connected to the MEDSLIK-II present formulation state variables.Even if three-dimensional Lagrangian trajectory equations have been used by different oceanographic communities, particular work will be required to adapt modelling assumptions to the specific oil transformation processes in the water column.This paper might offer the necessary detailed description of the present-day Lagrangian oil spill model assumptions so that the extension to threedimensional marine oil dynamics and transformation will be possible in the near future.
As it can be seen, Stokes drift velocity is a nonlinear quantity in terms of wave amplitude, and it decays exponentially with depth.
In MEDSLIK-II, the Stokes drift calculation is based on a discrete wave spectrum approach.We start from the following two expressions: the average of the wave spectrum, S, is equal to the variance of the surface displacement, ζ : and the wave energy is related to the variance of sea surface displacement by where ρ W is water density and g is gravity.Then, from Eqs. (C3) and (C4), we obtain the relation between the wave amplitude and wave spectrum: (C6) The wave spectrum, S, to be introduced into Eq.(C6), can be calculated using empirical parameterizations, that describe the wave spectrum as a function of wind speed.We have chosen to use the Joint North Sea Wave Project (JON-SWAP) spectrum parameterization (Hasselmann et al., 1973), taking the wind and fetch into account: where F is the fetch, which is the distance over which the wind blows with constant velocity, and W is the wind velocity intensity at 10 m over the sea surface.In practice, the fetch is calculated as the minimum distance between the oil slick centre and the coast in the opposite direction of the wind direction.

Fig. 1 .
Fig. 1.Schematic view of the oil tracer grids (the grey spheres represent the oil particles): (a) graphical representation of concentration classes; (b) 3-D view of one cell of the oil tracer grid for weathering processes: σ is the particle status index and H B indicates the bottom depth of the δx T , δy T cell ; (c) 2-D view of the oil tracer grid for weathering processes and coastline polygonal chain (red); and (d) 3-D view of the oil tracer grid for advection-diffusion processes.

Fig. 3 .
Fig. 3. Weathering processes using Mackay's approach.TK indicates the thick slick and TN the thin slick.V TK and V TN are the surface oil volumes of the thick and thin part of the slick and the suffixes indicate evaporation (E), dispersion (D) and spreading (S).

Fig. 4 .
Fig. 4. Results of the near coast extrapolation procedure: in red the original hydrodynamic current field and in black the extrapolated one.
parameters r, ε, ω p , γ , and φ were determined during the JONSWAP experiment and are expressed by the following formulae:

Table 1 .
Oil spill model state variables.Four are structural state variables or concentrations, eight are oil slick state variables used for the transformation processes, four are particle state variables used to solve for the advection-diffusion processes.
T TK (x, y, t) Slick Surface thickness of the thick part of the surface oil slick volume m T TN (x, y, t) Slick Surface thickness of the thin part of the surface oil slick volume m x k (t) = (x k (t), y k (t), z k (t)) Particle

1869, 2013 coastal
where T W (L i ) is the half-life of beached oil before it is washed off again.A value of T W (L segment depending on the coastal type.Example values are given in Table i ) is assigned to each www.geosci-model-dev.net/6/1851/2013/Geosci.Model Dev., 6, 1851-

Table 2 .
Model parameters (following the order of appearance in the text).