Geosci. Model Dev., 6, 1851-1869, 2013 www.geosci-model-dev.net/6/1851/2013/ doi:10.5194/gmd-6-1851-2013 © Author(s) 2013. CC Attribution 3.0 License.

Geoscientific

Model Development e

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

M. De Dominicis1, N. Pinardi2, G. Zodiatis3, and R. Lardner3

1Istituto Nazionale di Geofisica e Vulcanologia, Bologna, Italy 2Corso di Scienze Ambientali, University of Bologna, Ravenna, Italy 3 Oceanography Centre, University of Cyprus, Nicosia, Cyprus

Correspondence to: M. De Dominicis (michela.dedominicis@bo.ingv.it)

Received: 25 January 2013 - Published in Geosci. Model Dev. Discuss.: 8 March 2013 Revised: 12 July 2013 - Accepted: 17 July 2013 - Published: 1 November 2013

Abstract. 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.

1 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 three-dimensional physical-chemical transformation processes.

Some of the most sophisticated Lagrangian operational models are COZOIL (Reed et al., 1989), SINTEF OSCAR 2000 (Reed et al., 1995), OILMAP (Spaulding et al., 1994; ASA, 1997), GULFSPILL (Al-Rabeh et al., 2000), ADIOS

(Lehr et al., 2002), MOTHY (Daniel et al., 2003), MOHID (Carracedo et al., 2006), the POSEIDON OSM (Pollani et al., 2001; Nittis et al., 2006), OD3D (Hackett et al., 2006), the Seatrack Web SMHI model (Ambjern, 2007), MEDSLIK (Lardner et al., 1998, 2006; Zodiatis et al., 2008a), GNOME (Zelenke et al., 2012) and OILTRANS (Berry et al., 2012). In all these papers equations and approximations are seldom given and the results are given as positions of the oil slick particles and time evolution of the total oil volume. Moreover, the Lagrangian equations are written without a connection to the Eulerian advection-diffusion active tracer equations even though in few cases (Wang and Shen, 2010) the results are given in terms of oil concentration.

The novelty of this paper with respect to the state-of-the-art 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, 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 advection-diffusion 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 La-grangian 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.

2 Model equations and state variables

The movement of oil in the marine environment is usually attributed to advectionby 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

— + U ■ VC = V ■ (KVC) + ^rj(x, C(x, t), t), (1) dt j=1

