Oil Spill Modeling: A Critical Review on Current Trends, Perspectives, and Challenges

: Several oil spill simulation models exist in the literature, which are used worldwide to simulate the evolution of an oil slick created from marine trafﬁc, petroleum production, or other sources. These models may range from simple parametric calculations to advanced, new-generation, operational, three-dimensional numerical models, coupled to meteorological, hydrodynamic, and wave models, forecasting in high-resolution and with high precision the transport and fate of oil. This study presents a review of the transport and oil weathering processes and their parameterization and critically examines eighteen state-of-the-art oil spill models in terms of their capacity (a) to simulate these processes, (b) to consider oil released from surface or submerged sources, (c) to assimilate real-time ﬁeld data for model initiation and forcing, and (d) to assess uncertainty in the produced predictions. Based on our review, the most common oil weathering processes involved are spreading, advection, diffusion, evaporation, emulsiﬁcation, and dispersion. The majority of existing oil spill models do not consider signiﬁcant physical processes, such as oil dissolution, photo-oxidation, biodegradation, and vertical mixing. Moreover, timely response to oil spills is lacking in the new generation of oil spill models. Further improvements in oil spill modeling should emphasize more comprehensive parametrization of oil dissolution, biodegradation, entrainment, and prediction of oil particles size distribution following wave action and well blow outs.


Introduction
Due to the worldwide increase in oil/gas demand and the dwindling in onshore reserves, offshore oil/gas production has significantly increased its potential since the 1990s. In parallel, oil transportation technology evolved at the same pace as the oil production industry, with supertankers and pipelines of crude and product oil crossing the oceans. Offshore oil production and transportation threatens marine ecosystems with spillages, associated with immense environmental, social, and financial impacts, with the long-term effects being felt for decades [1].
Oil spills occur after natural releases, oil transportation, oil drilling and accidental collision or sinking of oil tankers, failures in pipelines and oil rigs, etc. Small spills are easier to handle effectively with existing technology. For this reason, size matters in oil spills and large spills are more important. On the other hand, recent studies have shown that major spill incidents have been fewer in number (Table 1), the broad public exhibits a "memory" on the major spills, but generally remains unaware that minor spills happen daily. Based on this figure, over 80% of the incidents recorded since 1970 were small spills (<7 tons). Following data from the European Space Agency (ESA), approximate losses of Table 1. Number of medium (7-700 tonnes) and large (>700 tonnes) spills per decade from 1970 to 2020 (source: [3]). Every time an oil spill occurs, the public loses faith in authorities and oil companies' capacity to implement preparedness and response decisions to mitigate impacts [4]. The severity of impacts typically depend on the quantity and type of oil spill, the ambient conditions, and the sensitivity of organisms and their habitats to the oil [2]. When crude oil is spilled on the sea, an oil slick is formed, i.e., a thin oily layer floating on the sea surface, affected by the large-scale advective processes dominated by currents, winds, and waves leading to center of mass slick transport (order of tens to hundreds of meters per day), and the slow, low-scale, diffusive processes reshaping the slick (order of centimeters to meters per day) responsible for modifying contaminants' concentration. The time scales and relative importance of the processes depend on spill-specific and environmental factors such as the quantity of oil spilled, the oil's initial physico-chemical characteristics, and meteorological and sea state conditions. In parallel, a series of natural, complex, and self-competing processes, referred to as "oil weathering processes" (OWPs), tend to degrade the slick [5,6]. As hydrocarbons are non-conservative pollutants, their physicochemical characteristics change over time as a result of OWPs. Spreading, evaporation, dispersion/diffusion, emulsification, and dissolution are the most crucial OWPs, acting at the early stages of the oil spill, while photo-oxidation, biodegradation, and sedimentation act in the longer term and determine the ultimate fate of the oil spilled. Oil density and viscosity are the parameters mostly altered by the OWPs after spillage [7].
Oil pollution may not only occur on the sea surface, but also in deeper waters, leading to even more extensive environmental impacts. The on-going exploitation of deep-water oil reserves and the installation of pipelines at high water depths increase the risks of accidental oil release from well blowouts and pipeline ruptures [8]. Major deep-water oil spill accidents caused by such occurrences are the 2010 Deepwater Horizon oil spill in the Gulf of Mexico, discharging approximately 492,000 to 627,000 tons of oil, and the 2011 Penglai 19-3 oil field spill in the Bohai Sea, China, discharging~200 tons.
Several oil spill models are extensively used at the global level to simulate the evolution of an oil spill. These models may range from simple vector-based calculations, such as the DHI oil spill model within MIKE [9], to the modern, new-generation, operational, three-dimensional (3D) numerical models, coupled to meteorological, hydrodynamic, and wave models, forecasting in high-resolution and with high precision the transport and fate of oil [1]. The simulation of the transport and fate of an oil spill at sea, appraising the physicochemical processes that occur between the oil phase and the water column, forms the basis for the evaluation of the engendered environmental, social, and economic impacts [10]. Thus, this study presents a state-of-the-art review on oil spill processes and their parameterization, and offers a critical comparison among the widely used oil spill models, in terms of their capacity to simulate the oil released from surface or submerged sources, the capacity to assimilate real-time field data to initiate model execution and correct the modelled forecasts in space and time, and the evaluation of uncertainty in the produced predictions. All of the above are crucial for the timely, efficient, and costeffective response to oil spills and should be considered in the real-time management of such incidents. Finally, this work will provide technical recommendations and will propose potential advancements and improvements in oil spill modeling. The paper is structured as follows: in Section 2, the physical transport and the physico-chemical OWPs are described; in Section 3, the state-of-the-art oil spill models are critically reviewed; Section 4 presents a critical comparison of existing oil spill models, while Section 5 proposes technical recommendations and modeling improvements. The present work expands and updates similar critical review papers, such as [10][11][12][13][14][15].

Oil Physical Transport and Weathering Processes
The behavior of an oil spill in the marine environment depends on a series of physical, chemical, and biological processes that are largely determined by both the properties of leaked oil and the environmental, hydro-meteorological conditions (wave, winds, currents, solar radiation, etc.), and discharge characteristics (instantaneous/continuous, surface/deep-water). The fate and behavior of an oil spill can be influenced by the physicochemical oil weathering processes: oil spreading, evaporation, emulsification, dissolution, photo-oxidation, biodegradation, and sedimentation, and the physical transport processes, like transport and turbulent mixing, dispersion, and resurfacing (see Figure 1). J. Mar. Sci. Eng. 2021, 9,181 3 of 39 models, in terms of their capacity to simulate the oil released from surface or submerged sources, the capacity to assimilate real-time field data to initiate model execution and correct the modelled forecasts in space and time, and the evaluation of uncertainty in the produced predictions. All of the above are crucial for the timely, efficient, and cost-effective response to oil spills and should be considered in the real-time management of such incidents. Finally, this work will provide technical recommendations and will propose potential advancements and improvements in oil spill modeling. The paper is structured as follows: in Section 2, the physical transport and the physico-chemical OWPs are described; in Section 3, the state-of-the-art oil spill models are critically reviewed; Section 4 presents a critical comparison of existing oil spill models, while Section 5 proposes technical recommendations and modeling improvements. The present work expands and updates similar critical review papers, such as [10][11][12][13][14][15].

Oil Physical Transport and Weathering Processes
The behavior of an oil spill in the marine environment depends on a series of physical, chemical, and biological processes that are largely determined by both the properties of leaked oil and the environmental, hydro-meteorological conditions (wave, winds, currents, solar radiation, etc.), and discharge characteristics (instantaneous/continuous, surface/deep-water). The fate and behavior of an oil spill can be influenced by the physicochemical oil weathering processes: oil spreading, evaporation, emulsification, dissolution, photo-oxidation, biodegradation, and sedimentation, and the physical transport processes, like transport and turbulent mixing, dispersion, and resurfacing (see Figure 1).   Spreading refers to the creation of a thin film, expanding over the sea surface, as soon as oil is being released [16,17]. Spreading algorithms in oil spill models provide an estimate of the spill thickness or surface area, used for modeling of many transport and fate processes such as evaporation, dispersion, and emulsification. Spreading rate and oil spill thickness depend on the sea surface temperature, oil viscosity, and density [18]. The most widely-used spreading algorithms have been developed by Fay [19,20] and Hoult [21]. The theory of gravitational spreading against viscous resistance is also followed in the Mackay's fate algorithms [22][23][24], modified versions of which are widely used in operational oil spill models (e.g., MEDSLIK, MEDSLIK-II). Advanced oil spreading algorithms consider processes such as wind shear stresses [25], turbulent mixing and wave breaking [26], and shear spreading. These processes may result in the break-up of the slick into patches and the dispersion and partial resurfacing of oil droplets [27], as well as into natural entrainment [28]. Some recent modifications and improvements in spill spreading estimation appear in the literature (e.g., [25,27,28]). The studies of Korinenko and Malinovsky [29] have shown that at different wind speeds the slicks have an elliptical shape and are oriented in the direction of air flow and that strong winds lead to an increase in the speed of spreading the slick along the main axis. Geng et al. [26] studied the effect of waves on the movement of oil droplets, illustrating that small eddy diffusivities decreasing rapidly with depth result in large horizontal spreading and vice versa. The work of [26] suggests that two-dimensional transport models could be overestimating the spreading of oil. The association of spreading with dispersion seemingly better illustrates the recognized physics of the dispersion process, once the initial gravity-viscosity spreading is accomplished [14]. Generally, spreading is a process with specific model limitations, as it depends on oil characteristics and ocean state, and existing algorithms only partially approximate the actual surface area of real spills. As oil is weathered, it is unevenly distributed into streamers and patches due to wave action and Langmuir circulation [25]. A rigorous solution to the problem requires sea state and oil data that might not be available in the initial stage of an operational spill response. Simecek-Beatty and Lehr [16] used Langmuir circulation models to approximate the merging of oil streaks and modify existing oil spreading parametrizations by estimating a spreading surface area correction due to Langmuir effects. Their model has been validated with measurements from the 1990s North Sea field experiments, and as it requires limited data, it could be successfully incorporated into operational response oil spill models. Another source of uncertainty is that, for computational purposes, oil spill models divide the slick into Lagrangian elements (LEs) or particles and track their movement, which does not directly provide an oil concentration or thickness at specific locations. Each model treats this with a different approach; in ADIOS2 [25] for example, each LE, representing a changing volume of oil, constitutes the center of a Thiessen polygon with a surface area relative to the local density of LEs, allowing the estimation of a variable local thickness, based on the polygon area and oil volume. The approaches followed by Lagrangian oil spill models to compute oil surface area or thickness adds further uncertainty in the spreading estimation.

