Spreading refers to the creation of a thin film, expanding over the sea surface, as soon as oil is being released ^[1][2][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 ^[3][18]. The most widely-used spreading algorithms have been developed by Fay ^[4][5][19,20] and Hoult ^[6][21]. The theory of gravitational spreading against viscous resistance is also followed in the Mackay’s fate algorithms ^[7][8][9][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 ^[10][25], turbulent mixing and wave breaking ^[11][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 ^[12][27], as well as into natural entrainment ^[13][28]. Some recent modifications and improvements in spill spreading estimation appear in the literature (e.g., ^[10][12][13][25,27,28]). The studies of Korinenko and Malinovsky ^[14][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. ^[11][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 ^[11][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 ^[15][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 ^[10][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 ^[1][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 ^[10][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 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 ^[16][17][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 ^[16][18][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 ^[19][31]. On the contrary, the viscosity of remaining “weathered patches” increases ^[17][30], leading to severe physical and chemical effects on the marine environment ^[16][10]. The most widely-used analytical method to assess the rate of evaporation is based on the work of Stiver and Mackay ^[20][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 ^[21][22][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 ^[10][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 ^[23][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 ^[10][25], while similar approaches of pseudo-component evaporation models are used in other models like in OSCAR ^[24][25][36,37] and SIMAP ^[26][27][38,39]. Contrary to the widely-used Mackay well-mixed boundary layer evaporation model, Fingas ^[28][29][30][40,41,42] suggested a different approach for oil evaporation modeling. Fingas ^[29][30][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 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” ^{[16][18][31][32][33][34]}[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 ^{[35][36][37][38]}[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 ^[17][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 ^[7][8][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 ^{[39][40][41][42]}[50,51,52,53]. Fingas ^[43][44][54,55] and Fingas and Fieldhouse ^[45][46][5,56] introduced a stability index (SI) according to density, viscosity, and type of oil in order to classify the emulsification tendency of oil ^[47][12]. In addition, SINTEF’s (Selskapet for INdustriell og TEknisk Forskning) data-based oil weathering model (OWM) ^[48][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.

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 ^[49][50][13,58]. Dissolution and evaporation are two competitive processes ^[51][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 ^[18][49][50][7,13,58]. Generally, dissolution is significant when evaporation is low ^[17][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 ^[47][12]. The algorithm developed by Mackay ^[52][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 ^[52][60]. In SIMAP ^[26][27][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 occurs when oil under the influence of the sunlight generates polar, water soluble, oxygenated compounds ^[35][11]. The process depends on the type of oil ^[53][1] and on the thickness of the oil slick ^[30][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 ^[53][1]. Thus, photo-oxidation is considered unimportant over the first few days of a spill but becomes visible after a week or more ^[16][52][54][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 ^[54][55][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 ^[56][57][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 ^[58][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 ^[35][11].

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 ^{[59][60][61][62][63]}[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 ^[59][63][64][66,70,71]. Furthermore, the rate of oil biodegradation ^[65][66][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 ^[53][1] is still debated, as studies during and after the Deepwater Horizon oil spill showed both enhancement and inhibition of biodegradation ^[67][68][74,75]. As knowledge is still limited on how dispersants affect microbial species and their ability to biodegrade oil ^[68][75], as well as on the effects of different dispersant-to-oil ratios and the chemistry of oil-dispersant mixtures ^[67][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 ^[55][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 ^[26][27][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., ^{[69][70][71][72][73][74]}[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. ^[75][76][77][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 ^[66][73]. Work performed by Kapellos ^[78][85] and Kapellos et al. ^[79][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 ^[65][72], employing the TAMOC model ^[80][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. ^[81][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 ^[82][83][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 ^[84][85][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) ^[69][86][76,93].

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) ^{[17][30][35][53][87][88][89][90]}[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 ^[90][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 ^[91][98]. In regions with high SPM concentrations, increased dispersion and removal of oil is accounted for due to ingestion and adhesion ^[17][24][30,36]. Several parameters (e.g., temperature, salinity, wave energy, and physio-chemical oil properties) may control the OSA formation ^[92][93][99,100]. Moreover, the properties and characteristics of sediments constitute a significant role in OSA formation ^[91][93][98,100]. Due to knowledge gaps in properly expressing the detailed dynamics of sedimentation in a quantitative parameterization scheme, data are limited ^[49][57][13,64]. Khelifa et al. ^[94][95][96][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 ^[15][14]. Recently, a conceptual development of oil–particle coagulation capability was developed by Zhao et al. ^[97][98][104,105]. Furthermore, the new term MOSSFA (marine oil snow sedimentation and flocculent accumulation) was identified in 2013 ^[99][106], after the Deepwater Horizon (DwH) accident, in order to assess the procedures affecting the formation and fate of oil-associated marine snow ^{[100][101][102]}[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. ^[102][109]. Finally, the MOSSFA process has not been incorporated into any existing operational oil spill model.

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 ^{[30][103][104][105][106]}[42,110,111,112,113]. The large droplets resurface to their primary region while the smaller spread and diffuse into the water column ^[103][110]. The rate of natural dispersion is influenced by environmental frameworks (i.e., the sea state), but also by oil properties and spill characteristics (oil-film thickness, density, viscosity, oil/water surface tension), developing rapidly with low-viscosity oils in the presence of breaking waves ^[15][107][14,114].

Mackay ^[52][60] and Mackay et al. ^[7][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. ^[7][22] has proven to succeed only at moderate wind speeds ^[49][108][13,115]. Delvigne and Sweeney ^[103][110] and Delvigne ^[109][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. ^[7][22] and Delvigne and Sweeney ^[7][103][22,110], as well as their modifications, are widely used in operational oil spill models, like ADIOS ^[10][21][25,33], SIMAP ^[26][110][38,117], OSCAR ^{[111][112][113][114]}[118,119,120,121], OILMAP ^[115][122], MEDSLIK ^{[116][117][118][119]}[123,124,125], and MEDSLIK-II ^[82][83][89,90]. In OpenOil ^[105][112], the method of Li et al. ^[120][126] is followed, which is a modification of the formulation of Delvigne and Sweeney ^[103][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 ^{[121][122][123]}[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 ^[10][25]. It is regularly parameterized via the wind speed (e.g., in ^[124][130], subsequently used in ^[103][110]). Currently, there are numerous formulations for the wave-breaking fraction (e.g., ^{[124][125][126]}[130,131,132]). However, vast uncertainty still exists in the linkage between the wind speed and wave breaking areal fraction ^[105][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 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 ^{[10][47][127]}[12,25,133]. Oil droplets resurfacing has been modeled by Tkalich and Chan ^[104][111] and used in OpenOil ^[105][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 ^{[105][106][128]}[112,113,134]. Droplet size distributions of dispersed oil may be declared either as a number size distribution or as a volume size distribution ^[127][133]. Although ^[103][110] noticed a power-law number size distribution, current experimental research indicates that the droplet size distribution is better expressed via a lognormal distribution ^[121][127][127,133] or as two regimes with various exponents of power-law ^[127][129][133,135]. Identified droplet size distributions depend on oil type, sea state, oil weathered state, oil-water interfacial tension, and initial oil slick thickness ^[121][127][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 ^[105][112].

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 ^{[47][100][104][127]}[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 ^[105][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 ^[130][131][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 ^[132][137] or be approximated by the wind speed (e.g., ^[133][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 ^[121][127]. Novel ocean models with real-time data present details about the vertical currents, stratification, and turbulent mixing, providing more sensible particle transport representations ^{[12][105][134]}[27,112,139]. Vertical mixing algorithms and parameterizations are provided by Galt and Overstreet ^[13][28], Röhrs et al. ^[105][112], and Nordam et al. ^[135][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 ^[136][141]. From this perspective, a random vertical displacement is estimated for each particle (e.g., ^[105][112]).

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 ^[137][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 ^[105][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 ^[108][138][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 ^[105][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 ^[139][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 ^[140][145], while examples of the sensitivity of current depth in oil spill modeling are provided by ^[83][90]. Improvements in parameterization of wind drag have been introduced by ^[141][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.