where d 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 rj(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, 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 La-grangian 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

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.

Variable Variable type Variable name Dimensions

Cs(x, y,t) Structural Oil concentration at the surface kg m- 2

Cd (x, y, t ) Structural Oil concentration dispersed kg m- 2

CC(x,y,t) Structural Oil concentration on the coast kg m- -1

CB(x,y,t) Structural Oil concentration at the bottom kg m- 2

Vs(x,y,t) Slick Oil slick surface volume 3 m

VD(x,y,t) Slick Oil slick subsurface (dispersed) volume 3 m

VTK(x, y,t) Slick Thick part of the surface oil slick volume 3 m

V TN (x,y,t) Slick Thin part of the surface oil slick volume 3 m

jTK(t) Slick Surface area of the thick part of the surface oil slick volume 2 m

jTN(t) Slick Surface area of the thin part of the surface oil slick volume 2 m

TTK(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) = (xk(t),yk(t),zk(t)) Particle Particle position m

UNE(nk,t) Particle Non-evaporative surface oil volume particle attribute 3 m

VE(nk,t) Particle Evaporative surface oil volume particle attribute 3 m

a(nk,t) = 0,1,2, < 0 Particle Particle status index (on surface, dispersed, sedimented, on coast) -

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:

dCi -r-^

— = 2_rj(x, Ci (x,o,0

is given by the advection-diffusion acting on C1. 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, CS, in the subsurface, CD, adsorbed on the coasts, CC, and sedimented at the bottom, CB (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, CS, is defined as

— = -U -VCi + V-(KVCi), at

where C1 is the oil concentration solution solely due to the weathering processes, while the final time rate of change of C

Cs(x, y,t) = j,

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

Cs(x,y,t) = j Vs.

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: a is the particle status index and Hg indicates the bottom depth of the Sxj, Syj 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.

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, CD, can then be written for all droplets composing the dispersed oil volume Vd as

CD(x,yJ) = — VD.

The weathering processes in Eq. (2) are now applied to Cs and Cd and in particular to oil volumes:

dCs p dVs ,

= "T~i— and df A df

dCD p dVD

The surface and dispersed oil volumes, Vs and Vd, 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\ inEq. (2), being split simply into Vs and VD 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 Lc, it can be adsorbed and the concentration of oil at the coasts, Cc, is defined as

Cq{x, v, t) = y-Vc, Lc

where Vc 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 Vc.

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, CB, 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, VTN, and a thick part, VTK. 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

VS = VTN + VTK

V TN(x,y,t) = ATN(t)T TN(x,y,t) and

VTK(x, y, t) = ATK(t)TTK(x, y, t).

where ATK and ATN are the areas occupied by the thick and thin surface slick volume and TTK and TTN are the thicknesses of the thick and thin surface slicks. VTN, VTK, ATN, ATK, TTN and TTK 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 CS, CD and CC, we define now the particle state variables. The surface volume VS is broken into N constituent particles that are characterized by a particle volume, t), by a particle status index, a(nk, t), and by a particle position vector:

Xfc(nt,t) = (*k("k,f), yk("k,f),Zk(nk,f)),

k = 1, N,

where is the particle identification number. The particle position vector xk(nk, i) 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" uE(nk,i) and "non-evaporative" uNE(«k,i) particle volume attributes:

u(»k, t) = UE(«k, t) + UNE(«k, t).

The particle volumes u(nk,i) 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, a(nk, i), identifies the four particle classes correspondent to the four structural state variables: forparticles at the surface, a(nk, i) = 0; for subsurface or dispersed particles, a(nk, i) = 1; for sedimented particles, a(nk, i) = 2; and for particles on the coasts, a(nk, i) = -L^ where L; 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.

3 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 r = (xT, yr), with a uniform but different grid spacing in the zonal and meridional directions, (5 xr, 5 yr) (see Fig. 1b). The unit area A defined in Eqs. (5) and (6) is then A = 5xr5 yr, and the spatially dis-cretized time evolution equations for the structural and oil slick state variables are

dCs ( t)

——(xr,yr,t) = dt

(xr,yr,t) =

P dVs 5xt5yT dt

P dVp 5xr5yr dt

(xT, yr, t) and

(Xt, yr,t)

The coastline is represented by a polygonal chain identified by a sequence of points connecting segments of length Li, identified by the coastline segment index, Li (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

Cc(Li,t) =

pVc (Li ,t) 5 L; '

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, 1S, 1D, of surface and subsurface particles contained in an oil tracer grid cell are defined as

Is(xr,yr,t) =

1d(Xt, yr, t) =

xr - ^ < xk(t) < xr + 5-xr «k yr - ^ < yk(t) < yr + a(«k,t) = 0

xr - < xk(t) < xr + «k; yr - ^ < yk(t) < yr + ^ a(«k,t) = 1

The discrete surface, CS, and dispersed, CD, oil concentrations are then reconstructed as

Í Cs(xr ,yr,t) = £»k*/s u(nk,t) [CD(xr, yr,t) = ^»k^/D u(nk,t)'

The oil concentration for particles on the coasts, Cc(L{, t), is calculated using /c(L;,i), which is the set of particles "beached" on the coastal segment L\.

Ic(LiJ) = {nk\a(nkJ) = -Li}.

The concentration of oil on each coastal segment is calculated by

C-c(L{, t) = -y- J2 v(*k,t).

O J-j'i j

njctlc

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 Bl, 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, /) and

i>i:(iii;. i). 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 advcction-difTusion time step. The oil spill centre position, arc = (xc(t), yc {t)), defined by

xc(t)= -; and yc(t) =

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 50km. 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.

3.1 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., 201 la), it may happen

initial conditions and environmental variables

evaporation

dispersion

spreading

update thick and thin slick state variables

update particl es oil volumes

update parti ;le positions

change particle status

calculate concentrations

Fig. 2. MEDSLIK-II model solution procedure methodology.

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, /Vs. consisting of a given part of the oil released during a time interval, 7'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 VsUc> 'o).

For a continuous oil spill release, every 7c a sub-spill is defined with the following oil volume:

Vs(xc,t0) = RcTc,

where Re is the oil spill rate in m3 s_1 and 7r is the time interval between each sub-spill release. The number of sub-spills released is equal to

where Dc (s) is the release duration.

During an instantaneous release, N particles are released at the beginning of the simulation, while for a continuous release Nc particles are released every T< -:

Nc = —. Ns

Each initial particle volume, u(nk, to), is defined as

NsVs(xcJo)

v(nkJ0) =

where in the case of an instantaneous release /Vs is equal to 1.

The initial evaporative and non-evaporative oil volume components, for both instantaneous and continuous release, are defined as

•pne

VE(nk, to) = (1 - — )v(nk, to) and une(nk, to) = to),

where </<\k 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 /'. which is the area ratio of the two slick parts, ATK and A™. Using Eqs. (10), (11) and (12), we therefore write

A™(t0) = FATK(t0) and

ATK(?o) =

Vs(XcJo)

TtHxcJo) + FT™(xcJo)'

The same formula is valid for both instantaneous or continuous release. The initial values 7'IK(xc. 'o). T™(xcJo) and F have to be defined as input. F can be in a range between 1 and 1000, standard TiK(xc.to) are between 1 x 10"4 - 0.02 m, while T™(xc, i0) lies between 1 x 10"6 and 1 x 10"5 m (standard values are summarized in Table 2). For a pointwise oil spill source higher values of 7' I K(xc. to) and 7'ix(a-('. /(>) 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™ + ATK, can be provided by satellite images and the thicknesses extracted from other information.

4 Time rate of change of slick state variables

Using Eq. (10), the time rate of change of oil volume is written as

dVs dt

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.

t dV™ z f dV™

- dt (e) - dt (e)

dVTK dV™

(s) x dt 7 (d) x dt 7

Fig. 3. Weathering processes using Mackay's approach. TK indi-

frT/ TNT

cates the thick slick and TN the thin slick. V and V 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).

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

dVTK dVTK dVTK dVTK

di di (e) di (d) di (s)

dV™ dV™ dV™ dV™

di di (e) di (d) di (s)

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, 1980). In Appendices Bl, B2 andB4, 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 emulsifi-cation is described in Appendix B3.

Following Mackay's assumptions, 7'IX does not change and T™(xc, t) = T™(xc, t0). Thus, A™ is calculated as

T™ di

where V™ is updated using Eq. (33). For the thick slick, on the other hand

The area of the thick slick, spreading, thus

only changes due to 6 Time rate of change of particle positions

where the time rate of change of the thick area due to spreading is given by Eq. (B20). VTK is updated using Eq. (32) and the thickness changes are calculated diagnostically by

T tk =

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 a(nk, t), the evaporative oil particle volume changes following the empirical relationship

A 100/ c -

u(ftk,io),

where f(E) is the fraction of oil evaporated defined as

/(E)(xc,t) =

V TK(x c, t) | (E) + VTN (x c, t)

V TK(to) + V™ (to)