Evaporation
Evaporation takes place when the volatile elements of the oil diffuse from the oil and entrain the gaseous stage, while the heavier components of oil remain at sea [10,30]. Evaporation removes most of the volatile fractions of oil to the atmosphere within a few hours, leading to the reduction of oil toxicity in the marine environment [7,10]. However, these compounds are transferred to the atmosphere and in some cases (e.g., large spills close to densely-populated areas), the effects of evaporation might be more toxic [31]. On the contrary, the viscosity of remaining "weathered patches" increases [30], leading to severe physical and chemical effects on the marine environment [10]. The most widelyused analytical method to assess the rate of evaporation is based on the work of Stiver and Mackay [32]. They assessed oil evaporation by means of a mass-transfer coefficient, expressed as function of wind speed, oil spill coverage, oil vapor pressure, and sea surface temperature. This parameterization has been included in the ADIOS1 model [33,34], developed by the National Oceanic and Atmospheric Administration (NOAA). However, this method treats the oil as a uniform element, whose features alter when the slick weathers, a procedure that may decrease the precision in the estimation of oil evaporation rates and is only valid for hydrocarbons with approximately linear distillation curves [25]. Oil is actually a complicated mixture of a large number of different types of chemical compounds; therefore it is vital to differentiate among the various chemical groups, in order to accurately estimate evaporation. A more elaborated and accurate model is the pseudo-component evaporation model of Jones [35]. In the pseudo-component technique, petroleum is assumed to constitute of a comparatively limited number of discrete, non-interacting components (pseudo-components or PCs). Each pseudo-component is handled as a single item with relative vapor pressure and relative mole fraction and molecular weight, and the total evaporation rate of the slick is the sum of the individual rates. A modified version of this model is incorporated in ADIOS 2 [25], while similar approaches of pseudo-component evaporation models are used in other models like in OSCAR [36,37] and SIMAP [38,39]. Contrary to the widely-used Mackay well-mixed boundary layer evaporation model, Fingas [40][41][42] suggested a different approach for oil evaporation modeling. Fingas [41,42] argued that oil evaporation is limited by the oil diffusion and therefore it is not an issue of the wind action on the slick thickness, but on the contrary, oil temperature is the main factor determining the evaporation rate.