and VTK(xc,t)|(E) and VTN(xc, 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 (a (nk, t) = 0 and a(nk,t) = 1), the non-evaporative oil component, uNE(nk, 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

UNE(nk,t) = UNE(«k,t0*)0.5 rs(ii), a(nk,t) = -i, (40)

where tf is the instant at which the particle passes from surface to beached status and vice versa, TS(Li) is a half-life for seepage or any other mode of permanent attachment to the coasts. Half-life is a parameter which describes the "ab-sorbency" of the shoreline by describing the rate of entrain-ment 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.

The time rate of change of particle positions in the oil tracer grid is given by uncoupled Langevin equations:

dxk(t) dt

= A(xk,t) + B(x k,t)f(t),

where the tensor A(xk, t) represents what is known as the deterministic part of the flow field, corresponding to the mean field U inEq. (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(xt), which characterizes random motion, and f (t), which is a random factor. If we define the Wiener process W(t) = /0 f (s)ds and apply the Ito assumption (Tompson and Gelhar, 1990), Eq. (41) becomes equivalent to the Ito stochastic differential equation:

dxk(t) =A(xk, t)dt + B(x k, t)dW (t),

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 e N(0,1), and multiplied by Vdi:

dxk(t) =A(xk, t)dt + B(x k, t)ZVdi.

The unknown tensors A(xk, t) and B(x k, t) inEq. (43) are most commonly written as (Risken, 1989):

dx k(t) = U(xk,t) ' V(xt,t) W(xk, t)

V2Kx 0 0

o o ' o

o V2KZ

(44) Z2 j Vdi

where A was assumed to be diagonal and equal to the Eule-rian field velocity components, B is again diagonal and equal to , , turbulent diffusivity coefficients in the three directions, and Z1, Z2, Z3 are random vector amplitudes. For particles at the surface and dispersed, Eq. (45) takes the following form:

dxk(t) =

/U(xk,yk, zk,t)\ V(xk,yk,zk,t)

/d4 (t)\

dyk (t)

\dzk(t)/

where for simplicity we have indicated with dx[ (t), dy[ (t), dzk(t) the turbulent transport terms written in Eq. (45). For particles at the surface, the vertical position does not change: zk = 0 and dzk (t) = 0. The zk 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:

' a = 0 dxk(t) = [Uc(xk, yk, 0, t) + Uw(xk, yk, t)

+U s(xk,yk,t)] dt + dx k(t) , (46)

.a = 1 dxk(t) = U c(xk,yk,zk,t)dt + dx k(t)

where Uc, 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 US, 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).

6.1 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 large-scale 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 Uc field of Eq. (46). In the past, oil spill modellers computed Uc(xk,t) from clima-tological data using the geostrophic assumption (Al-Rabeh et al., 2000). The ageostrophic Ekman current components were thus added by the term UW(xk,t). It is well known that Ekman currents at the surface UW = (UW, VW) can be parameterized as a function of wind intensity and angle between winds and currents, i.e.

UW = a (Wx cos p + Wy sin p) and

Vw = a (-Wx sinp + Wy cosp) (47)

where Wx and Wy are the wind zonal and meridional components at 10 m, respectively, and a and p 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 UW(xk,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 near-surface 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.

6.2 Wave current term

Waves give rise to transport of pollutants by wave-induced velocities that are known as Stokes drift velocity, US(xk, 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; Rohrs 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)

Ds(z = 0) = 2 j wk(w)S(w)dw, (48)

where m is angular frequency, k is wave-number, and S(m) 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

Us = Ds cos # and Vs = Ds sin # , (49)

where û = arctg ^ W^) is the wind direction, and Wx and Wy 7.1 Interpolation method

are the 10 m height wind zonal and meridional components. 6.3 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

dxk(i) = (2r - 1)d,

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

(dxk W) = /o [(2r - 1)d]2 dr = Id2

while the mean square displacement of the turbulent terms in Eq. (45) is simply dx^(t)2 = 2Kdt. Equating the mean square displacements, we have

d2 = 6Kx dt d2 = 6Ky dt d2 = 6K7di.

Finally, the stochastic transport terms in MEDSLIK-II are then written as

dxk(i) = ZiV2KXd? = [2r - 1] V6Khd/ dyk(t) = Z^y2KydF = [2r - 1] V6Khd? dzk(í) = Z3V2K¡dí = [2r - 1] V6KVd/,

where and Kv 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 m2 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 m2 s-1, while below it is 0.0001 m2 s-1 (see Table 2). This values is intermediate between the molecular viscosity value for water, i.e. 10-6 m2 s-1, usually reached below 1000 m, and the mixed layer values.

7 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.

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 (xE, yg, zE) the numerical grid on which the environmental variables, collectively indicated by q, are provided by the Eulerian meteorologi-cal/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 (xE(f),yE(f)) is considered to be a land grid node by the model, the current velocities component, qxE(i), yE0'), at the coastal grid point (xE(i), yE(i)), is set equal to the average of the nearby values, when there are at least two neighbouring points (NWP >= 2); that means

qXE(i),yE(!) =

gxE(i+1),yE(0+q*E(i-1),yE(0+q*E(0,yE(i-1)+q*E(0,yE(i+1) NWP '

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 (xk, yk), for a fixed depth zE, with the following interpolation algorithm:

q 1 = q*E(!),yE(i)[*E(i +1) - xk ] q 2 = qxE(!'+1),yE(!')[xk - xr(i)] q 3 = qxE(!'),yE(!'+1)[xE (i +1) - ]

q4 = qxE(i+1),yE(;+1)[xfc - xE(i)]

(q 1 + q2)[ye(i + 1) - yk] + (q3 + g4)[yt - ye(»)] axe Aye

where (xk, yk) is the particle position referenced to the oil tracer grid, (xE(i), yE(i)), (xE(i + 1),yE(i)), (xE(i + 1), yE(i + 1)), and (xE(i), yE(i + 1)) are the four external field grid points nearest the particle position and AxE, AyE are the horizontal grid spacings of the Eulerian model (oceano-graphic or meteorological). Using the same algorithm, the wind and sea surface temperature are interpolated to the oil slick centre, (xc(t), yc(t)), defined by Eq. (22).

A vertical interpolation of the currents at the particle positions is also needed and it is computed as follows:

= zE(î')-1E(i'-l^qXk,yk,ZE(! + l)[zE(i) — ] +

+q*k,»,ZE(!')[zk - ZE(i +1)]} where zE(i) and zE(i + 1) are the two Eulerian model levels nearest the particle depth.

7.2 Numerical time integration scheme

The Lagrangian horizontal particle motion Eq. (40) are solved using a Euler forward scheme. The particle position at time step i + Ai is calculated as follows:

x k(i + At) = x k(i) + U (x i)Ai + Axk(i),

where xk(i) represents the particle position at the current time step, Ai is the Lagrangian time step, normally taken to be 1800 s, U(xk, i) is the Eulerian ocean current velocity for the current time step at the particle position, and Ax^(i) 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 Si, i.e.

V TK(i + Si) = V TK(i) +

V TN(i + Si) = V TN(i) +

Si and

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 (Butenschon 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 Si = —.

7.3 Particle status updates

The particle oil volumes and the particle status are updated after the particles have moved for a Lagrangian time step (Ai). After this movement, the surface particle can become a dispersed particle if the probability function

f(D)(xC,i) - f(D)(xc,i - Ai)

P(D)(i) =

(1 - /(D)(xc,i - Ai)) becomes greater than a random number, r, defined to be be-

tween 0 and 1. In other words,

r < P(D)(i) a(«k,(i)) = 1.

_ V N \ \ \ \

S s N —» \ y l \ \ \ \

V 1 t \ \ * s —_ — N

11 E 30.00'

Fig. 4. Results of the near coast extrapolation procedure: in red the original hydrodynamic current field and in black the extrapolated

Here, /(D)(xC, i) is defined as

V TK(x c,i)|(D) + V TN(x c,i)

/(D)(xc ,i) =

V TK(x c,io) + V TN(x c,io)

VTK(xc, i) |(D) and V™(xC,i)|(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^ 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":

xk(i), yk(i)eLi a(«k,(i)) = -Li.

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 wash-back is given by

P (c)(Li,i) = 1 - 0.5 Tw(Li).

where Tw(L;) is the half-life of beached oil before it is washed off again. A value of Tw(L;) is assigned to each

coastal segment depending on the coastal type. Example values are given in Table 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

P(C)(i) a(«k,(i)) = 0.

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

HB(xk, yk) - Zk < 20 cm a(nk, (i)) = 2,

where HB(xk, yk) is the bottom depth below the particle position.

7.4 Oil tracer grid and number of particles

The oil tracer grid resolution, 5xT, and the total number of particles, N, used to discretize the oil concentration for ad-vection 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:

La = UAi « 180 m and LT = ^KAi « 60 m,

where LA is the advective scale (considering U = 0.1ms-and a model time step At = 1800 s), while LT is the diffusion scale (considering a diffusivity K = 2 m2 s-1). The oil tracer grid spatial resolution, 5xT, and the model time step must be chosen in order to have

Lt < 5xT < La.

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

Cs(Xt, yT,i) =

NsVS(x c,i0 ) P

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 (5xT, 5yT ). Thus

,,max NsVs(Xc,t0) ,..min Ns Vs(x C, i0)

N = ——-p and N = ——-p,

CSmm5xT5yT

CSmax5xT5yT

where NS 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.

8 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.

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

Model symbol Name

Default value

T ™(f0)

T TK(t0)

VMOL C(E)

■(D) ■(d)

T0 C. C C

(M) 1(M)

3TK 0(S) ,(s) ,|s)

Total number of particles released

Tracer grid cell size

Standard initial thin slick thickness

Standard initial thick slick thickness

Standard area factor between thin and thick slicks

Half-life for absorption on coast

Half-life for washing off coast

Horizontal diffusivity

Vertical diffusivity (above thermocline and below thermocline) Rate of change in vapour pressure with evaporated fraction Gas constant Molar volume Evaporation rate

Exponent of wind speed in evaporation rate Wind scale

Evaporated fraction on oil viscosity

Small droplet rising velocity

Thickness of the droplet cloud

Downward diffusion velocity of small droplets

Rate of dispersion of all droplets by waves

Dispersion from the thick slick

Viscosity emulsion scale

Dispersion from the thin slick (sheen)

Interfacial surface tension (oil/water)

Interfacial surface tension (oil/water) scale

Water fraction on mousse viscosity

Rate of increase of water fraction

Reciprocal of maximum water fraction

Thick oil thickness scale

Rate of spreading of thin slick

Rate of spreading of thick slick

Dependence of spreading of thin slick on thickness

Thickness offset

90 000 150m 1 x10-5 m 0.02 m 4

96 h (rocky) 24 h (sandy) 18 h (rocky) 24 h (sandy)

0.01 m s-1-1 x 10-4 m s-1 12.0 s-1

8.2 x 10-5 barm3 mol-1 K 2 x10-4 mol m-3 8 x10-4 m s-1 0.78

1ms 4.0

0.0003 ms-1 0.5 m

0.001 ms-1

10 m2 s-1 2000 24 kg s-2 24 kg s-2 0.65

1.333 1m

1.0s-1

150.0 s-1 0.0015 m

1 x10-5 m

0.8 x!0-5 s-1

1.6 10-6s-1

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 three-dimensional 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 three-dimensional marine oil dynamics and transformation will be possible in the near future.

Appendix A Oil density

The oil density depends on the oil type which is classified using the American Petroleum Institute gravity, or API gravity, which is a measure of how heavy or light a petroleum liquid is compared to water. From API gravity it is possible to calculate oil density. The conversion from API to density first requires conversion to specific gravity:

(API +131.5)

The specific gravity can subsequently be converted to density:

p = SGpw,

where pW is the water density assumed to be equal to 1026 kg m-3.

In MEDSLIK-II the density remains constant over time: temperature expansions and emulsification effects are not considered in the density calculation. Using this hypothesis, MEDSLIK-II concentrations are valid only for short-term forecasting and in the absence of abrupt changes of temperature.

Appendix B

Time rate of change of oil slick state variables B1 Evaporation

Evaporation changes the volume of the thick and thin parts of the slick, and is the major transformation process after the initial oil release at the surface. The volume of oil lost by evaporation is computed using Mackay's algorithm for evaporation (Mackay et al., 1980). Given the assumption that transformation processes are evaluated at the slick centre, the time rate of change of the volume lost by evaporation from the thick slick, VTK, is expressed as

[V TK(i0) + V TN(i0)]

where VTK (i0) and VTN (i0) are the initial thick and thin slick volumes, respectively, and is the time rate of change

of the fraction of oil evaporated. For the thick oil slick, the time rate of change of the fraction of oil evaporated is

di and

Poil =

-c/tki

(1 - /tk)

Vmol '

where Po;i (bar) is the oil vapour pressure, P0 is the initial vapour pressure (which depends on the oil type used), c (s-1) is a constant that measures the rate of decrease of vapour pressure with the fraction already evaporated, ATK (i) is the area of the thick part of the slick, KM (m s-1) is the evaporative exposure to wind, T (K) is the temperature, R (barm3 mol-1 K) is the gas constant and VMOL (molm-3) is the molar volume of the oil. For KM we assume

KM = C,(E) (3.6—)Y, M 1 ( W0 ) ,

where WW is the non-dimensional 10 m wind modulus (W0 is 1ms-1), y is a constant, and C(E) (ms-1) is the evaporation rate. The standard values of c, R, VMOL, y and C(E) are given in Table 2.

For the thin slick oil, the time rate of change of the volume is equal to

[VTK (to) + V TN(io)],

where f-

is the time rate of change of the oil fraction

evaporated from the thin slick.

The evaporative component in the thin slick is assumed to disappear immediately, but the thin slick, through the spreading process, is fed by oil from the thick slick that in general has not yet fully evaporated. Equating the oil content of the thin slick before and after the flow of oil coming from the thick slick, we obtain

(/MAX - fTK) VTN ,

where /MAX is the initial fraction of the evaporative component, which represents the maximum value that the oil fraction evaporated from the thin slick can attain. Evaporation leads to an increase in the viscosity of the oil, which is calculated using

n = n0eK(E)/TK,

where n0 (m2 s-1) is the initial viscosity (which depends on the oil type used) and K(E) is a constant that determines the increase of viscosity with evaporation (see Table 2).

B2 Dispersion

The oil dispersion processes occurring on an oil volume released at the surface were framed in empirical formulas developed by Mackay et al. (1979). Wave action drives oil into the water, forming a cloud of droplets beneath the spill. The droplets are classified as either large droplets that rapidly rise and coalesce again with the surface spill, or small droplets that rise more slowly, and may be immersed long enough to diffuse into the lower water column layers. In the latter case, they are lost from the surface spill and considered to be permanently dispersed. The criterion that distinguishes the small droplets is that their rising velocity under buoyancy forces is comparable to their diffusive velocity, while for large droplets the rising velocity is much larger.

The time rate of change of the thick slick volume due to water column dispersal of small droplets is given by

(D) .(D)

2 (c? - vs) CsA

i TK i dXS + di :

where C(D) (ms-1) is the downward diffusive velocity of the small droplets, and vS (m s-1) is the rising velocity of the small droplets. C(D) and vS are constant parameters listed in Table 2; cs is the fraction of the small droplets; and Xs is the volume of small droplets beneath the thick slick. The amount of small droplets is equal to

Xs = CswmA

where wm (m) is the vertical thickness of the droplet cloud (see Table 2). The large droplets are not regarded as dispersed since they eventually re-coalesce with the slick. The fraction of the small droplets is calculated using the following expres-

nEMo, and r0TK are listed in Table 2. For the thin slick dispersion only small droplets are considered. It is assumed that these droplets are all lost from the surface spill at the following rate:

= C3D)(— + 1 3 1 Wo +

TtnAtnStn and

stn =( 1 + c5d)-

(B12) (B13)

where C(D) is control dispersion from the thin slick (see Table 2), WW is the non-dimensional wind speed at the oil slick centre and Stn is the fraction of small droplets in the dispersed oil beneath the thin slick.

B3 Emulsification

Emulsification refers to the process by which water becomes mixed with the oil in the slick. The main effect of emulsifi-cation is to form a mousse with viscosity nEM given by

nEM = n exp

_1 - C(M)/ w

where n is defined by Eq. (B7), fW is the fraction of water in the oil-water mousse and C(M) is a constant which controls the effect of water fraction on mousse viscosity (see Table 2). Emulsification is assumed to continue until nEM reaches a maximum value nMAX corresponding to a mousse composed of floating tar balls. Mackay's model for the time rate of change of fW is (Mackay et al., 1979)

2c3D) (f + 1) ttkstk

vs + q

where c3d) (s 1) is a constant which controls the rate of dispersion of all droplets by waves (see Table 2), W is the non-dimensional wind speed at the oil slick centre, and STK is the fraction of small droplets in the dispersed oil beneath the thick slick, equal to

nEM A2 nEMO/

10-3T0TK

1 - C3m)/ w

where W is the non-dimensional wind speed calculated at the slick centre, C2M) (s-1) is a constant which controls the rate of water absorption in the mousse and c3m) is a constant which controls the maximum water fraction in the mousse (see Table 2).

According to Eqs. (B14) and (B15) emulsification influences the mousse viscosity, which, in turn, influences dispersion (see Eq. B11).

(B11) B4 Spreading

where c4d) controls the fraction of droplets below a critical size, t is the interfacial surface tension between oil and water (kgs-2) and nEM (m2 s-1) is the emulsified oil viscosity

that will be defined later. The standard values of C

T, To,

Spreading consists of two processes: the first is the area lost due to oil converted from the thick to the thin slick and the second corresponds to Fay's gravity-viscous phase of spreading (Al-Rabeh et al., 2000). The thin and thick slick volume

rates due to spreading are then written as

+ T TKFG and

where FG is defined later and correspond to Fay's gravity spreading. Mackay's model (Mackay et al., 1979, 1980) approximates the thin slick area increment by

= C(S) (atn)1/3 (t0tk)4/3

Ttk + e

where c(s) (s-1) is the constant rate of spreading of the thin

slick, and C3 ) (m) controls the dependence on thickness of the spreading of the thin slick, and e is a constant parameter. C(S), C3S) and e standard values are listed in Table 2. For the thick slick, Fay's spreading is assumed to be written as

(S) ( , TK)1/3 (TTK)4/3

FG = C2S) (A

where c2s) (s-1) is a constant rate of spreading of the thick slick (see Table 2).

The time rate of change of the area of the thick slick due to spreading is

TTK di

Mechanical spreading is considered to occur for an initial period of 48 h after the oil release or until the thickness of the thick part of the slick, TTK, determined by Eq. (35), becomes equal to that of the thin slick, TTN. If this occurs the model terminates all further spreading and from that point the slick is modelled as a thin slick only.

Appendix C Stokes drift

Stokes drift velocity is the difference between the average Lagrangian flow velocity of fluid particles and the average Eulerian flow velocity of the fluid at a fixed position (the average is usually taken over one wave period). Stokes drift velocity is given by (Stokes, 1847)

2 , cosh [2k (Z -

Ds («, z) = a «k-=-

2sinh2(KH)

where m is the angular frequency, k is the wave number, a is the wave amplitude and H is the depth of the ocean. The horizontal component of Stokes drift velocity, DS, for deep-water waves (H ^ to) is approximately

DS («,z) ~ «ka2e2kz.

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, Z:

(z2) = j S(«)d«,

and the wave energy is related to the variance of sea surface displacement by

E = pwg (z 2) = 2 Pwga2,

where pW is water density and g is gravity. Then, from Eqs. (C3) and (C4), we obtain the relation between the wave amplitude and wave spectrum:

a2 = 2 (z2j = 2 j S(«)d«.

Introducing Eq. (C5) into Eq. (C2), the Stokes drift formulation becomes

Ds(z) = 2 j «k(«)S(«)e2k(®)zd«.

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 (JONSWAP) spectrum parameterization (Hasselmann et al., 1973), taking the wind and fetch into account:

CV £g

S(«) = —5 exp

5 / «p > 4 V « v

The parameters r, e, «p, y, and 0 were determined during the JONSWAP experiment and are expressed by the following formulae:

r = exp

(« — «p)2 202 «2

,'W 2\ e = 0.076 -

«P = 22 ( — ) ; Y = 3.3; and 0 =

0.07 0.09

« < «p « > «p

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.

From the wave spectrum, the significant wave height can be also calculated. It is defined as the average height of the highest third of the waves, during a fixed sampling interval, and it can be written as

Hs = 4// S(«)d«. (C9)

Appendix D Technical specifications

The oil spill model code MEDSLIK-II is freely available and can be downloaded together with the User Manual, test case data and output example from the website http://gnoo. bo.ingv.it/MEDSLIKII/. MEDSLIK-II is available under the GNU General Public License (Version 3, 29 June 2007). The code is written in Fortran77, Python and Shell scripting. The model can run on any workstation and laptop. The architecture currently supported is Linux (tested on Ubuntu 10.04 LTS). The software requirements are a Fortran compiler (gfortran is fully compatible) and NetCDF libraries.

Acknowledgements. This work was funded by MyOcean Project and MEDESS4MS Project.

Edited by: D. Ham

References

Al-Rabeh, A. H., Lardner, R. W., and Gunay, N.: Gulfspill Version 2.0: a software package for oil spills in the Arabian Gulf, Environ. Modell. Softw., 15, 425-442, 2000.

Ambjorn, C.: Seatrack Web, Forecasts of Oil Spills, a New Version, Environ. Res. Eng. Manage., 3, 60-66, 2007.

ASA: OILMAP for Windows (technical manual), Narrangansett, Rhode Island: ASA Inc, 1997.

ASCE: State-of-the-Art Review of Modeling Transport and Fate of Oil Spills, J. Hydraulic Eng., 122, 594-609, 1996.

Berry, A., Dabrowski, T., and Lyons, K.: The oil spill model OIL-TRANS and its application to the Celtic Sea, Mar. Pollut. Bull., 64,2489-2501,2012.

Butenschon, M. and Zavatarelli, M., and Vichi, M.: Sensitivity of a marine coupled physical biogeochemical model to time resolution, integration scheme and time splitting method, Ocean Model., 52, 36-53,2012.

Carracedo, P., Torres-Lopez, S., Barreiro, M., Montero, P., Balseiro, C., Penabad, E., Leitao, P., and Perez-Munuzuri, V.: Improvement of pollutant drift forecast system applied to the Prestige oil spills in Galicia Coast (NW of Spain): Development of an operational system, Mar. Pollut. Bull., 53, 350-360,2006.

Castanedo, S., Medina, R., Losada, I. J., Vidal, C., Mendez, F. J., Osorio, A., Juanes, J. A., and Puente, A.: The Prestige oil spill in Cantabria Bay of Biscay). Part I: Operational forecasting system for quick response, risk assessment, and protection of natural resources, J. Coast. Res., 22, 1474-1489,2006.

Coppini, G., De Dominicis, M., Zodiatis, G., Lardner, R., Pinardi, N., Santoleri, R., Colella, S., Bignami, F., Hayes, D. R., Soloviev, D., Georgiou, G., and Kallos, G.: Hindcast of oil-spill pollution during the Lebanon crisis in the Eastern Mediterranean, July-August 2006, Mar. Pollut. Bullet., 62, 140-153, 2011.

Daniel, P., Marty, F., Josse, P., Skandrani, C., and Benshila, R.: Improvement of drift calculation in Mothy operational oil spill prediction system, in: International Oil Spill Conference (Vancouver, Canadian Coast Guard and Environment Canada), vol. 6, 2003.

De Dominicis, M., Leuzzi, G., Monti, P., Pinardi, N., and Poulain, P.: Eddy diffusivity derived from drifter data for dispersion model applications, Ocean Dynam., 62, 1381-1398, doi:10.1007/s10236-012-0564-2, 2012.

De Dominicis, M., Pinardi, N., Zodiatis, G., and Archetti, R.: MEDSLIK-II, a Lagrangian marine surface oil spill model for short-term forecasting - Part 2: Numerical simulations and validations, Geosci. Model Dev., 6, 1871-1888, doi:10.5194/gmd-6-1871-2013,2013.

Griffa, A.: Applications of stochastic particle models to oceano-graphic problems, in: Stochastic modelling in physical oceanography, Progress in Probability, 39, 113-140, 1996.

Gurney, K., Law, R., Denning, A., Rayner, P., Baker, D., Bousquet, P., Bruhwiler, L., Chen, Y. H., Ciais, P., Fan, S., Fung, I. Y., Gloor, M., Heimann, M., Higuchi, K., John, J., Maki, T., Maksyutov, S., Masarie, K., Peylin, P., Prather, M., Pak, B.

C., Randerson, J., Sarmiento, J., Taguchi, S., Takahashi, T., and Yuen, C.-W.: Towards robust regional estimates of CO2 sources and sinks using atmospheric transport models, Nature, 415, 626630, 2002.

Gurney, K., Law, R., Denning, A., Rayner, P., Pak, B., Baker, D., Bousquet, P., Bruhwiler, L., Chen, Y., Ciais, P., Fung, I. Y., Heimann, M., John, J., Maki, T., Maksyutov, S., Peylin, P., Prather, M., and Taguchi, S.: Transcom 3 inversion intercom-parison: Model mean results for the estimation of seasonal carbon sources and sinks, Global Biogeochem. Cy., 18, GB1010, doi:10.1029/2003GB002111, 2004.

Hackett, B., Breivik, 0., and Wettre, C.: Forecasting the Drift of Objects and Substances in the Ocean, Ocean Weather Forecasting, Springer, Netherlands, 507-523, 2006.

Haidvogel, D. B. and Beckmann, A.: Numerical ocean circulation modeling, Imperial College Pr, 318 pp., ISBN 9781860941146, 1999.

Hasselmann, K., Barnett, T., Bouws, E., Carlson, H., Cartwright,

D., Enke, K., Ewing, J., Gienapp, H., Hasselmann, D., Kruse-man, P., Meerburg, A., Mller, P., Olbers, D., Richter, K., Sell, W., and Walden, H.: Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP), Ergänzungsheft zur Deutschen Hydrographischen Zeitschrift Reihe, A8-12, 1973.

Huntley, H. S., Lipphardt, B. L., and Kirwan, A. D.: Surface drift predictions ofthe Deepwater Horizon spill: The Lagrangian perspective, in: Monitoring and Modeling the Deepwater Horizon

Oil Spill: A Record-Breaking Enterprise, Geophys. Monogr. Ser., 195, 179-195,2011

Lardner, R., Zodiatis, G., Loizides, L., and Demetropoulos, A.: An operational oil spill model for the Levantine Basin (Eastern Mediterranean Sea), in: International Symposium on Marine Pollution, 1998.

Lardner, R., Zodiatis, G., Hayes, D., and Pinardi, N.: Application of the MEDSLIK Oil Spill Model to the Lebanese Spill of July 2006, European Group of Experts on Satellite Monitoring of Sea Based Oil Pollution, European Communities, 2006.

Lehr, W., Jones, R., Evans, M., Simecek-Beatty, D., and Overstreet, R.: Revisions of the ADIOS oil spill model, Environ. Modell. Softw., 17, 189-197,2002.

Lenn, Y. D. and Chereskin, T. K.: Observations of Ekman currents in the Southern Ocean, J. Phys. Oceanogr., 39, 768-779, 2009.

Liu, Y., MacFadyen, A., Ji, Z.-G., and Weisberg, R. H.: Introduction to Monitoring and Modeling the Deepwater Horizon Oil Spill, in: Monitoring and Modeling the Deepwater Horizon Oil Spill: A Record-Breaking Enterprise, Geophys. Monogr. Ser., 195, 1-7, 2011a.

Liu, Y., Weisberg, R. H., Hu, C., and Zheng, L.: Trajectory forecast as a rapid response to the Deepwater Horizon oil spill, in: Monitoring and Modeling the Deepwater Horizon Oil Spill: A Record-Breaking Enterprise, Geophys. Monogr. Ser., 195, 153-165,2011b.

Lorimer, G.: The kernel method for air quality modelling-1, Mathematical foundation, Atmos. Environ. (1967), 20, 1447-1452, 1986.

Mackay, D., Buist, I., Mascarenhas, R., and Paterson, S.: Oil spill processes and models. Report to Research and Development Division, Environment Emergency Branch, Environmental Impact Control Directorate, Environmental Protection Service, Environment Canada, Ottawa, 1979.

Mackay, D., Paterson, S., and Trudel, B.: A mathematical model of oil spill behaviour, Report to Research and Development Division, Environment Emergency Branch, Environmental Impact Control Directorate, Environmental Protection Service, Environment Canada, Ottawa, 1980.

Nittis, K., Perivoliotis, L., Korres, G., Tziavos, C., and Thanos, I.: Operational monitoring and forecasting for marine environmental applications in the Aegean Sea, Environ. Modell. Softw., 21, 243-257, 2006.

Noye, J.: Numerical methods for solving the transport equation, Num. Model. Application Mar. Syst., 145, 195-229, 1987.

Pedlosky, J.: The buoyancy and wind-driven ventilated thermocline, J. Phys. Oceanogr., 16, 1077-1087, 1986.

Pinardi, N. and Coppini, G.: Preface "Operational oceanography in the Mediterranean Sea: the second stage of development", Ocean Sci., 6,263-267, doi:10.5194/os-6-263-2010,2010.

Pinardi, N., Allen, I., Demirov, E., De Mey, P., Korres, G., Las-caratos, A., Le Traon, P.-Y., Maillard, C., Manzella, G., and Tziavos, C.: The Mediterranean ocean forecasting system: first phase of implementation (1998-2001), Ann. Geophys., 21,3-20, doi:10.5194/angeo-21-3-2003, 2003.

Pollani, A., Triantafyllou, G., Petihakis, G., Nittis, K., Dounas, C., and Christoforos, K.: The Poseidon operational tool for the prediction of floating pollutant transport, Mar. Pollut. Bull., 43,270278, 2001.

Price, J. F., Weller, R. A., and Schudlich, R. R.: Wind-driven ocean currents and Ekman transport, Science, 238, 1534-1538, 1987.

Pugliese Carratelli, E., Dentale F., and Reale F.: On the effects of wave-induced drift and dispersion in the Deepwater Horizon oil spill, in: Monitoring and Modeling the Deepwater Horizon Oil Spill: A Record-Breaking Enterprise, Geophys. Monogr. Ser., 195, 197-204, 2011.

Reed, M., Gundlach, E., and Kana, T.: A coastal zone oil spill model: development and sensitivity studies, Oil Chem. Pollut., 5,411-449, 1989.

Reed, M., Aamo, O. M., and Daling, P. S.: Quantitative analysis of alternate oil spill response strategies using OSCAR, Spill Sci. Technol. Bull., 2, 67-74, 1995.

Reed, M., Johansen, 0., Brandvik, P. J., Daling, P., Lewis, A., Fiocco, R., Mackay, D., and Prentki, R.: Oil Spill Modeling towards the Close of the 20th Century: Overview of the State of the Art, Spill Sci. Technol. Bull., 5, 3-16, 1999.

Risken, H.: The Fokker-Planck equation, Springer, 475 pp., ISBN 354061530X, 9783540615309, 1989.

Rohrs, J., Christensen, K. H., Hole, L. R., Brostrom, G., Driv-dal, M., and Sundby, S.: Observation-based evaluation of surface wave effects on currents and trajectory forecasts, Ocean Dynamics, 62, 10-12, 1519-1533, 2012.

Schreurs, P., Mewis, J., and Havens, J.: Numerical aspects of a Lagrangian particle model for atmospheric dispersion of heavy gases, J. Hazardous Materials, 17, 61-80, doi:10.1016/0304-3894(87)85042-2, 1987.

Shen, H. T., Yapa, P. D., and Petroski, M. E.: A Simulation Model for Oil Slick Transport in Lakes, Water Resour. Res., 23, 19491957, 1987.

Sibert, J. R., Hampton, J., Fournier, D. A., and Bills, P. J.: An advection-diffusion-reaction model for the estimation of fish movement parameters from tagging data, with application to skipjack tuna (Katsuwonus pelamis), Canadian J. Fish. Aquatic Sci., 56, 925-938, 1999.

Sobey, R. J. and Barker, C. H.: Wave-driven transport of surface oil, J. Coastal Res., 13, 490-496, 1997.

Sotillo, M., Alvarez Fanjul, E., Castanedo, S., Abascal, A., Menen-dez, J., Emelianov, M., Olivella, R., Garcia-Ladona, E., Ruiz-Villarreal, M., Conde, J., Geomez, M., Conde, P., Gutierrez, A., and Medina, R.: Towards an operational system for oil-spill forecast over Spanish waters: Initial developments and implementation test, Mar. Pollut. Bull., 56, 686-703, 2008.

Spaulding, M., Kolluru, V., Anderson, E., and Howlett, E.: Application of three-dimensional oil spill model (WOSM/OILMAP) to hindcast the Braer spill, Spill Sci. Technol. Bull., 1,23-35,1994.

Stohl, A.: Computation, accuracy and applications of trajectories - A review and bibliography, Atmos. Environ., 32, 947-966, doi:10.1016/S1352-2310(97)00457-3, 1998.

Stokes, G.: On the theory of oscillatory waves, Transactions of the Cambridge Philosophical Society, 8, 441-473, 1847.

Tompson, A. and Gelhar, L.: Numerical Simulation of Solute Transport in Three-Dimensional, Randomly Heterogeneous Porous Media, Water Resour. Res., 26, 2541-2562, 1990.

Wang, J. and Shen, Y.: Development of an integrated model system to simulate transport and fate of oil spills in seas, Sci. China Technol. Sci., 53, 2423-2434, 2010.

Wang, S., Shen, Y., Guo, Y., and Tang, J.: Three-dimensional numerical simulation for transport of oil spills in seas, Ocean Eng., 35,503-510,2008.

Woods, J.: Primitive equation modelling of plankton ecosystems, Ocean Forecasting, Conceptual Basis and Applications, edited by: Pinardi, N. and Woods, J., Springer-Verlag Berlin Heidelberg, 2002.

Wunsch, C.: The work done by the wind on the oceanic general circulation, J. Phys. Oceanogr., 28,2332-2340, 1998.

Zelenke, B., O'Connor, C., Barker, C., Beegle-Krause, C. J., and Eclipse, L. (Eds.): General NOAA Operational Modeling Environment (GNOME) Technical Documentation, US Dept. of Commerce, NOAA Technical Memorandum NOS OR&R 40. Seattle, WA: Emergency Response Division, NOAA, 105 pp., available at: http://response.restoration.noaa.gov/gnome_ manual, 2012.

Zodiatis, G., Lardner, R., Hayes, D., Georgiou, G., Pinardi, N., De Dominicis, M., and Panayidou, X.: The Mediterranean oil spill and trajectory prediction model in assisting the EU response agencie, in: Congreso Nacional de Salvamento en la Mar, Cadiz, 2-4 October, libro de actas, 535-547, 2008a.

Zodiatis, G., Lardner, R., Hayes, D. R., Georgiou, G., Sofianos, S., Skliris, N., and Lascaratos, A.: Operational ocean forecasting in the Eastern Mediterranean: implementation and evaluation, Ocean Sci., 4, 31-47, doi:10.5194/os-4-31-2008, 2008b.

Zodiatis, G., Lardner, R., Solovyov, D., Panayidou, X., and De Dominicis, M.: Predictions for oil slicks detected from satellite images using MyOcean forecasting data, Ocean Sci., 8,1105-1115, doi:10.5194/os-8-1105-2012, 2012.

Copyright of Geoscientific Model Development is the property of Copernicus Gesellschaft mbH and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.