Emulsification
Emulsification is the process by which water is being mixed into the oil. This water-in-oil emulsion in the form of suspended small droplets is often referred to as "mousse" [7,10,[43][44][45][46]. It occurs due to wave breaking, inducing sea surface turbulence, while oil composition, temperature, and viscosity play a significant role in the process [11,[47][48][49]. As oil viscosity increases, a higher amount of oil emulsifies and this additionally disrupts the rate of evaporation. In parallel, the rate of emulsification expands with increasing wind speed and turbulence at the sea surface [30]. Emulsified droplets may remain in the water column for longer time periods (from months to years). The main effect of emulsification is that it creates an emulsion of considerably increased viscosity, compared to the oil initially spilled, resulting in serious implications for treatment methods. Another important negative effect of emulsification is that it increases the volume of the slick; this means that the cleanup costs are greatly increased. Thus, emulsification is a process with specific model limitations and a crucial role on the impact assessment and response in oil spill modeling. A simple emulsification algorithm has been developed by Mackay [22,23], modified versions of which are currently included in several oil spill models (e.g., MEDSLIK, MEDSLIK-II, SIMAP). A literature overview of emulsification algorithms is given in [50][51][52][53]. Fingas [54,55] and Fingas and Fieldhouse [5,56] introduced a stability index (SI) according to density, viscosity, and type of oil in order to classify the emulsification tendency of oil [12]. In addition, SINTEF's (Selskapet for INdustriell og TEknisk Forskning) data-based oil weathering model (OWM) [57] can simulate emulsification quite well for certain types of hydrocarbons, for which laboratory or field data exist, by interpolating available data sets. However, a reliable emulsification forecasting algorithm based on environmental conditions and oil properties is not currently available to be incorporated into oil spill models.

Dissolution
Oil contains very small amounts of soluble compounds (<1 mg/L), which may dissolve in water, but is still considered an important process, since the lower molecular weight aromatic hydrocarbons (monoaromatic and polynuclear aromatic hydrocarbons, MAHs and PAHs), which are both highly volatile and soluble, are the most toxic elements of oil to aquatic organisms. Therefore, dissolution plays a significant role on environmental impact assessment and on response support models. Oil can dissolve in the water column from the surface slick or from dispersed oil droplets [13,58]. Dissolution and evaporation are two competitive processes [59], although evaporation exhibits faster rates and affects larger parts of the spill. Hydrocarbon components of lower molecular weight are highly soluble in seawater and relatively more volatile. Examples include the light hydrocarbons of benzene and toluene, which can dissolve within a few hours [7,13,58]. Generally, dissolution is significant when evaporation is low [30], therefore dissolution is substantial for subsurface oil spills and dispersed oil droplets. This is due to the lack of atmospheric exposure and the higher available oil surface area per unit of volume [12]. The algorithm developed by Mackay [60] is usually applied in oil spill modeling for estimating dissolution from the surface slick. This treats dissolution as a mass flux connected to oil solubility and temperature. Considering dispersed oil, dissolution is usually handled as a mass flux across the surface area of a droplet [60]. In SIMAP [38,39], the pseudo-component approach is followed for modeling oil weathering processes, including evaporation, therefore the dissolution and the toxic effects of lower molecular weight induced by the aromatic compounds to ecosystems can be more accurately estimated.

Photo-Oxidation
Photo-oxidation occurs when oil under the influence of the sunlight generates polar, water soluble, oxygenated compounds [11]. The process depends on the type of oil [1] and on the thickness of the oil slick [42]. Thick slicks may partially oxidize, generating tar balls, which are accumulated in bottom sediments or leach off the coast long after a leak. Generally photo-oxidation has long been considered a very slow process, with thin oil films dissolving, even in bright sunlight, at rates lower than 0.1% per day [1]. Thus, photo-oxidation is considered unimportant over the first few days of a spill but becomes visible after a week or more [10,60,61]. Therefore, photo-oxidation is not contained in modern oil spill response models. However, studies following the Deepwater Horizon oil spill have indicated that under conditions of high UV light exposure, photo-oxidation can be fast, affecting a considerable fraction of the spilled oil, and can contribute to emulsification [61,62]. Moreover, photo-oxidation alters the physicochemical characteristics of the oil and its related elements, with the oxygenated parts being more polar, expanding the dispersibility and dissolution and ultimately changing the toxicity biodegradability of the oil [63,64]. Currently, the existing numerical models do not consider photo-oxidation, since limited knowledge exists on this process, parametric expressions are lacking, and the importance and rates of the process have not yet been fully studied. Kolpack [65] developed a formulation for the rate of photo-oxidation in terms of the extrapolation of laboratory works to open ocean slicks, although this concept is yet not validated [11].

Biodegradation
Biodegradation of oil by native microorganisms is one of the most significant natural processes that can attenuate the environmental effects of marine oil spills in the long term. A number of in-depth reviews on this process are cited in the literature [66][67][68][69][70]. The biodegradation rate of oil depends on the type of petroleum hydrocarbons, temperature, species of micro-organisms, and the availability of oxygen and nutrients [66,70,71]. Furthermore, the rate of oil biodegradation [72,73] increases with the available water-oil interface, which for dispersed oil droplets increases as the size of droplets decreases. The application of chemical dispersants to enhance biodegradation by increasing the interfacial region available for biological activity [1] is still debated, as studies during and after the Deepwater Horizon oil spill showed both enhancement and inhibition of biodegradation [74,75]. As knowledge is still limited on how dispersants affect microbial species and their ability to biodegrade oil [75], as well as on the effects of different dispersant-to-oil ratios and the chemistry of oil-dispersant mixtures [74], further in-depth studies are needed to evaluate the application of chemical dispersants as a response option.
Biodegradation was generally considered a long-term oil weathering process, having a significant impact only after the first seven days of a spill and with time scales reaching up to several months; however, studies following the Deepwater Horizon oil spill have shown that under specific conditions it can become a significant process at shorter time scales, within the first week of an oil spill [62]. Therefore, until recently, biodegradation was typically not included in operational oil spill models, and when it was considered, it was treated as a first-order decay process, depending only on oil composition (e.g., SIMAP [38,39]).
Although several studies have been conducted measuring biodegradation rates and modeling the kinetics of dissolved oil and dispersed oil droplets, under different environmental conditions (e.g., [76][77][78][79][80][81]), these have not yet been integrated into oil spill models. Of particular importance is the recent work of Brakstad et al. [82][83][84] on oil droplet biodegradation, taking into account the droplet size distribution. This research has been incorporated into NOAA's GNOME oil spill model, by including a new biodegradation algorithm, which is dependent on the surface area of droplets [73]. Work performed by Kapellos [85] and Kapellos et al. [86] on the effect of biofilm formation, including the development of a shrinking core model, has not yet found its way into operational oil spill models. A comparison of different modeling approaches for oil droplets biodegradation following a deep sea blowout is documented in [72], employing the TAMOC model [87], while the importance of initial oil droplet size distribution and biodegradation for the subsurface transport of oil spills is highlighted in the work of North et al. [88] using LTRANS. Generally, there is a need for a more realistic description of biodegradation kinetics in oil spill models, including oil composition, dispersed oil droplets-water interface, but also other important parameters that may limit biodegradation rate such as microbial population, biofilm formation, and availability of dissolved oxygen and nutrients, to enable a more accurate prediction and evaluation of possible bioremediation scenarios and risk assessment in the mid-and long-term. Such an attempt was recently carried out by modifying the MEDSLIK-II [89,90] model, adding modules describing biodegradation of oil dispersed or dissolved in the water column and improving existing oil transport and weathering subroutines. In this modified version of MEDSLIK-II, the pseudo-component approach has been adopted for simulating weathering processes [91,92]. Biodegradation of petroleum is modelled via Monod kinetics. The kinetics of oil droplets size reduction due to the microbe-mediated degradation at the water-oil particle interface are described by the shrinking core model (SCM) [76,93].

Sedimentation
Sedimentation of oil droplets occurs as a result of three processes: increased density of the entrained oil and surface slicks due to weathering processes; incorporation of fecal pellets by means of zooplankton or benthic organisms' ingestion; and oil adherence or flocculation and agglomeration with suspended particulate matter (SPM) aggregates (OSA) [1, 11,30,42,[94][95][96][97]. For this reason, offshore OSA is a vital process to limit the transport of oil to nearshore benthic areas [97]. Generally, sedimentation of oil causes several impacts on the marine environment and for this reason it is a fundamental process for biological impact analysis and response oil spill modeling. On the other hand, oil sediments in OSA processes may or may not fall to the bottom. Recent works indicated that interactions among oil and sediments are critical in the dispersion and degradation of oil spills [98]. In regions with high SPM concentrations, increased dispersion and removal of oil is accounted for due to ingestion and adhesion [30,36]. Several parameters (e.g., temperature, salinity, wave energy, and physio-chemical oil properties) may control the OSA formation [99,100]. Moreover, the properties and characteristics of sediments constitute a significant role in OSA formation [98,100]. Due to knowledge gaps in properly expressing the detailed dynamics of sedimentation in a quantitative parameterization scheme, data are limited [13,64]. Khelifa et al. [101][102][103] introduced a Monte Carlo scheme, in terms of collision concept between particles of oil droplets and SPM to simulate the formulation of oil-mineral aggregates [14]. Recently, a conceptual development of oilparticle coagulation capability was developed by Zhao et al. [104,105]. Furthermore, the new term MOSSFA (marine oil snow sedimentation and flocculent accumulation) was identified in 2013 [106], after the Deepwater Horizon (DwH) accident, in order to assess the procedures affecting the formation and fate of oil-associated marine snow [107][108][109]. A well-defined schematic diagram of the process of MOSSFA into the water column and its driving parameters is presented in Quigg et al. [109]. Finally, the MOSSFA process has not been incorporated into any existing operational oil spill model.

Dispersion
Dispersion occurs when the waves or other turbulent events break over the oil slick surface and generate droplets of several sizes into the water column [42,[110][111][112][113]. The large droplets resurface to their primary region while the smaller spread and diffuse into the water column [110]. The rate of natural dispersion is influenced by environmental frameworks (i.e., the sea state), but also by oil properties and spill characteristics (oilfilm thickness, density, viscosity, oil/water surface tension), developing rapidly with low-viscosity oils in the presence of breaking waves [14,114].
Mackay [60] and Mackay et al. [22] developed an early model of wave entrainment, based on the fraction of the sea surface subjected to dispersion and the fraction of the entrained oil containing fairly small droplets to be constantly dispersed in the water column. Such parameterization considers both oil properties and oil film thickness. The fractional area of the surface slick dispersed at each time step depends on the sea state and is parameterized proportionally to the square of the wind speed. The formulation of Mackay et al. [22] has proven to succeed only at moderate wind speeds [13,115]. Delvigne and Sweeney [110] and Delvigne [116] later developed empirical formulations based on the experimental investigation of natural dispersion due to breaking waves. These commonly used models are empirical relations of the entrainment rate as a function of the dissipation of wave energy per unit area, the fractional area of the sea surface enclosed through breaking waves, and the volume of oil entrained per unit of water volume. The formulations of Mackay et al. [22] and Delvigne and Sweeney [22,110], as well as their modifications, are widely used in operational oil spill models, like ADIOS [25,33], SIMAP [38,117], OSCAR [118][119][120][121], OILMAP [122], MEDSLIK [123][124][125], and MEDSLIK-II [89,90]. In OpenOil [112], the method of Li et al. [126] is followed, which is a modification of the formulation of Delvigne and Sweeney [110], parameterizing the entrainment rate via the dimensionless Weber (We) and Ohnesorge (Oh) numbers. More recent parameterizations of entrainment from wave breaking incorporated the effect of viscosity, density, and oil-water interfacial tension [127][128][129]. It should also be noted that the fractional area of the sea surface encased by breaking waves, used to describe the sea state, has specific model limitations and is subject to large uncertainty [25]. It is regularly parameterized via the wind speed (e.g., in [130], subsequently used in [110]). Currently, there are numerous formulations for the wave-breaking fraction (e.g., [130][131][132]). However, vast uncertainty still exists in the linkage between the wind speed and wave breaking areal fraction [112]. The algorithms that are currently employed in oil spill models for natural dispersion, do not handle the knowledge gap of the process well, by assuming wave-averaged Eulerian velocities or mean dissipation rates. Future models should include the wave spectrum and white capping to improve the dispersion parameterizations and droplet formation. This is expected to improve the estimation of dissolution and biodegradation in the water column, as these weathering processes are influenced by oil droplets formation. Such an approach could also improve surface processes such as evaporation and distinguish oil partitioned between evaporation and dispersion.

Resurfacing of Submerged Oil
Resurfacing of the entrained oil droplets has as a result the movement of oil between the sea surface and the water column. The submerged oil droplets increase by virtue of their buoyancy, which is forced by means of their droplet size and the density difference among oil and water. It quickly proceeds for larger oil droplets, although the small-scale droplets remain in the subsurface for an extended period of time and can only resurface when wave turbulence decreases [12,25,133]. Oil droplets resurfacing has been modeled by Tkalich and Chan [111] and used in OpenOil [112]. The final vertical velocity relies on the Reynolds number of the flow over the droplet, according to Stokes' law, for low Reynolds numbers, and an experimental definition for the larger Reynolds numbers. Entrainment of surface oil and the associated droplet size spectra for the submerged oil, naturally affect the estimation of the subsequent oil resurfacing [112,113,134]. Droplet size distributions of dispersed oil may be declared either as a number size distribution or as a volume size distribution [133]. Although [110] noticed a power-law number size distribution, current experimental research indicates that the droplet size distribution is better expressed via a lognormal distribution [127,133] or as two regimes with various exponents of power-law [133,135]. Identified droplet size distributions depend on oil type, sea state, oil weathered state, oil-water interfacial tension, and initial oil slick thickness [127,133]. The diameter of oil droplets, used via the droplet size distribution in oil spill models, directly affects oil droplets resurfacing through the calculation of the advective flux due to buoyancy [112].

Turbulent Mixing
Turbulent mixing moves oil and mixes it into the water column. While buoyancy moves oil droplets in one direction, turbulent mixing transfers oil particles upwards and downwards. It affects mainly smaller oil droplets, diminishing their opportunity to resurface [12,107,111,133]. This process has a significant role in the vertical interchange among the surface oil spill and the vertical layers of the water column [112]. The main source of uncertainty in oil spill simulations arises from uncertainties in the forcing of models, i.e., ocean circulation, wave and atmospheric coupled models, therefore reliable forecasts are essential for accurately determining the advective transport [15,136]. The volume of turbulent mixing is widely represented via an eddy diffusivity coefficient. The eddy diffusivity can be provided by ocean circulation models [137] or be approximated by the wind speed (e.g., [138]). Using eddy diffusivities provided by ocean circulation models, the exerted wind forcing is considered, as well as the advection and inertia of turbulence and buoyancy and inhibition through seawater stratification. When turbulent mixing levels are properly increased, the oil particles are maintained in the water column [127]. Novel ocean models with real-time data present details about the vertical currents, stratification, and turbulent mixing, providing more sensible particle transport representations [27,112,139]. Vertical mixing algorithms and parameterizations are provided by Galt and Overstreet [28], Röhrs et al. [112], and Nordam et al. [140]. In Lagrangian particle tracking oil spill models (e.g., MEDSLIK, MEDSLIK-II, OpenOil), the turbulent flux can be expressed via a random walk process, according to Visser [141]. From this perspective, a random vertical displacement is estimated for each particle (e.g., [112]).

Transport
Horizontal and vertical transport of oil spilled at sea are separate processes, which are vital for the circulation of oil spills in the sea water [142]. Horizontal transport includes spreading and advection while vertical transport involves vertical dispersion and wave entrainment, turbulent mixing, and resurfacing. Horizontal transport mainly depends on advection due to ocean currents, waves, and winds, while vertical transport has a crucial aspect, affecting the horizontal transport of oil slicks and generating a mixing layer at the top of the water column via breaking waves [112]. Wind resistance is generally considered to affect only the surface slick, while ocean currents and the wave-induced Stokes drift range according to depth and are therefore also important for the movement of discrete oil parcels in the subsurface [115,143]. Therefore, in order to simulate the transport and fate of oil spill, a well-described expression of surface slicks is demanded, together with the vertical distribution of submerged oil [112]. The common modeling technique, employed by nearly all oil spill models, to account for the effect of wind on the oil slick floating on the sea surface is to use a "wind factor" approach, i.e., the effect of wind will move oil at a certain fraction of the wind speed and at a certain angle to the wind direction [144]. However, there is considerable dispute as to what are the most effective options for the values of the drift factor and angle in combination with the sufficient vertical resolution of ocean forecasting models to resolve the vertical structure of the current flow, so that the motion of the surface layer is computed accurately. A comprehensive review on this is given in [145], while examples of the sensitivity of current depth in oil spill modeling are provided by [90]. Improvements in parameterization of wind drag have been introduced by [146]. In addition, in the majority of Lagrangian oil spill models (e.g., MEDSLIK, MEDSLIK-II, OpenOil) oil particles/parcels are assigned an advective displacement according to currents, wind and Stokes drift, and a diffusive displacement given by a random walk model.

Oil Spill Models-The State-of-Art
Oil spill models are numerical tools capable of (a) forecasting the trajectory of a spill, (b) estimating the time needed for the spill to reach specific areas of interest, and (c) assessing its state when it arrives at the modeled locations. The first two issues require accurate data on winds, currents, and waves in the broader area of the oil spill accident, while the third issue requires deep understanding and reliable algorithms of the oil weathering processes. Oil spill models may be used by authorities for contingency planning and emergency response to a crisis occurring due to accidental oil releases. Such planning, in conjunction to an oil spill model, may lead to the deeper understanding on the effects of oil weathering processes on oil spillage, at the surface and within the water column, and thus to improved methods to monitor and to clean it up [8].
Generally, oil spill models may be categorized in two types: Eulerian and Lagrangian. The first approach deals with the mass and momentum conservation equations applied to the oil slick or with a convection-diffusion equation. In this latter, the diffusive part of the equation illustrates the spreading of oil and the convective terms describe the advection of oil through currents and wind [145]. On the other hand, the Lagrangian models discretize oil slicks as a large number of particles advected by the merged result of winds, waves, and currents, but also being transported via dispersion. Researchers have shown that Lagrangian models are more appropriate for prompt simulations, when oil spill accidents occur, and consequently are easier, more efficient, and computationally more cost-effective than the Eulerian approaches [10,15,147]. The precision of the particle position will depend on the precision of its initial position, on the capacity of the coupled metocean models to guarantee reliable forecasting, and on the involvement of suitable physical mechanisms taking action on the tracked particles [145]. The transport and fate of oil particles is estimated by solving the general equation for an active tracer concentration, C = C → x , t , with units of mass per volume, mixed in the marine environment, and is widely referred to as the advection-diffusion-reaction equation:

∂C ∂t
where ∂ ∂t is the local time-rate-of-change operator, → U is the three-dimensional ocean current mean field, → K is the diffusivity tensor, which parameterizes the turbulent effects, and x , C, t represent the transformation of oil via the physico-chemical processes.
Overall, there are various models in literature aiming to predict the trajectory and fate of hydrocarbon spills occurring at the sea surface and the water column. The main inputs in each case are not only oil spill data, such as the type of oil and the initial location of spillage, but also metocean variables, such as the three-dimensional flow field, sea temperature, salinity and density profiles, atmospheric winds, and bathymetry. These latter data may be obtained from different operational oceanographic forecasting systems, such as CMEMS (Copernicus Marine Environmental Monitoring Service, http://marine. copernicus.eu/ [148]) and NOAA (National Oceanic and Atmospheric Administration). Furthermore, advanced oil spill models have the capacity to use satellite synthetic aperture radar (SAR) imagery to detect oil spills [149][150][151]. Generally, to run an oil spill model, the following data are required as a minimum: (a) oil spill scenario details, (b) oil properties, (c) metocean data, and (d) output requirements.

Surface Oil Spill Models and Blowout/Buoyant Plume Models
Oil spill models have largely emphasized surface spill behavior and tracking, thus requiring the highlighting in the construction of comprehensive blowout/buoyant plume models. Subsurface oil spillages may result from single well tapping in a unique tank or from a platform that connects multiple tanks and other sources of submarine spills, such as wrecks, accidents, or blow-ups [152]. On the grounds that the combined vertical and horizontal transport of oil occurs during a massive deep spill, the potential ecological and human effects are far more substantial and convoluted than in ordinary surface oil spills [153]. A comprehensive table comparing the weathering processes during surface and subsurface releases is given in [154].
Earlier research has demonstrated that an underwater oil spill is primarily controlled by three parameters: initial jet momentum, plume buoyancy, and ambient current and turbulence [8]. Therefore, the problem of modeling subsea releases tends to become a problem of any underwater outfall, determined by the dynamics of buoyant jets and oil plumes, describing the eventual fate and transport of the spill [155]. Latest achievements in underwater oil spill modeling are presented in a review paper from Yapa et al. [156].
Overall, several elements like the restrictions in accurately measuring the prevailing oceanographic conditions and the limited available data on released volumes, exact location, fluxes, etc., make the modeling of deep subsea releases more demanding than shallowwater and surface releases. In blowouts/buoyant jets, this is attributed mostly to the expanded interaction among oil and the water column, in advance of atmospheric exposure, the presence of strong ocean currents, high pressures and low temperatures near the seabed, the interaction with the sub-bottom rocky layers and submarine sediments, and the existence of high pressures and temperatures in oil and gas reservoirs [152][153][154][157][158][159]. A detailed description of the convoluted thermodynamic processes, which take place in the near-field, and the hydrodynamic processes in the far-field is presented, for example, in [159]. Several works are cited in the literature that identify various flow regimes, varying in the composition of individual droplets and the output geometry to the complete spray of a jet [160][161][162][163]. A schematic diagram aiding the understanding of deep-water oil spill fate and effects is introduced by Murawski et al. [153].
In these deep-water releases, highly affected by the hydrostatic pressure, jet breakup events may occur [159], simulated by two methods: the first involves experimental equations predicting a droplet size at the end of the dynamic breakup location [135,164,165], while the second approach addresses the competing physiochemical processes by adjusting decomposition and aggregation, permitting a dynamic result of the complete size distribution of bubbles and droplets via the jet splitting area [166][167][168][169]. Furthermore, as for the buoyant plume, the first technique includes integral models, resolving the cross-sectional averaged flow along the center line trajectory of the plume [170][171][172], and the second method involves the use of three-dimensional computational fluid dynamics models [159,173].

The New Generation of Oil Spill Models
Broadly speaking, state-of-the-art oil spill models produce not only oil spill predictions, but also the assessment of ambiguity of such forecasts, which is crucial and urgent for up-to-date, beneficial, and cost-effective responses. This uncertainty in the forecasting of oil transport and transformation arises mostly from uncertainties in the input fields (errors in initial conditions, environmental data, and in the predictions of metocean models), internal model dynamics (e.g., numerical scheme, parameterization of transformation processes), and sparse observational data [174,175]. Due to this large number of uncertain sources introduced in oil spill models, ensemble forecasts are important to improve the quality of predictions (e.g., [176,177]). ASTM (American Society for Testing and Materials) [178] has established a standard for oil spill models requiring uncertainty estimates for oil spill trajectory forecasts to support response operators. However, the methodology for uncertainty estimation is not well-established in oil spill models [179] and presents a field of future research. Nevertheless, NOAA's GNOME [180,181] model, for example, includes uncertainty algorithms regarding the perturbation of current and wind fields. De Dominicis et al. [176] used an ensemble of metocean models to improve oil trajectory forecasts with MEDSLIK-II model. Liubartseva et al. [179] introduced an uncertainty module in WITOIL Decision Support System, which includes MEDSLIK-II for oil spill forecasting, to automatically estimate prediction uncertainties related to the initial conditions of the spill, based on a parametric analysis methodology, employed in atmospheric pollution models [182]. Oil spread probability maps are produced as an indication of predictions uncertainty.
In parallel, state-of-the-art oil spill models use satellite SAR images/data to identify potential oil slicks and implement spill and drifter surveillance to improve slick forecasting. In detail, existing oil spill remote sensing techniques are presented in the review papers of Fingas and Brown [183,184]. The attention of the scientific community has been focused on enhancing 4D predictions by simulating oil spills backward in time to track the slick to its source [145]. These back-propagation approaches, when correlated with the operation of the AIS (automatic identification of ships) system, could track down the sources of world-wide oil spills.
CDOG (comprehensive deepwater oil and gas model) is a three-dimensional model, developed by Yapa and Li [206] and modified by Zheng et al. [186]. The model simulates the aspect of oil and gas released from deep water accidents [172,186]. Moreover, in CDOG, hydrate formation and disintegration, gas dissolution, non-ideal behavior of gas, and potential gas partition from the basic plume, in virtue of strong cross-currents, are connected to the jet/plume hydrodynamics and thermodynamics [172]. CDOG includes unsteady-state 3D fluctuation of ambient currents, density stratification, salinity, and water temperature [207]. Although CDOG has been implemented for response purposes, its main objective is research. Recently, the US government agencies (MMS (Minerals Management Service), NOAA) and oil companies have started using the CDOG model. OILMAPDEEP (deep water oil spill model and analysis system) [188][189][190][191]208] has been developed by Applied Science Associates (ASA) in order to estimate the fate and transport of subsea releases. OILMAPDEEP estimates the near-field plume characteristics and oil droplet size distributions for a specified release [189,208]. Oil droplet size distribution predictions are in accordance with the study of [171,209]. Moreover, the trap height and droplet size distribution are used as initial conditions for the SIMAP [38,39,210], which computes the transport and fate oil processes in accordance with the near-field buoyant plume [208,210]. The model provides simulations of both the near-field and far-field environment [127,191,208,211]. Fate processes are included for both gas and oil, however, the details of the modeling algorithms are unpublished and rely on a database of oil chemical features, according to ASA. The model provides a subsurface dispersant treatment module and a Lagrangian particle tracking module, incorporating 2D and 3D hydrodynamic model flow fields [190,191]. Output data contain plan and section views of plume, in-water, and on-surface model forecasting [189][190][191]. The model has global capacity and includes RPS (Rural Planning Services) ASA's own GIS.
SIMAP (integrated oil spill impact oil system) [38,39,208,210] also developed by ASA provides simulations of the three-dimensional trajectory, fate, and transport, as well as biological effects and other impacts of spilled oil and fuels [39,192]. Moreover, the model may be run in both stochastic and deterministic modes and includes a buoyant plume transition stage to the far field. The model has a Lagrangian particle tracking module in the far field [212]. It includes oil processes with specific model limitations, such as dissolution and sedimentation of oil, sinking, evaporation, dispersion, and spreading, complex oil and ice interaction, together with sediment and shoreline contamination. Some applications of the SIMAP model [39,192,212] include the environmental impact assessment of oil spills, hindcast/forecast simulations, natural wealth damage evaluation, contingency planning, environmental risk assessment, and cost-effective study. SIMAP has been validated against data of more than 20 large spills, such as the Exxon Valdez [38,39,213].
OSCAR (oil spill contingency and response model) is an advanced, three-dimensional model for planning and response to oil spills, developed by SINTEF [37,115]. It calculates the fate and effects of surface releases or blowout/buoyant plume of oil or gas [121]. The chemical fates sub-model allows multiple separate pseudo-components, which are transported across all environmental segments [37]. The transport and fate of oil spills at the surface are described not only by virtue of currents, winds, and turbulent diffusion, but also by means of oil-weathering algorithms, such as spreading, evaporation, natural dispersion, emulsification, dissolution, and volatilization. Moreover, in the water column, horizontal and vertical dispersion of entrained and dissolved hydrocarbons are represented via random walk approaches. Finally, the degradation and sedimentation processes of oil are described as first order decay process [121]. Essential elements of the model are SINTEF's data-based oil weathering model [214][215][216], the three-dimensional oil trajectory and chemical fates model [118], an oil spill battle model [120], and exposure models for fish and ichthyoplankton [119], birds, and sea mammals [217]. Overall, OSCAR has been used in oil spill risk assessment, as well as in response planning and operations [121]. The model has been applied for hindcast and forecast of accidental releases in locations such as the North Sea, the Baltic Sea, the Gulf of Mexico, and the Mediterranean basin [121]. In the UK, OSCAR is routinely used for operational forecasting of oil spills, forced by ocean circulation models such as the U.S. Navy global hybrid coordinate ocean model (HYCOM) or the Copernicus system and wind forecasts from NOAA's GFS (global forecast system) or CFS (climate forecast system).
OILMAP [122] has been developed by ASA as well as SIMAP, and both of them share the same code base. However, OILMAP is a three-dimensional oil spill response and contingency planning model. It deals with both surface and subsurface hydrocarbon releases and provides algorithms for oil spreading, evaporation, emulsification, entrainment, and oil-shoreline, oil bed, and oil-ice interaction [122,218]. The stochastic module predicts an extensive number of trajectories from a single site for producing probability statistics [218]. The distribution and mass balance of oil over time are simulated per type of oil spilled. The model has been applied in Dubai and Gulf region in 2006 [219]. It is used operationally by Oil Spill Response Limited (OSRL) in United Kingdom.
TAMOC (Texas A&M oil spill calculator) [193][194][195] is an open-source model, written in Python and Fortran, which simulates subsea oil spills and blowout plumes. Furthermore, its code is available for users in Github: http://github.com/socolofs/tamoc. It computes near-field plume dynamics, dissolution, particle tracking, transport of oil droplets, and phase equilibrium of hydrocarbons and incorporates an all-inclusive fate module. An extensive report of the TAMOC model and its mathematical background and equations are mentioned in the works of Gros et al. [193,195]. The oil fate and transport are expressed according to the formulation of McGinnis et al. [220] and as for jet and plume schemes, these are described by an integral model method [170,186,221]. A key feature of this model is the combination among the extended hydrodynamic behavior and the dynamic equations of motion, such as plume and intrusion formation. Finally, TAMOC has been validated via several experimental studies of bubble plumes, such as [221].
MOTHY (modèle océanique de transport d'hydrocarbures), developed by Météo-France [197], is a 3D Lagrangian pollutant drift model predicting the fate and transport of oil slicks on the ocean surface. MOTHY has been operational since 1994 and it has been used and validated during major real oil spill incidents, such as the Erika [197,222] and the Prestige [223]. The mixed layer is expressed via a combination of a shallow water model relative to the wind and the atmospheric pressure, in cooperation with a well-described turbulent viscosity model, while hydrodynamics are provided by CMEMS (Copernicus Marine Service) Med MFC (Marine Fisheries Commission) models and wind forcing provided by European Center for Medium-Range Weather Forecasts (ECMWF) [224]. National, higher-resolution ocean forecasting systems nested in CMEMS Med MFC are used in sev-eral cases to resolve coastal scale processes in various areas of the Mediterranean. The water column is described by a continuous profile from surface to bottom [224,225]. Turbulent diffusion is modeled via a three-dimensional random walk scheme [226]. This oil spill model provides some additional capacities: beaching, sedimentation, and backtracking, while pollutants can be either oil or floating objects [225].
OILTOX is a Lagrangian oil spill model [200] adapted to the Black Sea environment. It includes hundreds of oil types that are transported via the Black Sea and their fundamental physical-chemical features. This model simulates oil transport and fate according to [200,227] in five phases: oil-on-surface, oil-in-water, oil-on-bottom, oil-on-suspended sediments, and oil-at-shoreline. The model incorporates the basic transport and weathering processes, such as spreading by virtue of gravity and surface tension, advection due to wind and surface currents, evaporation, emulsification, oil-shore interaction, wave entrainment, resurfacing of entrained oil, and sedimentation [200]. Moreover, the model incorporates horizontal and vertical turbulent diffusion processes, which are represented by means of a random walk method.
The MOHID (Modelo Hidrodinâmico) Lagrangian oil spill model [201] is a sub-model of the MOHID water modeling system [228], developed by the Technical University of Lisbon. The movement of the tracers is caused by the surface flow field, the atmospheric winds, the spreading velocity from the dispersion module, and a randomly produced velocity via a random walk approach. MOHID Lagrangian transport module includes the following features: oil transport in water column, sedimentation, and beaching; oil weathering processes such as evaporation, dispersion, entrainment, sedimentation, dissolution, emulsification, and dispersion; and Eulerian concentration result.
The POSEIDON OSM is an oil spill model generated by the Hellenic Centre for Marine Research (HCMR), implemented and operational in the Aegean and Ionian Seas ( [202,229] since 2000. It is a completely 3D oil spill model with the capacity not only to predict the transport, spreading, and weathering of the oil particles in the 3D space, but also to provide various oil weathering processes, such as evaporation, emulsification, beaching, and sedimentation [174,225]. Oil advection and dispersion is illustrated via a vast number of particles, each of which expresses a group of oil droplets of similar size and composition [174]. Oil transport is calculated using two modules: the circulation module and the wind generated wave module [174,225]. Moreover, the horizontal movement considering advection and the vertical transport of the oil are described through the results of the POSEIDON ocean forecasting system [225]. Stokes drift is also provided by the coupled wave model of POSEIDON ocean forecasting system [174,225]. Recently, some additional characteristics were integrated in the model via a dedicated web-based application (https://poseidon.hcmr.gr/components/forecasting-components/oil-spillmodel [230]), where the user can determine the parameters of a real or hypothetical scenario and submit this event to the system, receiving the model output in Google Earth file format for a more real-time geospatial simulation. GNOME (general NOAA operational modeling environment) is an oil spill model that predicts the fate and transport of pollutants and oil movement caused by winds, currents, tides, and spreading [180]. GNOME was developed by the Hazardous Materials Response Division and it is an open-source model, freely available in Github: https: //github.com/NOAA-ORR-ERD [231]. The model is publicly available for use by the broader academic, response, and oil spill planning communities. GNOME provides the following elements: 3D particle transport, able to work with virtually any hydrodynamic model and measured field data, 1, 2 or 3rd order Runga-Kutta algorithm, with droplet rise velocity depending on density and droplet size; "leeway" wind surface transport: randomly adjustable with various user-adjustable values, providing a configurable downward spread; open-source code; backward running; oil weathering algorithms from the integrated open source ADIOS oil database, which is currently getting updated, with a beta version available for testing at: https://adios-stage.orr.noaa.gov [232]; sea ice interaction according to ice concentration and velocity; shoreline interaction (beaching) with configurable half-life based re-float; includes the TAMOC deep-water blowout model; comprehensive script for stochastic analysis and other batch processing; configurable for use on other drifting objects, such as for SAR and marine debris; integrated response options calculator (ROC) to evaluate the performance of spill response systems such as skimming, burning, and application of chemical dispersant; the PyGNOME, a Python setting, to build the web GNOME interface to the model, that performs batch processing and testing; and includes a GIS system for the model outputs visualization. In addition, GNOME is extremely configurable and tunable to adjust to field conditions and it can be driven via numerous data: measured point data, met models, and hydro models with a variety of meshes (structured, triangular). Finally, it has been used to support spill response for oil spills in the USA for almost twenty years. As GNOME can be integrated with any ocean circulation and meteorological model providing forecasts at different file formats, as well as observational data, NOAA has developed the GNOME Operational Oceanographic Data Server (https://gnome.orr.noaa.gov/goods) [233], a publicly available system to provide access to all the driver models and data sources available. Another important feature of the operational use of GNOME is the assimilation of available observations of oil spill locations in each forecasting cycle. Model parameters are fine-tuned to match the observations, and subsequently a new forecasting cycle and analyses are produced. Observational data assimilation improves the accuracy of forecasts for response authorities.
The OILTRANS particle transport model [203] is based on the LTRANS v.2 particle transport model, developed by North et al. [234]. The oil fates module of OILTRANS simulates the transport, fate, and oil weathering processes coupled to state-of-the-art operational metocean model [203]. The model provides the oil fate processes of spreading, advection, diffusion, evaporation, emulsification, and dispersion in order to estimate the horizontal movement of surface oil slick, the vertical entrainment of oil into the water column and the oil mass balance [203]. OILTRANS can be applied in any ocean or coastal field. The minimum data required are: bathymetric data, tidal current, and wind fields, together with information on the location, quantity, and type of spilled oil [235,236]. OILTRANS has been used, for example, for an accidental release in the Celtic Sea in February 2009 [203].
OSERIT, oil spill evaluation and response integrated tool, is an oil spill model that is capable of predicting the 3D drift and the fate of an oil spill at the surface and into the water column [237]. It contributes to the forecasting service of EMSA CleanSeaNet and has been used in the North Sea. The Lagrangian module expresses the independent movement of each parcel due to the winds, currents, and waves. Furthermore, OSERIT contains the buoyancy effect, turbulent diffusive transport, vertical dispersion of oil from surface to the water column, horizontal spreading, and beaching [237]. Moreover, it is able to calculate the drift of chemically dispersed oil and forecasts oil weathering processes, such as evaporation and emulsification, and their effects on oil features. Biodegradation and oil sedimentation are not included in OSERIT. The oil database of OSERIT is based on the oil types included in the ADIOS database [54,238].
BLOSOM (blowout and spill occurrence model) has been developed by the National Energy Technology Laboratory (NETL) of the USA. The model is written in Java programming language [239] and it is an extensive modeling suite that displaces the fate and transport of both subsurface oil blowouts and surface spills [240]. Moreover, this model is developed to predict offshore oil spills resulting from deep water (>150 m) and ultra-deepwater (>1500 m) well blowouts [155]. The model simulates oil spills from source to final fate and degradation stage. BLOSOM is flexible in its construction and utility from using it for basic particle tracking to applying advanced weathering modules and modules for jet/plume modeling [155]. BLOSOM supports risk evaluation and provides a comprehensive tool for response planning. It is designed to handle deep-water blowouts, such as Deepwater Horizon [241]. The jet plume element of BLOSOM has been assessed via experimental studies, which took place in the North Sea [242,243]. BLOSOM integrates various oil types from the ADIOS oil library [25].
Delft3D-PART, developed by Deltares, is a module of the Delft3D modeling suite that estimates the transport and simple water quality processes via a particle tracking method, implementing the 2D or 3D flow data by the Delft3D-FLOW (hydrodynamic module) [244]. Some test cases of the model are included in [245,246]. The particle tracking scheme follows a random-walk approach, referred to as the "Monte Carlo method" [247]. Since the simulated behavior is stochastic, the number of particles is limited [248], however it is the only stochastic model in the full Delft3D modeling suite. Moreover, Delft3D-PART provides two modules: (a) the tracer module, simulating conservative or first order decaying substances, and (b) the oil spill module, simulating oil spills with floating and dispersed oil fractions. Furthermore, the processes that are involved in the oil module are advection of floating oil via wind and ocean currents, the dispersion of oil caused by waves, evaporation of floating oil, emulsification, sticking of petroleum to the coastline or seabed [244]. In addition, this oil module contains variations of oil features (density, viscosity, water content) on account of the above processes. It includes vertical dispersion for well-mixed systems and horizontal dispersion resulting from turbulence, being enhanced in time in accordance to the turbulence theory [244]. In Delft3D-PART it is feasible to illustrate a maximum of 30 different oil types.
MEDSLIK-II [89,90] is based on its precursor MEDSLIK oil spill model. It is an opensource oil spill model for surface oil spills in the marine environment. It is designed to forecast the transport and weathering of an oil slick and to express the displacement of a floating particle, using a Lagrangian formalism, in conjunction with an Eulerian ocean circulation model. Moreover, MEDSLIK-II predicts the transport of the surface slick due to the water currents and the wind. Oil particles are also dispersed by turbulent fluctuation elements [23], being formulated via a stochastic approach [13,249] using a random walk scheme [250,251]. For the Stokes drift parameterization, MEDSLIK-II uses the experimental Jonswap wave spectrum in terms of wind speed and fetch [252], while in MEDSLIK the forecasting wave parameters of SWH and wave period are used.
The necessary oil spill data to define initial conditions include the oil spill location, time and areal coverage of the spill, rate and duration of spillage, type of oil, and age of the oil spill from initial arrival in the sea. These data can be simply included to MEDSLIK-II via satellite monitoring systems [184]. In addition to oil spill data, MEDSLIK-II requires as input the wind field, the sea surface temperature, and the three-dimensional sea currents. MEDSLIK-II is closely coupled in terms of input format with atmospheric fields, provided by the European Center for Medium-Range Weather Forecasts (ECMWF) and with oceanographic fields (currents, temperature, salinity, and density), provided by CMEMS Med MFC. For several coastal applications in the Mediterranean basin, local high-resolution forecasting systems, which are nested in CMEMS Med MFC products, are used for providing met ocean forcing for operational oil spill forecasting with MEDSLIK-II. MEDSLIK-II produces as output the oil properties evolution and the position of the surface, dispersed oil, and of the oil arrived on the coasts. Furthermore, MEDSLIK-II calculates the mass balance components of the oil, with respect to time, providing time-effective tracking of oil weathering processes [174,179]. Figure 2 includes a schematic diagram of Medslik-II model with input and output data.
A comprehensive description of MEDSLIK-II with the elaborate mathematical concept and the corresponding basic parameters of the model is given by De Dominicis et al. [89]. Moreover, the model has been validated with in situ data, surface drifters data, and with satellite data [90,176]. Interesting applications of the model are included in [176,177,179,[253][254][255][256][257]. MEDSLIK-II is used operationally in the Mediterranean region, allowing also support to REMPEC [258] for oil spill emergencies in the entire Mediterranean basin. cept and the corresponding basic parameters of the model is given by De Dominicis et al [89]. Moreover, the model has been validated with in situ data, surface drifters data, and with satellite data [90,176]. Interesting applications of the model are included in [176,177,179,[253][254][255][256][257]. MEDSLIK-II is used operationally in the Mediterranean region, al lowing also support to REMPEC [258] for oil spill emergencies in the entire Mediterranean basin. The model has the following extra features: it incorporates a built-in oil database (from REMPEC [259]) with over 220 oil types which are widely used in the Mediterranean and the Black seas; it has been applied to forecast oil spill fate and transport during numerous emergency cases in the Mediterranean Sea (Lebanese oil pollution crisis in 2006, Und Adriyatik in 2008, and the Costa Concordia emergency in 2012). The sensitivity of the oil spill predictions to several model parameterizations is examined and the outputs are validated by means of surface drifters, SAR (synthetic aperture radar), and optical satellite images.
Due to a vast range of parameters that handle the oil movement and transformation in MEDSLIK-II, recently, implemented straightforward and effective algorithms to evaluate uncertainties may resulting from the initial oil spill conditions [179], ocean currents, and winds [174].
The oil drift module OpenOil is based on the open source [260], python-based, trajectory framework of OpenDrift [205], and it is a newly-integrated oil spill transport and fate model [112]. A conceptual flow diagram of the OpenDrift model is included in Figure 3. OpenOil (Figure 4) has been implemented operationally in Norway, as an oil spill contingency and search and rescue model [261,262] and for drifter and oil slick observations in the North Sea [112,263]. MET Norway uses in-house high-resolution ocean circulation and meteorological models as forcing for providing operational oil spill forecasts with OpenOil. If needed, coarser resolution forecasts from CMEMS for hydrodynamics and ocean state and NOAA's GFS for wind fields can be also used. This model integrates algorithms with several physical processes, such as wave entrainment of oil [126], vertical mixing by virtue of oceanic turbulence [126,141], resurfacing of oil on account of buoyancy [111], and emulsification of oil properties [25,126]. Resurfacing is parameterized based on oil density and droplet size by means of Stokes Law, and for this reason the model's physics are very sensitive to the specification of the oil droplet size distribution [261,262]. In OpenOil, the oil properties are obtained from the ADIOS Oil Library. The ADIOS oil database [25] is also open-source, written in Python, and contains measured properties of almost 1000 oil types across the world [205]. In contrast to the above, dissolution, which is an important oil weathering process for blowout/buoyant plume models, is not yet implemented.
ocean state and NOAA's GFS for wind fields can be also used. This model integrates algorithms with several physical processes, such as wave entrainment of oil [126], vertical mixing by virtue of oceanic turbulence [126,141], resurfacing of oil on account of buoyancy [111], and emulsification of oil properties [25,126]. Resurfacing is parameterized based on oil density and droplet size by means of Stokes Law, and for this reason the model's physics are very sensitive to the specification of the oil droplet size distribution [261,262]. In OpenOil, the oil properties are obtained from the ADIOS Oil Library. The ADIOS oil database [25] is also open-source, written in Python, and contains measured properties of almost 1000 oil types across the world [205]. In contrast to the above, dissolution, which is an important oil weathering process for blowout/buoyant plume models, is not yet implemented.

A Comparative Assessment of State-of-the-Art Oil Spill Models
Accurate prediction of transport and fate of spilled oil is challenging due to the difficulties in understanding the oil behavior and its interaction with the marine environment. Through the simulation of these oil fate and transport processes, with specific model limitations, oil spill models can be used for different purposes: to support operational marine oil spill response, in planning and preparing for oil spill response operations, in environmental hazard analysis of (potential) marine spills and in impact assessment on humans and the wildlife after a spill. Oil spill models range from open-source code (e.g., MED-SLIK-II, OpenOil) to commercial software (e.g., OILMAP, OSCAR) and differ in capabilities from predicting surface transport by currents, winds, and oil drift to including sophisticated weathering algorithms. In the present work, eighteen widely used oil spill models are compared. The examined models are categorized depending on the different purposes and use cases that they support and are critically reviewed regarding the input data re-

A Comparative Assessment of State-of-the-Art Oil Spill Models
Accurate prediction of transport and fate of spilled oil is challenging due to the difficulties in understanding the oil behavior and its interaction with the marine environment. Through the simulation of these oil fate and transport processes, with specific model limitations, oil spill models can be used for different purposes: to support operational marine oil spill response, in planning and preparing for oil spill response operations, in environmental hazard analysis of (potential) marine spills and in impact assessment on humans and the wildlife after a spill. Oil spill models range from open-source code (e.g., MEDSLIK-II, OpenOil) to commercial software (e.g., OILMAP, OSCAR) and differ in capabilities from predicting surface transport by currents, winds, and oil drift to including sophisticated weathering algorithms. In the present work, eighteen widely used oil spill models are compared. The examined models are categorized depending on the different purposes and use cases that they support and are critically reviewed regarding the input data required for the different processes considered, forcing models and data used, results provided, parameterization of processes incorporated, requirements in terms of the processes that need to be accurately modeled for different purposes. Deep sea blowout models are separately reviewed as the processes, the driving models and the data required differentiate, whether these are coupled or not with surface spill models.

Operational Response Models
Operational oil spill modeling has the goal of providing support to response authorities in case of an oil spill, by forecasting the transport and fate of the spill, within a short time period (a few hours) of notification of the spill occurrence. Due to the limited data and information typically available in the initial period of a spill and the large number of uncertainties introduced in oil spill models, forecasts are provided typically with a time span of 2-3 days, updated regularly as more data on spill conditions and updated ocean circulation and wind forecasts become available. For operational response, the capabilities of an oil spill model or modeling system to provide an evaluation of the uncertainty of model's forecasts are also very important. After an oil spill occurs, oil is quickly transported in the marine environment by currents, winds, and the action of waves, often in long distances. Therefore, the most important aspect in operational oil spill modeling is to capture transport, as accurately as possible, in order to direct response efforts in a timely way and to the right place. Oil weathering is also important, especially evaporation and emulsification, as oil properties are altered and different response measures are required.
All oil spill models require metocean fields as forcing, however for operational oil spill response, forecasting and monitoring systems need to be set in place, properly calibrated and validated to reduce uncertainty, and be operational in order to quickly respond to oil spill emergencies. Lagrangian oil spill fate and transport models used for operational purposes, such as GNOME, MOTHY, POSEIDON-OSM, MEDSLIK, MEDSLIK-II, OpenOil, OSCAR, OILMAP, SIMAP, MOHID, OILTRANS, OSERIT, OILTOX, and Delft3D-Part, have been implemented operationally and have been validated in several real oil spill cases. Some of these models are closely coupled to metocean forecasting systems. MOHTY, MEDSLIK, MEDSLIK-II, and POSEIDON-OSM, for example, are operational in the Mediterranean Sea region. MOTHY and MEDLSIK-II are coupled to CMEMS Med MFC models output and ECMWF wind forecasts. For POSEIDON-OSM driver models are the meteorological, ocean circulation and wave models of POSEIDON forecasting system in the Aegean Sea, while MEDSLIK is closely coupled to the CYCOFOS forecasting system for the Eastern Mediterranean, and is also coupled with CMEMS data and used in the Black Sea, Sea of Azov, Baltic Sea, Red Sea, and the Arabian Gulf. OILMAP and OSCAR are, among other regions, operational in the UK by Oil Spill Response Limited (OSR), coupled with Copernicus oceanographic fields and NOAA GFS wind fields. OpenOil is applied by MET Norway forced by CMEMS for hydrodynamics and ocean state and NOAA's GFS for wind fields, while high-resolution local forecasting models are also easily coupled and implemented. OSERIT is operational in North Sea by the UK Met Office, acquiring meteorological conditions from the global model of the UK Met Office, hydrodynamic conditions by MUMM's operational hydrodynamic models OPTOS-NOS and OPTOS-BCZ, and sea state by MUMM's operational version of WAM model. OILTOX is configured for the Black Sea, while Delft3D-part, directly coupled to Delft3D modeling suite, is used by Spill Response Group Holland (SRGH). Although virtually any of these models can be forced by any available meteorological and oceanographic model in the area of operation, after certain modifications in input format and requirements, some are more easily configurable, like GNOME, which can be integrated with any ocean circulation and meteorological model providing forecasts in different file formats, via the GNOME Operational Oceanographic Data Server (GOODS). Open source-code models are only GNOME, MEDSLIK-II, OILTRANS, MOHID, and OPENOIL, providing freely available software for users. For some cases, wave forecasts, additional to oceanographic fields, may be needed, to derive wind information from wind-wave relations. In MOHID, MOTHY, and MEDSLIK-II, wave forecasts are not currently integrated, while in MEDSLIK the forecasting SWH (Significant Wave Height) and wave period are used. Stokes drift, expressing the variation among the average Lagrangian velocity of fluid parcels and the average Eulerian velocity of fluid at fixed positions is included in OSCAR, MOTHY, POSEIDON-OSM, OILTRANS, MOHID, OSERIT, MEDSLIK-II, and OPENOIL, provided either from coupled wave models, like in POSEIDON-OSM, or, like in MEDSLIK-II, by calculating the wave spectrum using the Joint North Sea Wave Project (JONSWAP) spectrum parameterization taking the wind and fetch into account. All operational oil spill models considered here are 3D models and their bathymetry usually come from the driver ocean circulation. This coarse model provides the 3D fields of currents and temperature for interpolating the velocity field required for the Lagrangian particles' advective and diffusive displacements and for accounting for sedimentation (if included).
Concerning oil spills back-tracking that simulates scenarios in reverse to predict where an oil spill may have originated from, a limited number of models support this feature: OILMAP, GNOME, SIMAP, MOTHY, MEDSLIK, MOHID, MEDSLIK-II, OSERIT, OPENOIL, and BLOSOM. These models appear ideal to investigate a spill of unknown origin. Furthermore, GNOME, MEDSLIK, MEDSLIK-II, OSERIT, OILTRANS, MOHID, and OPENOIL include a built-in oil database, with several oil types.
As for oil transport and weathering processes, these operational response models include all the basic processes such as advection, spreading, evaporation, dispersionentrainment, diffusion, evaporation, and emulsification. MEDSLIK, MEDSLIK-II, and OPENOIL follow the Mackay's modeling approach for emulsification parameterization. GNOME has very limited evaporation and oil-shoreline interaction capacity [181]. A relatively sophisticated natural dispersion algorithm is included in OPENOIL, principally computed via vertical mixing approach in combination with wave-breaking entrainment and the option of different droplet size-distributions. Beaching or oil-shoreline interaction is a vital process in oil spill modeling, and it is included in all operational response models, considering the significant biological, social, and economic impacts on coastal areas. On the other hand, dissolution is supported by a limited number of models, namely, SIMAP, GNOME, OSCAR, OILTRANS, MEDSLIK-II, and MOHID. However, since dissolution is mainly significant for deep sea blowout oil spill modeling, as well as estimating the toxicity effects of oil dissolved in the water column, not including the parameterization of this process in operational oil spill models is not considered a drawback. Sedimentation (OSA) algorithms are included in OILMAP, SIMAP, OSCAR, MOTHY, OILTOX, POSEIDON-OSM, MEDSLIK, MEDSLIK-II, MOHID, and Delft3D-Part. Biodegradation, which is a longterm process with specific modeling limitations can only be simulated by SIMAP, OSCAR, GNOME, and MEDSLIK-II. This process, however, is more important for environmental impact analysis and damage assessment, as well as for deep sea blow out models, in order to apply bioremediation strategies. For such cases, biodegradation algorithms, taking into account droplet size distribution and oxygen and nutrients concentration, in addition to oil composition, should be considered. Photo-oxidation, a long-term process, with knowledge gaps and not well-defined parameterization, is only included in Delft3d-Part model. The vertical movement of oil droplets and the vertical mixing of oil into the water column is supported by most oil spill models (OILMAP, SIMAP, OSCAR, MOTHY, POSEIDON, OILTOX, MOHID, OILTRANS, OSERIT, MEDSLIK, MEDSLIK-II, OPENOIL, and Delft3D-Part). Ultimately, the buoyancy effects, in which large droplets rise faster than small droplets, are supported by OPENOIL, MEDSLIK, MEDSLIK-II and OSERIT.
The evaluation of uncertainty of oil spill forecasts is an important feature to be included in an operational oil spill model. Such stochastic components are included in GNOME, SIMAP, OILMAP, OSCAR, MOTHY, OILTOX, MEDSLIK, MEDSLIK-II, OIL-TRANS, OSERIT, MOHID, OPENOIL, and Delft3D-part. In GNOME, for example, uncertainty algorithms regarding the perturbation of current and wind fields are included, providing an uncertainty evaluation for the model's best estimate. In MEDSLIK-II, an uncertainty module has been included to automatically estimate prediction uncertainties related to the initial conditions of the spill. Uncertainty in oil spill forecasts can also be reduced by performing ensemble oil spill simulations using different meteorological and oceanographic forcing data available for an area, typically by different external providers. Such ensemble oil spill forecasts were a basic feature of the Mediterranean Decision Support System for Maritime Safety (MEDESS4MS) developed in 2010 (Zodiatis et al. [225]) where MEDSLIK, MEDSLIK II, MOTHY and POSEIDON OSM are used. Finally, using oil spill drifting observations, e.g., from drifters and satellites, and assimilating them into operational oil spill models is essential in reducing the uncertainty of operational forecasts.

Deep Sea Blowout/Buoyant Plume Models
Blowout/buoyant plume models (e.g., CDOG, OILMAPDEEP, TAMOC) simulate oil spills originating from the sea floor or at various depths below surface and depend on complex physiochemical processes. Models which can be used operationally for subsea releases, without being too computationally expensive, are integral plume and Lagrangian particle tracking models, like CDOG, DeepBlow, which is SINTEF's multiphase integral plume model based on the Lagrangian concept and OILMAPDEEP, developed for evaluating the nearfield dynamics of a blowout plume. TAMOC provides comprehensive integral plume models (Stratified Plume Models or SPM and a Bent Plume Model or BPM) and Single Bubble Model (SBM), which tracks the fate of a single bubble or droplet as it rises through the water column, advected by the three-dimensional ambient currents, and undergoing dissolution and heat transfer. Computational fluid dynamics (CFD) models of 3D multiphase plume flow field are recent developments, which could be in the future incorporated in operational models. The MEDSLIK plume module [144] incorporated a modified model originally proposed by Yapa and Zheng [185], Zheng and Yapa [264], and Malačič [265] for entrainment of sea water into the oil plume and a revised model of Yapa and Zheng for the entrainment of oil into the water body. The mass and momentum equations in the MEDSLIK plume module differ from those provided by Yapa and Zheng [185], and no gas-oil mixture is considered. MEDSLIK has been used for risk assessment from hypothetical blowouts of existing and planned offshore platforms in the Eastern Mediterranean Levantine basin [266,267].
Near field models are usually coupled with far field models for the seamless simulation of both environments. OILMAPDEEP is coupled with OILMAP or SIMAP, DEEPBLOW with OSCAR, while TAMOC is now integrated with the GNOME modeling suite. BLOSOM is an integrated modeling system considering offshore oil spills in deep water until they reach the surface and including also surface transport and weathering processes. All these integrated models provide fate/weathering and transport process except from CDOG, which is an integral plume model. Moreover, the only open-source model in this category is GNOME. A back-tracking component is considered in GNOME and BLOSOM and a stochastic component is included in all models apart from CDOG. CDOG is the only model running in the 2D domain, without bathymetry input requirement. An integrated oil database is included in GNOME and TAMOC. Ocean currents, density stratification, biodegradation of oil droplets (included in GNOME, SIMAP, OSCAR, BLOSOM and TAMOC), oil dissolution (considered in all models reviewed here except for CDOG and BLOSOM), and the oil droplet size distribution from a deep sea blowout greatly influence subsea oil fate and trajectories. In the Deepwater Horizon accident, the significant role of biodegradation of dissolved oil and oil droplets was indicated in several studies (e.g., [85]). Finally, blowout models face challenges, due to the fact that there are gaps in measurements of deep sea currents. These currents are weaker, but have crucial impact on the oil transport.

Spill Response Planning and Environmental Impact Assessment
As described above, all analyzed oil spill models are deterministic trajectory models that predict the movement of an oil spill based on a specific scenario and varying over time under the prevailing metocean conditions. These models may be used in both response and contingency planning scenarios. Except from CDOG, OILMAPDEEP, and POSEIDON-OSM, all models include a stochastic component that demonstrates the probability that an oil spill may impact a specific area within predefined time periods. This is done by performing a series of model runs to produce multiple trajectories under various met ocean conditions based on historic records. These outputs illustrate the marine regions and shorelines that are mostly at risk from oiling during various seasons. Oil-shoreline interaction or beaching processes are essential for response planning and environmental impact assessment, considering the significant ecological, and socio-economic impacts of oil spills on coastal areas. These processes are included in all eighteen models, apart from CDOG, which is a plume model, and DELFT3D-PART. On the other hand, dissolution increases the toxicity of water, causing ecological impacts on marine life. Therefore, dissolution is important for environmental impact analysis, risk and hazard assessment, and response planning modeling. Emulsification, sedimentation, and biodegradation are processes with specific model limitations and play a significant role on risk assessment, biological impact analysis, and response planning modeling. Models with response support are GNOME, OILMAPDEEP, SIMAP, OSIS, OSCAR, OSERIT, OILMAP, MOTHY, MED-SLIIK, MEDSLIK-II, OILTRANS, MOHID, OPENOIL, and BLOSOM. As for environmental impact analysis, which identifies the spill impact on marine species, such as fish, plankton, macroinvertebrates, and macrophytes, this feature is included in SIMAP, MOHID, and OSCAR. Moreover, the models which are mainly implemented for research purposes are CDOG, OSCAR, MOHID, POSEIDON-OSM, MEDSLIK-II, and DELFT3D-PART. Although the main purpose of CDOG is research, it has been implemented for response purposes and the USA government agencies (MMS, NOAA) and oil companies have recently started using the CDOG model. Models with risk and hazard assessment are SIMAP, OSCAR, OSERIT, and DELFT3D-PART.

Models Performance against Field Data
In the Deepwater Horizon (DWH) spill, extensive validation of oil spill models have been performed by Spaulding et al. [190,191] for the nearfield subsurface blowout and by French Mc-Cay et al. [192] for the far field transport. These validations used state-of-the-art models (OILMAPDEEP for the nearfield and SIMAP for the far field) and were based on observations of the oil spill openly accessible to the public. In parallel, OILMAP has been validated on a global scale [218] and it has been implemented in Persian Gulf War spill and Braer spills [122], in Dubai and Gulf region [219], in Southeast Monsoon and in Indian Ocean [268], in Mersin Bay in Mediterranean Sea [269], and in Bay of Samsun in Black Sea [270]. Moreover, SIMAP has been validated against data of more than 20 large spills, such as the Exxon Valdez [37,38,208] and [117]. In parallel, the simulation of SIMAP model in the pack ice provides results same as the experimental observations [271].
Socolofsky et al. [241], by means of Deepwater Horizon oil spill study, presents a comparison of OSCAR, BLOSOM and OILMAPDEEP. Moreover, this study identifies that OSCAR computes a lower surface flux when biodegradation and dissolution are included. In parallel, only OSCAR simulates the plume trajectory into a few oscillations of the intrusion layer. The majority of these models predict the large and small oil droplets to enter the far field together, either at the first neutral buoyancy level or the maximum rise height. In OSCAR, the oil reaches the surface almost immediately via the plume, thus OSCAR generally predicts surfacing to occur furthest downstream. In addition, OSCAR has been widely applied in oil spill risk evaluation and for response planning and operations [121] and for numerous hindcasts and predictions of oil spill accidents in the Northern and Baltic Sea, Gulf of Mexico, Mediterranean Sea [272,273], and in the Caspian Sea [274]. Further effort to upgrade the OSCAR model has been provided by Nordam et al. [275]. The comparison of OSCAR with the Prestige's first and second oil slick evolution with drifter buoys trajectories have demonstrated the capacity of the model to identify the areas with the strong possibility of being influenced by a specific oil spill and the arrival time of the oil spilled at a specific coastal region. Thus, further study is needed to incorporate these improvements into the OSCAR model, including a hindcast of ocean currents for the case study, increasing the spatial resolution in the numerical field near the shore and assigning a different probability of occurrence for launch points considering the most frequent transport routes. In addition, there are many hypothetical functions in the dispersant algorithm and subsequently more empirical data are needed to overcome this uncertainty. BLOSOM predicts oil surfacing over an order of magnitude farther downstream since it simulates the oil rise passively (i.e., outside the plume) over the depth of the water column. It highlights the importance of the rising plume in shallow water and near the blowout source. The jet plume element of BLOSOM has been evaluated by means of experimental cases, which took place in the North Sea [242,243]. Moreover, the ASA OILMAP DEEP model predicts a farther downstream surfacing zone due to the assumed smaller droplets in size distribution. Lower oil surfacing fluxes either result from dissolution and biodegradation (e.g., OSCAR) or from a large quantity of small droplets that have not yet surfaced at the end of the simulation (e.g., OILMAPDEEP). Furthermore, TAMOC has been validated through various experimental works of bubble plumes, such as [72].
Socolofsky et al. [72], present a review of biodegradation algorithms and implementations in oil spill models (OSCAR, SIMAP, TAMOC, BLOSOM, GNOME). Moreover, this study applies the TAMOC model to evaluate how differences in biodegradation formulations affect the predicted fate of oil droplet size distribution on subsurface transport. This study highlights that the approach of pseudo-components for various oils is not the same in SIMAP and the other oil spill models. For dissolved hydrocarbons, OSCAR, SIMAP, BLOSOM, and GNOME parameterize biodegradation as a first-order decay process, with various rates for different pseudo-components. For undissolved oil slicks and droplets, these models also use a simple first-order decay. First-order decay is usually considered not only for computational efficiency but is also to some extent justified by the limited solubility of most hydrocarbons. Effects of non-linear kinetics (e.g., Monod kinetics) become important at hydrocarbon concentrations comparable to the half-saturation constant. Socolofsky et al. [72] also considered the effect of the oil-water interface (e.g., droplets size distribution) on biodegradation and employed the TAMOC model, which includes such a surface-area-dependent degradation algorithm, to find that droplet diameter had the most significant effect on subsurface oil biodegradation and transport. According to the Scotian Shelf case study, the trajectory results of the MOHID and OSCAR models are remarkably similar. Differences in oil mass balance among MOHID and OSCAR indicates the following key challenges in the weathering process of MOHID: (1) the evaporation algorithm is purely empirical; (2) the dispersion algorithm defined via Delvigne and Sweeney [110] is purely empirical and may not well reflect the dispersion rate depending on oil viscosity; (3) the lack of a biodegradation algorithm creates uncertainties in the oil mass balance for the long-term simulation, and may not simulate a reduction in the dispersed oil.
Duran et al. [240] have implemented and compared GNOME and BLOSOM for the hindcast of the 2003 Point Wells oil spill. BLOSOM provides more accurate results due to the inclusion of an internally calculated deformation angle for the formation of wind energy. GNOME uses a random walk with a constant diffusion coefficient in order to simulate spreading of oil due to turbulent diffusion, while BLOSOM considering both the diffusion scheme and the diffusion coefficient provides another option as described in Duran [240]. A random walk with constant diffusion coefficient was also used for BLOSOM, in order to provide the predictions absolutely comparable, affecting the way oil beaches along the coast. Moreover, BLOSOM provides an equidistant spatial reference system with units of meters, while GNOME computes oil particle trajectories via longitude and latitude directly, converting these computations compatible with the metric system. These different formulations lead to some deviation in orbits that can be estimated by calculating an orbit, while maintaining a constant ocean velocity in space and time. Both models incorporate ocean currents using a Euler integration approach. In parallel, GNOME includes the process of beaching of oil in coastline, while BLOSOM does not include a refloat option [240].
Moreover, OILTOX was implemented in the North-Western Shelf of Black Sea in September 2002 and in Dnipro-Boog Estuary in 2002 [227], so limited test cases have been implemented in order to highlight the accuracy and reliability of the model. Additionally, the OILTRANS model has been validated through an accidental release during a shipto-ship fuel transfer in the Celtic Sea in February 2009 [203]. Comparisons with aerial observations of the accidental oil spill, and subsequent model simulations, indicate that the OILTRANS model has the capability of accurately predicting the transport and fate of the oil spill [203]. Furthermore, validation of the oil fates component of OILTRANS was made against the ADIOS oil model. The validation exercise verified that OILTRANS algorithms show a good agreement when compared against results from the ADIOS oil model application [203]. The CDOG model, Zheng et al. [186], has been compared with large scale and unique field experiments [170] performed in Norway [172]. After that, CDOG model was revised and improved by other scholars, and the underwater oil spill was simulated and predicted successfully [156,234,276]. Recently, the USA government agencies (MMS, NOAA) and oil companies have started using the CDOG model.
Some limited test cases of the DELFT-PART model are included in [245,246]. These studies predict the surface oil trajectories of the Penglai 19-3 oil spill. Comparing the oil-spill simulation results with remote sensing information, that the model with variable temperature and salinity exhibits a better fit with the actual observations. According to the result of the variable thermohaline model, density currents should be taken into account in future oil-spill studies to create more realistic outputs [246].
In parallel, OSERIT model has been validated against various academic and real case studies, including the Gannet platform accident in August 2011 [204]. The model outputs have been compared against satellite images of the oil spill taken from radar sensors. It showed that OSERIT simulates the position of the observed oil slick accurately.
Furthermore, the OpenOil model has recently indicated a great agreement with satellite observations of the DWH oil slick [261,262] for two seven-day simulations forced via the metocean results from the high resolution GoM-HYCOM 1/50 (in the GoM), FKEYS-HYCOM 1/100, and ECMWF models [262]. The performance of the integrated oil spill model is assessed through the comparison of model simulations with airborne observations of an oil slick. The outputs demonstrate that an accurate description of fate and transport processes, in particular vertical mixing and oil weathering, is needed to represent the horizontal spreading of the oil spill [112]. This study found that the difference between two different oil droplet size distribution is negligible, demonstrating that the effect is strong in terms of the configuration option [262]. On the other hand, according to Röhrs et al. [112] and the study of OPENOIL model, small dispersed droplets should not be ignored in oil spill modeling as they can appear as resurfaced oil. Furthermore, this study highlights that initial weather conditions at the time of release affect long-term transport by determining whether the oil is first emulsified and thus retained on the surface or first submerged and therefore protected from further emulsification [112].
During MEDESS-4MS project, four comprehensive stand-alone oil spill models (MED-SLIK, MEDSLIK-II, POSEIDON-OSM, MOTHY) implemented in the Mediterranean Sea were interconnected into an integrated multi-model oil spill forecasting network [225]. In order to accomplish this same integrative step, each stand-alone system should be available to the oil spill systems, environmental data from the CMEMS, the national ocean forecasting systems, and the oil spill information via monitoring platforms (REMPEC, EMSA CSN). POSEIDON-OSM and MEDSLIK both used the same ocean and hydrodynamic forcing data: the low resolution ocean forcing from the Copernicus Med-MFC product with 6.5 km horizontal resolution, the atmospheric forcing from ECMWF winds. Results have shown that the values of the fate processes differ between the two models: POSEIDON-OSM marks for the evaporation and emulsification are lower than the MEDSLIK model (25.6% and 66.7% compared to 35.8% and 70.8%, respectively) and the percentage of the oil remained in the sea is 74.4% in POSEIDON-OSM and 63.5% in MEDSLIK. These results stem from the different parametrization of these oil spill models (evaporation, emulsification and hydrodynamic processes). For this reason, further study is needed along with in-situ data in order to produce more accurate results.
MEDSLIK-II has been validated with in situ data, surface drifters data, and with satellite data [89,123,172]. Interesting applications of the model are included in [176,177,179,[253][254][255][256][257]. MEDSLIK-II is used operationally in the Mediterranean region, allowing also support to REMPEC for oil spill emergencies in the entire Mediterranean basin. MEDSLIK-II has been applied to forecast oil spill fate and transport during numerous emergency cases in the Mediterranean Sea (the MS Und Adriyatik oil spill in 2008 and the Costa Concordia oil spill in 2012), while MEDSLIK has been applied operationally during the biggest so far oil pollution in the Eastern Mediterranean, i.e., the Lebanese oil pollution crisis in the summer of 2006, following a request from the European Civil Protection and REMPEC). MOTHY has been successfully compared with observations of Erika pollution incident [198]. Table 2 represents eighteen oil spill models according to their general features and the transport and oil weathering processes that are included.  Table 3 presents five of the most widely-used operational oil spill models (MOTHY,  POSEIDON-OSM, MEDSLIK, MEDSLIK-II, and OPENOIL), according to their input data requirements and their characteristics.

Conclusions and Recommendations
Over the latest decades, the actual number of oil spills has been decreasing. However, pollution from oil spills still remains a crucial issue for assessing the environmental footprint of the oil and gas industry. Clearly, nature conservation demands efforts to safeguard that oil spill accidents will be minimized, as far as possible. Thus, government and industry have to collaborate in order to decrease the risk of oil spills via the insertion of strict new legislation and strict operating codes [54]. Industries have enforced new operating and maintenance measures to diminish accidents resulting in spillages. In parallel, comprehensive training programs have been developed in order to decrease the possibility for human error. Oil spill modeling involves several parameters and includes processes that require further study and intensive research [54,184], especially in the field of improving the parameterization of OWPs and the capacity to operationalize model execution in real-time.
It is apparent from the present work that the simulation of oil spills and the parameterization of oil fate and transport processes vary considerably among these eighteen oil spill models reviewed, but almost all of them are critically dependent on metocean forcing. Currents, waves, and winds are the most crucial elements; to accommodate this need most state-of-the-art oil drift models have the ability to import data from external databases (e.g., NOAA and CMEMS) for currents, waves, wave-induced drift, air temperature, water temperature, salinity, and turbulent kinetic energy. Nevertheless, the method with which these external data are used varies according to the formulation indicated by the specific oil spill model [145]. Furthermore, the movement of oil across the water column and on the sea surface is a process with specific model limitations, as chemical and biological fluctuations in the oil influence its physical properties, which subsequently impact the fate and transport of oil spill in the marine environment.
Oil spill models advanced significantly over the recent decades, with major improvements in trajectory forecasting and oil weathering processes [10][11][12]14,277]. However, few comparisons [72,225,241] of oil spill models have taken place with actual oil spills to assess their performance. More broadly the models are compared against drifter trajectories [177]. Overall, modern oil spill models follow the Lagrangian approach and incorporate the processes of spreading, advection, and diffusion of oil, together with a standard set of transformation (fate or weathering) processes [145]. The majority of existing Lagrangian oil spill models are surface models. These models use the same or similar semi-empirical relationships obtained from laboratory and field experiments [14]. The most crucial oil weathering processes are spreading, evaporation, emulsification, dissolution, biodegradation, and photo-oxidation. Evaporation, emulsification, and dissolution obtain a well-defined parameterization and are formerly included in operational oil spill models. Nevertheless, the separation between these three fate/weathering processes is lacking in operational response models. In addition, biodegradation, depending on oil droplet size distribution, plays a vital role in biological impact analysis and in long-term response planning modeling. Nevertheless, as a general rule, processes such as dissolution, vertical mixing, photo-oxidation, and biodegradation are lacking in the majority of existing oil spill models. Therefore, users of the new generation of oil spill models demand not only oil spill predictions but also an assessment of ambiguity of such forecasts, which is crucial and urgent for up-to-date, beneficial, and cost-effective responses. Due to this lack in numerical oil spill modeling, there is still a need for further study and enhancement in existing oil spill models. In particular, the development of a biodegradation algorithm by means of optimum parameterization via real-time laboratory results and the incorporation of this algorithm in an advanced existing oil spill model, providing an efficient and cost-effective response on oil spills, is required.
Moreover, as mentioned above, there are several wide-spread and comprehensive Lagrangian oil spill models, such as SIMAP, OSCAR, OILLMAP, MOHID, MEDSLIK, MEDSLIK-II, and OPENOIL. Two open-source code Lagrangian oil spill models (OPENOIL and MEDSLIK-II) are widely used by the scientific community. Firstly, the all-inclusive oil drift model OPENOIL for the transport and weathering of oil spill in the marine environment considered significant mechanisms for vertical mixing of oil and includes the effect of different parameterizations for the size distribution of dispersed oil droplets. On the other hand, the horizontal transport of oil spills relies on oil type, meteorological data, and the levels of turbulence. As for MEDSLIK-II, it is also a freely available open source Lagrangian particle tracking model and has been extensively used operationally in many oil spill accidents in recent years.
In parallel, complications in comprehensive oil spill modeling maintain when the time of oil release is unidentified or when the oil spill model is formed from observed oil slicks. Regularly, the size of oil spill is not identified, thus the process of emulsification and the whole amount of oil has to be predicted. In order to minimize the model sensitivity derived from uncertainties in the emulsification rates and slick thickness, the use of the droplet size distribution formulation [127] could be proven helpful in operational models. This approach is less sensitive to the emulsification rate and it does require knowledge of oil slick thickness [12,112]. Mixing and surface interaction approaches utilize the vertical variation of eddy diffusivity as prevalent by means of the mixed layer [112]. Ocean circulation models with high resolution and data assimilation schemes contain details on the vertical profile of currents, stratification, and turbulent mixing, allowing more accurate real-time particle transport predictions [139].
Taking everything into account, further improvements in oil spill modeling should focus on (a) the drift and evolution of oil spills, (b) the inclusion of small dispersed droplets that can appear as resurfaced oil, (c) the comprehensive parameterization of dissolution, a process very important for blowout plume and subsea oil spill modeling, (d) the development and incorporation of a biodegradation algorithm in an operational oil spill model, with optimum parameterization via real-time laboratory results, and (e) further potential advances with more detailed approaches for the vertical mixing of oil particles and improved resolution and downscaling of metocean models. Finally, more studies should be carried out to determine whether oil is firstly emulsified and remains at the sea surface or firstly submerged and protected from further emulsification [112].
Further improvements are required on: 1. the integration of blowout and droplet size distribution model, since the majority of operational oil spill modes are surface. The merging of blowout/ buoyant plume model with the surface oil spill model is essential, considering an oil droplet size distribution algorithm and integrating it in operational models; 2.
the improvement in the parameterization of oil transport, since the accuracy of the transport process depends on the accuracy of the circulation and atmospheric models. Furthermore, a droplet size distribution algorithm should be taken into account for transport processes in the future operational response models; 3.
the improvement in the parameterization of entrainment, since evaporation and entrainment are highly complex processes, and generally inconsistently handled in spill models. Entrainment is a function of white capping/wave breaking, which is an intermittent process. Droplets are formed, then are entrained, they rise at various rates to the surface. This is not well handled by assuming wave-averaged Eulerian velocities or average dissipation rates. In addition, treating the oil as a film for the evaporation and as droplets for other processes is inconsistent. For this reason, future operational model should include the wave spectrum, the white capping and entrainment of oil, depending on oil droplet distribution; 4.
the parametrization of photo-oxidation and integration in operational oil spill models, since the existing operational models do not include photo-oxidation, as there is limited knowledge about the process and parametric expressions are missing. The importance of this process has not yet been properly understood, thus an algorithm should be developed; 5.
the parametrization of MOSSFA process and its integration in the operational oil spill models, since the current operational models do not include the MOSSFA process, which has significant role in risk and hazard assessment, in biological impact analysis and in response planning modeling. Informed Consent Statement: Not applicable.

Data Availability Statement:
This study is a review paper, and it does not report any data.