Modes of Operation and Forcing in Oil Spill Modeling: State-of-Art, Deficiencies and Challenges

1. Introduction

When crude oil is accidentally released in the marine environment, an oil slick is formed appearing as a thin oily layer, floating on the sea surface [1]. The slick is shaped by the slow, low-scale, diffusive processes, responsible for changing the contaminants' concentration, and is transported by the large-scale advective processes, advancing the center of the oil slick mass to the direction of background currents, winds, and waves [2]. This implies that the ambient environment of the spill significantly determines its movement and fate. The amount of oil spilled in the ocean, the oil's initial physico-chemical properties, the prevailing weather and sea state conditions, and other spill-specific and environmental factors affect the timing, duration and relative importance of each physical and biochemical oil weathering process (known as OWP), affecting the slick [3,4,5]. Since hydrocarbons are non-conservative pollutants, OWPs cause long-term changes in their physico-chemical properties, like oil density and viscosity [6]. The most important OWPs are spreading, evaporation, dispersion/diffusion, emulsification, and dissolution. Photo-oxidation, biodegradation and sedimentation act over longer time periods and determine oil’s ultimate behavior [7].

Planning for and responding to an oil spill requires rigorous comprehension of baseline environmental characteristics and processes [8]. Oil spill models help the response agencies lessen the damaging impacts on the environment by predicting the path of at-sea oil spills. Predicting the spillage trajectory is the main outcome of oil spill models, highlighting the potential for ecosystem harm as an incident develops, while in parallel assisting in the optimization of the cleanup efforts [9,10,11]. Any guidance that oil spill modeling may offer could be extremely important for the authorities, given the tremendous effect and costs associated with oil spills. Risk evaluation, readiness planning and analysis of the environmental effects of the oil industry infrastructure, heavily rely on oil spill modeling [12]. When models are run, a wide range of input variables and actual met-ocean conditions might result in multiple alternative trajectories [13]. Following analysis, these trajectories are plotted on maps to create reaction strategies. Emergency responders must be knowledgeable about the type of oil, the location, and the marine and coastal habitats the spillage may affect. Thus, governments, oil exploration and production firms, insurance companies, and other stakeholders may evaluate whether the adequate resources, tools, and procedures are in place to respond in oil spill incidents. Simulating different scenarios may allow to assess the potential environmental impacts and device plans on the movement of the response supplies to the necessary locations [14,15]. This procedure could lead to the assessment of the efficacy of various response strategies, as well as their benefits and drawbacks [10,12]. Additionally, it is expected to aid responders to organize and mobilize socioeconomic resources to limit environmental impacts along the oil's potential course [8].

As explained above, met-ocean conditions, i.e., currents, wind and waves, represent the fundamental components influencing the spreading of oil in the marine environment [7,16,17,18]. For this reason, it is crucial to be able to illustrate that oil spill forecasts are accurate and reliable, as well as that the constraints of a model are well-understood when evaluating model’s predicting capacity and performance [19,20]. An assessment of the ocean currents, water characteristics over the water column and waves at a particular time and location is provided by the three-dimensional ocean circulation models [21]. These models aid in determining how these factors will affect the transport of oil, once it reaches the sea surface. Meteorological and atmospheric models provide information on air properties, like temperature, relative humidity, barometric pressure, as well as on the surface winds that might transport and affect the evaporation rate of floating oil [22]. In parallel, wave models provide information about the significant wave height and Stokes drift fields, affecting wave turbulence, vertical mixing and oil dispersion within the water column [23].

Therefore, cutting-edge, quality-controlled, high-resolution meteorological, hydrodynamic, and wave models should be combined with oil spill models to more precisely predict the fate and weathering of spilled oil at sea. Furthermore, once the oil is discharged into the marine environment, the chemical and physical changes it will undergo could be predicted by the fate models [17,24]. Hence, biogeochemical sub-modules capable to describe the fields of nutrients, carbon and plankton are significant in determining the dominant physico-chemical processes, like biodegradation, sedimentation and photo-oxidation [7]. All of the aforementioned factors must be taken into account in operational oil spill modeling to respond to oil spills in a timely, effective, and economical manner.

This study presents a state-of-the-art review on oil spill modeling advancements, emphasizing on the met-ocean data used for their forcing. In order to demonstrate the significance of high-resolution and accurate forcing produced from background models in the accuracy of the operational oil spill model results, the current review will provide a critical comparison of the widely used hydrodynamic, meteorological, wave, and biogeochemical models commonly utilized in oil spill modeling. This effort will provide technical guidance as well as future directions and advancements for oil spill simulation. The paper is organized as follows: in Section 2, the widely-used oil spill models are shortly described; Section 3 presents the most widely-used circulation models; Section 4 provides a critical comparison of coupling oil spill models with ocean and meteorological models, through several case studies, while Section 5 suggests technical recommendations and modeling improvements and considerations.

2. The Oil Spill Models

In this survey nine state-of-the-art oil spill models, namely: OpenOil [23,25], MEDSLIK [18,26,27], MEDSLIK-II [28,29], SIMAP [30,31], GNOME [32,33,34], OILMAPDEEP [35,36,37], BLOSOM [38,39,40], MIKE 21/3 OS FM [41,42] and STFM [43,44] are examined in terms of their meteorological, hydrodynamic, wave and biogeochemical forcing in twenty-four oil accidental release cases, studied worldwide. An analytical description and comparison of these oil spill models in terms of their characteristics, capabilities, and simulated processes is presented in the comprehensive review of Keramea et al. [7].

OpenOil is a newly-developed, open-source oil spill transport and fate model [45], part of python-based trajectory framework of OpenDrift [25]. To reach operational oil spill forecasts with OpenOil, MET Norway employs in-house, high-resolution ocean circulation and meteorologic models [23]. However, the model allows the coupling with the coarser resolution forecasts from CMEMS (Copernicus Marine Environmental Service), FVCOM, SHYFEM, CYCOFOS, HYCOM, Norshelf for hydrodynamics and ocean state, and NOAA, ECMWF and SKIRON wind fields with netcdf and many different files format. The OpenOil has been applied in several cases worldwide, such as the Norwegian Sea [23], the Gulf of Mexico and the Cuban coast [46,47,48,49], the Thracian Sea [50] and the Caribbean Sea [51].

MEDSLIK-II [28,29] is a version of MEDSLIK oil spill model [18,26]. MEDSLIK-II employs the experimental JONSWAP wave spectrum in terms of wind speed and fetch for the Stokes drift parameterization [52], while MEDSLIK use directly the wave heigh and period to estimate the Stokes drift. MEDSLIK-II is also coupled in terms of input format with the forecasted atmospheric fields, provided by the European Center for Medium-Range Weather Forecasts (ECMWF), and the oceanographic fields provided by CMEMS (currents, temperature, salinity and density), while MEDSLIK is coupled in addition to the CMEMS with the downscaled CYCOFOS sea currents and the SKIRON winds. MEDSLIK-II has been implemented in many case-studies in recent years, such as the Northern Atlantic [53], the Northwestern Mediterranean Sea [54], the Aegean Sea [55], the offshore of southern Italy [56] and the Brazilian coast [57]. In addition, MEDSLIK has been implemented in real oil spills incidents in the Eastern Mediterranean Levantine basin [58] and in numerous test cases in the Levantine basin [59,60,61,62], in the Black Sea [63] and in the Red Sea [64].

SIMAP, the integrated oil spill impact oil system, developed by ASA [30,31] simulates the three-dimensional trajectory, fate and transport of spilled oil and fuels, as well as the biological effects and other impacts [30]. SIMAP has been validated against data from over 20 major spills, including the Exxon Valdez [30]. The analytical description of the SIMAP oil trajectory and fate model is presented in McCay [30]. Wind-driven wave drift (i.e., Stokes drift) and Ekman transport at the surface can be modeled, based on the results of Stokes drift and Ekman transport formula produced by Youssef and Spaulding [65]. Moreover, the model has the capability to couple with 3-Dimensional hydrodynamic models (HYCOM (3-4 km), POM (10 km), SABGOM (5 km)) and with wind data from NOAA and ECMWF [66]. Currently, SIMAP has been implemented in the Gulf of Mexico [66,67]. Moreover, OILMAPDEEP, the Oil Model Application Package for deep water releases [35,36,37] has been created by Applied Science Associates (ASA) to evaluate the fate and transport of subsea emissions. Moreover, it couples with several 2D and 3D hydrodynamic model flow fields [36,37]. The results of the OILMAPDEEP model could be used as input in the SIMAP model, which is a far field fate and transport model. The model has worldwide capability and contains RPS (Rural Planning Services) using ASA’s own GIS and recently has been implemented in the Gulf of Mexico by French-MacCay et al. [67].

GNOME, the general NOAA operational modeling environment, is an oil spill model that forecasts the fate and transport of pollutants, as well as the movement of oil due to winds, currents, tides and spreading [32,33]. Furthermore, this model is highly configurable and tunable to field conditions and it can be driven by a variety of data: measured point data, meteorologic and hydrodynamic models with a variety of meshes (structured, triangular) (NOAA, ECMWF, CMEMS, POM, CROCO, RTOFS, GLB-HYCOM, FVCOM, Salish Sea model). Since GNOME can integrate any ocean circulation and meteorologic model that supports forecasts in various file formats, as well as observational data, NOAA has created the GNOME Operational Oceanographic Data Server (https://gnome.orr.noaa.gov/goods), a publicly-accessible system that provides access to all available driver models and data sources. Moreover, GNOME has been applied in many regions over the latest years, like the Indonesia [68], the Gulf of Suez in Egypt [69], offshore Odisha in India [70] and the Red Sea in Egypt [71].

BLOSOM, the blowout and spill occurrence model has been generated by the National Energy Technology Laboratory (NETL) of the USA (https://edx.netl.doe.gov/offshore/blosom/) [38,40,72]. The model may be coupled to wind and current data from different models (Salish Sea model, FVCOM, NOAA, NCOM AMSEAS) with multiple flexible file formats and output types. Finally, it incorporates a number of oil types from the ADIOS oil library [73]. Recent, the BLOSOM has been applied in the Gulf of Paria in Venezuela [74]. Furthermore, the MIKE 21/3 OS FM model contains the MIKE 21 ECO Lab/Oil Spill module, which makes use of the oil particle method to develop the oil spill model [41,42]. The MIKE 21/3 OS FM model, is coupled with ECMWF and MIKE 21 or MIKE 3. Currently, this model has been implemented in the Chinese Bohai Sea [42].

Finally, the STFM (Spill, Transport, and Fate Model) is a transport and weathering model of spilled oil based on Lagrangian elements, created by the Institute of Astronomy, Geophysics and Atmospheric Sciences, University of Sao Paulo (IAG/USP) for operation in marine and environmental fields [43,44]. Moreover, STFM is a fully three-dimensional model that uses the Weather Research and Forecasting (WRF) atmospheric model and the hydrodynamic Hybrid Coordinate Ocean Model (HYCOM), feeding the oil spill module with current speed and direction, water temperature, salinity and bathymetry data. Also, it has the capability to couple with ADIOS oil database. It has been recently applied at the Brazilian coast by Zacharias et al. [44].

3. The Oil Spill Forcing Models

3.1. Wind Fields

The National Oceanic and Atmospheric Administration (NOAA) and United States Navy government website provides several datasets, that have been widely be used as model wind data inputs in the several oil spill test-cases. Firstly, at global scale, the NOAA and the National Center for Environmental Prediction (NCEP) supports the Climate Forecast System Reanalysis (CFSR) model, which was created and implemented as a worldwide, high-resolution, linked atmosphere-ocean-land surface-sea-ice system to properly predict the conditions of these coupled domains during a 32-year period (January 1979—March 2011) [75,76]. The CFSR has a spatial horizontal resolution of 0.5° (~56 km) with hourly time step (https://www.ncei.noaa.gov/products/weather-climate-models/climate-forecast-system). The CFSR has been applied for oil spill modeling studies by French-McCay et al. [66] and Meza-Padill et al. [77]. Moreover, NOAA and the US Navy provide meteorological model outputs of the Navy Operational Global Atmospheric Prediction System (NOGAPS) with horizontal resolution of 0.5^o (~56 km) and temporal resolution of 6 hours, globally (ttps://www.ncdc.noaa.gov/data-access/model-data/model-datasets/navyoperational-global-atmospheric-prediction-system). The NOGAPS has been integrated in several oil spill test cases, as in King et al. [78], Brushett et al. [79], Le Hénaff et al. [80], Vaz et al. [81] and French-McCay et al. [66]. Similarly, the NCEP (National Centers for Environmental Prediction) provides atmospheric data from the GFS (Global Forecasting System) for dozens of atmospheric and land-soil variables, including water temperature, winds, precipitation, soil moisture, and atmospheric ozone concentration [82]. NCEP-GRF produces forecasts at three spatial resolutions of 0.25°, 0.5°, and 1° (https://www.nco.ncep.noaa.gov/pmb/products/gfs/) covering the whole globe [83]. Most oil spill models have been forced with the 0.25° horizontal resolution [44,50,70,83]. The temporal resolution of GFS is 3 hours and the NCEP contains wind velocity 10 meters above sea level, for the entire planet from May 6, 2011 to the present [84].

At regional scales, the National Center for Atmospheric Research (NCAR) and the National Centers for Environmental Prediction (NCEP) maintains the Weather Research and Forecasting (WRF) Model, which is a cutting-edge mesoscale numerical weather prediction system, intended for both atmospheric research and operational forecasting. In oil spill modeling, the WRF has been implemented by Zacharias et al. [44] in a region from 20° S, 50° W to 10° N, 20° W, with 1-hour time interval and 0.15° horizontal resolution [85]. Moreover, the NOAA NCEP system (https://www.ncei.noaa.gov/products/weather-climate-models/north-american-mesoscale) provides atmospheric forecasts for the North America through the North American Mesoscale Forecast System (NAM). NAM creates different grids (or domains) of weather forecasts with varying horizontal resolutions [86]. Temperature, precipitation, light intensity, and turbulent kinetic energy are just a few of the weather elements estimated at each grid cell. NAM employs additional numerical weather models to create high-resolution predictions over fixed regions and, on occasion, to track major weather events, such as hurricanes. In oil spill modeling NAM has been applied at 12 km horizontal resolution, with 1-hour time step, by French-McCay et al. [66]. In addition, the NOAA NCEP provides the North American Regional Reanalysis (NARR) system for weather reanalysis, (http://www.emc.ncep.noaa.gov/mmb/rreanl/) having 3-hr time step and 0.3° (~ 32 km) spatial resolution over the North America [87,88]. NARR project has been applied in oil spill modeling by French-McCay et al. [66,89]. Moreover, NARR system is based on a version of the Eta Model and its 3D variational data assimilation system (EDAS) has been operational since April 2003 [90].

Furthermore, the NOAA and FNMOC (Fleet Numerical Meteorology and Oceanography Center) Regional Navy Coastal Ocean Model (NCOM) (https://www.ncei.noaa.gov/products/weather-climate-models/fnmoc-regional-navy-coastal-ocean) provides the hindcast wind data, through the dataset “AmSeas Prior”, with a spatial resolution 1/36° (~ 3 km) (https://www.ncei.noaa.gov/thredds-coastal/catalog/ncom_amseas_agg_20130405_20201216/catalog.html). The dataset covers a time period from 2013-04-05 to 2020-12-16, with 3-hour time step. The system supports only the broader region of the Gulf of Mexico and the Caribbean Sea. Using the Navy Coupled Ocean Data Assimilation (NCODA) system, this model assimilates all available satellite and in-situ observations within the domain [91]. The NCOM model has been coupled with the oil spill model BLOSOM, as in Grubesic and Nelson [74].

In addition, the short-term model results produced by the European Centre for Medium-Range Weather Forecasts [92] have been widely used as forcing in oil spill modelling. ECMWF provides reliable daily global atmospheric forcing with three-hourly winds and spatial resolution of 0.125° (approx. 27 km). Specifically, ERA5 contains wind forcing reanalysis data (https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5) and is derived from fifth-generation ECMWF atmospheric reanalysis of the global climate, which integrates multi-source measurements with numerical simulations, using an assimilation model. This dataset has been produced using the 4D-Var data assimilation scheme and model forecasts in CY41R2 of the ECMWF Integrated Forecast System (IFS) [48]. It has a high temporal-spatial resolution (1 hour - 0.25°) and a long-time span from January 1, 1979 to the present [93]. Data can be obtained by visiting https://cds.climate.copernicus.eu/. Recent model upgrades have improved the overall performance of the forecasting system throughout the medium range. Further details on model description and verification can be found in the works of Ehard et al. [94], Haiden et al. [95] and Hersbach et al. [96]. ERA5 has been used as wind boundary forcing in various oil spill scenarios, as in Zhang et al. [83], Abdallah & Chantsev [69,71], Davis Morales et al. [51] and Liu et al. [42]. In parallel, in case of simulating retrospective oil spills, the ERA-Interim dataset could be used [97], a global reanalysis data product covering the data-rich period since 1979. Originally, ERA-Interim was run from 1989, but the 10-year extension for 1979-1988 was produced in 2011, providing data every 6 hours, with a 1/8° spatial resolution [57,98] [99]. ERA-Interim has been applied in oil spill simulations in several test cases all over the world [46,47,48,49,53,54,56,57]. Moreover, Kampouris et al. [55] has used the ECMWF ensemble prediction system at ∼ 9 km and ~18 km horizontal resolution for wind forcing.

The Eta/NCEP model [100,101] has been in operational use at the Hellenic National Meteorological Service and at the University of Athens (http://forecast.uoa.gr). Moreover, the high frequency winds from the SKIRON non-hydrostatic forecasting model [102], with 5 and 10 km spatial resolution has been utilized during real oil spill pollution [58,62] and in numerous oil spill test-cases, such as in Zafirakou-Koulouris et al. [103], Ribotti et al. [104], Zodiatis et al. [60,63], Goldman et al. [105], De Dominicis et al. [19] and Sepp-Neves et al. [53], both in the Mediterranean and the Black Sea. In parallel, the HCMR (Hellenic Centre for Marine Research) provides meteorological forecasts via the POSEIDON weather forecasting system [106], also based on SKIRON/Eta model [102], which covers an area broader than the Mediterranean basin, with a horizontal resolution of ~5 km. POSEIDON has been coupled with oil spill models, as in Annika et al. [107], and Zodiatis et al. [108] . Finally, The University of Malta (UoM) provides meteorological forecasts through the MALTA Maria ETA Model (http://www.capemalta.net/maria/regional/results.html) with a horizontal resolution of ~4 km, covering the Central Mediterranean Sea and the Maltese Islands [109]. The model has temporal resolution of 3 hours, providing forecasts for only 1 day in advance. MARIA/Eta High Resolution Atmospheric Forecasting System is based on the atmospheric limited area NCEP/Eta model [101,110] and it has been applied in oil spill case studies, like in Drago et al. [109].

The Spanish Met Office, AEMET (Spanish State Meteorological Agency) (https://www.aemet.es/es/portada), produces meteorological forcing forecasts using the HIRLAM (High Resolution Limited Area Model) [111,112]. This forecasting system runs with a horizontal resolution at 1/7° (~15 km) over the whole Western Mediterranean, providing hourly data every 6 hours, for the next 72-hours [113]. HIRLAM has been coupled to the oil spill model TESEO, and has been applied in the Prestige oil spill accidental release in the Bay of Biscay [113,114]. Météo-France contributes with atmospheric data through the ARPEGE model (Action de Recherche Petite Echelle Grande Echelle) for the entire Mediterranean basin (http://www.umr-cnrm.fr/spip.php?article121&lang=en). ARPEGE is a numerical model for global general circulation. Météo-France developed it in collaboration with ECMWF (Reading, UK) for operational numerical weather forecasting [115,116,117]. The ARPEGE model has incorporated the four-dimensional variational assimilation (4D-Var). The spatial resolution of the ARPEGE model is ~10 km in the Mediterranean Sea and the temporal resolution is 3 hours [118]. Recently, the model was upgraded in its vertical grid, composing of 105 levels, with horizontal grid of ~5 km over Europe and 24 km elsewhere [119]. This fine resolution 5 km edition has not yet been applied to oil spill modeling. Oil spill models have been coupled only to the 10 km resolution version.

Figure 1. General concept of Operational oil spill modeling.

Figure 1. General concept of Operational oil spill modeling.

Table 1. Atmospheric models with the corresponding domains and horizontal resolutions used in oil spill modeling.

Wind	Provider	Geographical area	Spatial Resolution	Data Type	Reference
Poseidon	HCMR	Mediterranean	~5 km	Forecast	[106]
HIRLAM	AEMET	Western Mediterranean	~ 15 km	Forecast	[108,112,113,114]
ARPEGE	Meteo-France	Mediterranean	~10 km	Forecast	[108,118]
SKIRON	UOA	Mediterranean and Black Sea	~5 and 10 km	Forecast	[53,102,120] [19,58,61,103,104,121]
MALTA/Maria ETA model	UOM	Central Mediterranean	~4 km	Forecast	[108,109]
NAM	NOAA/NCEP	North America	12km	Forecast	[66,86]
NARR	NOAA/NCEP	North America	0.3° (32 km)	Reanalysis	[66,89,90]
NCOM AMSEAS	NOAA/FNMOC	Gulf of Mexico and Caribbean	1/36° (~ 3km)	Hindcast	[74]
WRF	NCAR/NCEP	Regional	0.15° (~16 km)	Forecast	[44]
NOGAPS	NOAA/United States Navy	Global	0.5° (~56 km)	Forecast	[66,81]
CFSR	NOAA/NCEP	Global	0.5° (~56 km)	Reanalysis	[66,77]
GFS	NOAA/NCEP	Global	0.25° (~27 km)	Forecast	[50,69,70,71,83]
ERA5	ECMWF	Global	0.25° (~27 km)	Reanalysis	[42,51]
Era-Interim	ECMWF	Global	0.125° (~12.5 km)	Reanalysis	[46,47,48,49,53,54,56,57]

Table 1. Atmospheric models with the corresponding domains and horizontal resolutions used in oil spill modeling.

Wind	Provider	Geographical area	Spatial Resolution	Data Type	Reference
Poseidon	HCMR	Mediterranean	~5 km	Forecast	[106]
HIRLAM	AEMET	Western Mediterranean	~ 15 km	Forecast	[108,112,113,114]
ARPEGE	Meteo-France	Mediterranean	~10 km	Forecast	[108,118]
SKIRON	UOA	Mediterranean and Black Sea	~5 and 10 km	Forecast	[53,102,120] [19,58,61,103,104,121]
MALTA/Maria ETA model	UOM	Central Mediterranean	~4 km	Forecast	[108,109]
NAM	NOAA/NCEP	North America	12km	Forecast	[66,86]
NARR	NOAA/NCEP	North America	0.3° (32 km)	Reanalysis	[66,89,90]
NCOM AMSEAS	NOAA/FNMOC	Gulf of Mexico and Caribbean	1/36° (~ 3km)	Hindcast	[74]
WRF	NCAR/NCEP	Regional	0.15° (~16 km)	Forecast	[44]
NOGAPS	NOAA/United States Navy	Global	0.5° (~56 km)	Forecast	[66,81]
CFSR	NOAA/NCEP	Global	0.5° (~56 km)	Reanalysis	[66,77]
GFS	NOAA/NCEP	Global	0.25° (~27 km)	Forecast	[50,69,70,71,83]
ERA5	ECMWF	Global	0.25° (~27 km)	Reanalysis	[42,51]
Era-Interim	ECMWF	Global	0.125° (~12.5 km)	Reanalysis	[46,47,48,49,53,54,56,57]

3.2. Hydrodynamics

The Copernicus Marine Environmental Monitoring Service (CMEMS) provides several hydrodynamic datasets at global and six regional seas. In the present study, only the data products relevant to oil spill modeling will be discussed. Firstly, the Global Ocean 1/12° Physics Analysis and Forecast model, provides daily and monthly-mean data for sea temperature, salinity, currents, sea level, mixed layer depth and ice parameters, from the top to the bottom of the global ocean [122]. In addition, it produces the hourly-mean surface fields for sea level height, temperature and currents. The global ocean output files have a horizontal resolution of 1/12° (~ 9 km) and a regular longitude/latitude equirectangular projection. This dataset has been widely applied in oil spill simulations, like in the studies of Devis Morales et al. [51], Sepp Neves et al. [53] and Siqueira et al. [57]. Moreover, CMEMS provides the dataset Mediterranean Sea Physics Analysis and Forecast (MEDSEA_ANALYSISFORECAST_PHY_006_013) [123] which is being produced from a coupled hydrodynamic-wave model, implemented over the entire Mediterranean Basin. It consists the physical component of the Mediterranean Forecasting System (Med-Currents), with a horizontal grid resolution of 1/24° (approximately 4 km) and it has 141 unevenly spaced vertical levels. This dataset has been implemented in oil spill boundary forcing, e.g., in Liubartseva et al. [54], Kampouris et al. [55] and Keramea et al. [50]. The hydrodynamics are provided by the Nucleus for European Modelling of the Ocean (NEMO v3.6). The model solutions are corrected by a variational data assimilation scheme (3DVAR) of temperature and salinity vertical profiles, as well as along track satellite Sea Level Anomaly observations [124]. Finally, CMEMS supports the GLO-CPL dataset (Global Ocean 1/4° physics analysis and prediction) which is a data assimilation and forecast system that provides 10 days of 3D global ocean forecasts at ~27 km spatial resolution. The system employs the Met Office Unified Model v10.6 atmosphere configuration (at 40 km resolution) that is hourly coupled to NEMO v3.4 [125] ocean configuration and the CICE v4.1 multi-thickness category sea ice model (both on the ORCA025 grid) [126]. The GLO-CPL dataset has been used as forcing input in the GNOME oil spill simulations of Abdallah & Chantsev [71] for the Red Sea.

The NOAA National Ocean Service (NOAA/NOS) Coast Survey Development Laboratory (CSDL) runs the NOS GOM Nowcast/Forecast Model (NGOM) [127], which is an application of the POM model [128] in the Gulf of Mexico. Moreover, the spatial resolution of NGOM is 5–6 km in the northeastern and central GoM, with 37 levels in the vertical (https://www.bco-dmo.org/dataset/831523). Furthermore, the forecasts are obtained every 3 hours. NGOM has been used as forcing data in oil spill simulations, as in the case of the Deepwater Horizon buoyant plume simulation in combination with the OILMAP DEEP model [66,67]. In parallel, the NOAA and FNMOC (Fleet Numerical Meteorology and Oceanography Center) provides operational ocean predictions using the Navy Coastal Ocean Model (NCOM), with horizontal resolution of 1/36° (~3km) and 40 levels in the vertical. The model is capable to produce 4-day forecasts at 3-hour time steps. French-McCay et al. [67] and Grubesic et al. [74] have applied NCOM in oil spill simulations of OIL-MAPDEEP and BLOSOM, respectively. The AMSEAS ocean prediction system assimilates all quality-controlled observations, including satellite sea surface temperature and altimetry, as well as surface and profile temperature and salinity data, using the NRL-developed Navy Coupled Ocean Data Assimilation (NCODA) system [129].

The hydrodynamic model, Hybrid Coordinate Ocean Model (HYCOM; hycom.org), uses as outer model the operational GLoBal HYCOM (GLB-HYCOM) with horizontal resolution 1/12° (approximately 9 km) and 32 vertical layers (https://www.nrl.navy.mil/) [130]. The GLB-HYCOM model has been used in oil spill simulations, as in the case of a Brazilian oil spill model implementation, using the Spill, Transport, and Fate Model (STFM) [44], and in offshore India, coupled with GNOME [70]. In the Gulf of Mexico, the Gom-HYCOM model has 1/25° horizontal resolution, vertical resolution of 20 hybrid layers, and current predictions every 3 hours [131]. Retrospective model results are included in the reanalysis dataset of GoM-HYCOM, i.e., the HYCOM-NRL Reanalysis product (GOMu0.04/expt_50.1) produced by the US Naval Research Laboratory’s (NRL). The product has 1/25°spatial resolution (~ 3.5 km) at mid-latitudes, 36 vertical layers, and contains current predictions for the Gulf of Mexico every 3 hours. This dataset can be downloaded from these links: http://tds.hycom.org/thredds/catalog/datasets/GOMu0.04/expt_50.1/data/netcdf/catalog.html, http://hycom.org/data/gomu0pt04/expt-50pt1. Similarly, the GoM-HYCOM includes the Real-time dataset, the HYCOM-NRL Real-time [130] that uses the product HYCOM + NCODA GOM 1/25° with spatial horizontal resolution of 1/25° and 36 vertical layers, producing hourly 3D outputs in netCDF format (https://www.hycom.org/data/goml0pt04/expt-31pt0). These two datasets, HYCOM-NRL Reanalysis and HYCOM-NRL Real-time, have been used as forcing inputs in several oil spill case-studies, as in French-McCay et al. [66,67]. On the other hand, the GLB-HYCOM 1/12° is used to provide boundary conditions to a regional implementation for the GoM, having higher horizontal resolution (1/50°) with 32 hybrid vertical layers (GoM-HYCOM 1/50°) for the Atlantic Ocean areas, over the Southeastern US Continental Shelf. The GoM-HYCOM model has been implemented in a near-real time mode, by the Coastal and Shelf Modeling Group at the Rosenstiel School of Marine and Atmospheric Science (RSMAS), University of Miami (https://coastalmodeling.rsmas.miami.edu/), together with the Ocean Modeling and OSSE Center (OMOC) between RSMAS and the NOAA Atlantic Oceanographic and Meteorological Laboratory (AOML). The model covers the entire Gulf of Mexico, as well as a portion of the Caribbean Sea, the Florida Straits, and a portion of the Atlantic Ocean, along Florida, Georgia, and the Bahamas Islands. Le Hénaff and Kourafalou [22] and Androulidakis et al. [132] conducted detailed descriptions of the technical characteristics of the GoM-HYCOM 1/50° simulation (parameterizations, initial, boundary, and atmospheric forcing) and extended evaluations against non-assimilated in situ and satellite observations. The model user's manual contains additional information about the HYCOM model (www.hycom.org). The GoM-HYCOM with 1/50° horizontal resolution has been used in oil spill simulations, like in the studies of Hole et al. [47], Androulidakis et al. [46], and Kourafalou et al. [49]. Moreover, the FKeys-HYCOM model, based on HYCOM, is a high-resolution forcing model for oil spill simulations, covering the Southern Florida coastal and shelf areas and the Straits of Florida, with horizontal resolution 1/100° (~ 1 km) and 26 vertical layers. Besides, it has enabled new findings in eddy variability, with Kourafalou et al. [133,134] presenting more detailed information and data-based evaluation of its simulations. FKeys-HYCOM has been integrated in the oil spill simulations of Hole et al. [48] and Androulidakis et al. [46].

The North Carolina State University (NCSU) developed the South Atlantic Bight and Gulf of Mexico (SABGOM) hydrodynamic model, based on the Regional Ocean Modeling System (ROMS). A model implementation for the GoM exists [135,136] with horizontal resolution ~5 km and 36 vertical layers. French-McCay, et al. [66] used SABGOM in their oil spill model. Moreover, SABGOM has now been replaced by the Coupled Northwest Atlantic Prediction System (CNAPS) (http://omgsrv1.meas.ncsu.edu:8080/CNAPS/), covering a larger area than SABGOM [137]. CNAPS is a three-dimensional marine environmental nowcast and forecast model, developed by OOMG (Ocean Observing and Modeling Group). The model computes the daily fields of ocean circulation, wave, and atmospheric variables. In addition, the SABGOM developed the Intra-Americas Sea Regional Ocean Modeling System (IAS ROMS) with a horizontal resolution of ~ 6 km and 30 vertical levels. Chao et al. [138] developed an IAS ROMS simulation (version "4C") for year 2010, including a 2-km nested grid within the coarser and extended IAS ROMS domain, as part of the trustees' NRDA program.

SANIFS (Southern Adriatic Northern Ionian coastal Forecasting System) is an operational coastal-ocean model, developed by the CMCC-OPA (Euro-Mediterranean Centre for Climate Change), producing short-term forecasts. The operational chain is based on a downscaling approach that begins with a large-scale system for the entire Mediterranean basin (MFS, Mediterranean Forecasting system, e.g., Oddo et al. [139]; Tonani et al. [140]) for the derivation of the open-sea fields. SANIFS is based on the finite-element three-dimensional hydrodynamic SHYFEM model, using an unstructured grid [141,142]. The horizontal resolution ranges from 3 km in open-sea to 500-50 m in coastal areas. The model configuration has been outlined to provide reliable hydrodynamics and active tracer forecasts in the mesoscale-shelf-coastal waters of South Eastern Italy (Apulia, Basilicata, and Calabria regions). The model is forced in two ways: (a) at the lateral open boundaries, using a full nesting strategy, directly imposed by the MFS (temperature, salinity, sea surface height, and currents) and the OTPS (tidal forcing) fields; and (b) at the sea surface using two alternative atmospheric forcing datasets (ECMWF-12km and COSMOME-6km) through the MFS-bulk-formulae [143,144]. SANIFS open-sea features were validated by comparing model results to observed data, such as Argo floats, CTDs, XBTs, and satellite SSTs, as well as MFS operational products. The model's large-scale oceanographic dynamics are completely consistent with MFS [145]. SANIFS model results have been imported as sea surface boundary condition in MEDSLIK-II model [56].

Moreover, the NorShelf model, developed by the Norwegian Meteorological Institute, provides forecasted ocean currents for the Norwegian Shelf Sea, based on the Regional Ocean Modeling System (ROMS), with a 4D-variational (4D-Var) DA assimilation scheme (MET Norway). To accommodate the scale of the available observations and to compromise the need to resolve high resolution eddy dynamics, while confining non-linearities that limit the 4D-Var DA capabilities, a horizontal model resolution of 2.4 km was chosen. The model is intended to be used for the forecasting of ocean circulation and hydrography beyond the coastal area, including the entire shelf sea and the dynamics of the North Atlantic current at the shelf slope [146]. Röhrs et al., [23] used the NorShelf model outputs in OpenOil simulations.

In parallel, the IRD (French Institute of Research for the Development) provides the oceanic modeling system, CROCO (Coastal and Regional Ocean COmmunity model), which is based on the ROMS-UCLA model [147] and ROMS AGRIF model [148]. CROCO is a free surface hydrostatic C-grid model with a terrain-following coordinate system and an efficient split-explicit approach for distinguishing between barotropic and baroclinic terms. It is the oceanic component of a complex coupled system that includes the atmosphere, surface waves, marine sediments, biogeochemistry, and ecosystems, among others [68,149]. CROCO also offers MATLAB-based pre-processing tools [150] for creating a model grid, interpolating atmospheric and oceanic data as boundary and forcing input, and setting tides and rivers in the grid model. MPI and OPENMP computations are supported by CROCO. In GNOME oil spill simulations, CROCO has been implemented with 1 km horizontal resolution and 25 vertical layers in the study of Nugroho et al. [68].

In addition, the DHI (Danish Hydraulic Institute) developed the MIKE 21 model (https://www.mikepoweredbydhi.com/areas-of-application), a comprehensive simulation system for hydraulics and hydraulic-related phenomena in estuaries, coastal waters, and seas [151]. It simulates unsteady two-dimensional flows in a single layer of fluid (vertically homogeneous) [152]. In oil spill modeling MIKE 21 has been applied by Liu et al. [42]. Hydrodynamics, advection-dispersion, short waves, sediment transport, water quality, eutrophication, and heavy metals are all simulated by the system's various modules. The system has many engineering and environmental applications, including coastal hydraulics, oceanography, wave dynamics, harbors, rivers, environmental hydraulics, and sediment processes. A comprehensive pre- and post-processing module allows for the analysis and graphical presentation of data and model results [153]. The versatile, interactive menu system simplifies data entry, model input, and program execution.

SHYFEM (Shallow water HYdrodynamic Finite Element Model) is 3D hydrodynamic model, developed by ISMAR-CNR (Institute of Marine Sciences - National Research Council). The model solves the system of primitive equations by vertically integrating them across each vertical layer, using the Boussinesq approximation horizontally and the hydrostatic approximation vertically. It employs the generic ocean turbulence model [154] to calculate vertical diffusivity and viscosity. It was integrated with a transport simulation module and has a spatial resolution ranging from 25 m in extremely coastal areas or shallow waters, to a few kilometers offshore [142,155]. In oil spill modeling SHYFEM has been used by Cucco and Daniel [156], integrating the Lagrangian trajectory and weathering module (FEMOIL) into the operational forecasting system (BOOM), and Ribotti et al. [104], coupling it with the MEDSLIK-II model.

The POSEIDON System was established by the Hellenic Center for Marine Research (HCMR) and provides its hydrodynamic forecasting products [157] also in the MEDESS4MS format to suit the needs of oil spill models for the entire Mediterranean and the Aegean Sea. The POSEIDON hydrodynamical forecasts are released through the implementation of the Princeton Ocean model (POM) [128] with a spatial resolution of 10 km and 24 vertical sigma layers [19,108]. The POSEIDON Mediterranean model provides boundary conditions at the POSEIDON Aegean model’s western and eastern open boundaries [107,157]. Every week, the POSEIDON Aegean Sea model is re-initialized using the HCMR Mediterranean model analysis at 3.5 km horizontal resolution [108].

The CYCOFOS (Cyprus Coastal Ocean Forecasting System) provides hydrodynamic data [158] for the Eastern Mediterranean, covering the Aegean Sea and the Levantine basin. CYCOFOS hydrodynamic model is based on a modified version of POM [128,159] with a spatial resolution of 2 km and 30 vertical sigma layers, while for oil spill models produce 15 vertical z-layers in MEDESS4MS format [108]. The CYCOFOS hydrodynamical model is nested to the Copernicus marine service, while the surface forcing is provided by the ECMWF. CYCOFOS generates daily, 6-hourly mean forecasts for the following four and a half days in a dedicated NetCDF files, designed specifically to cover the needs of oil spill models, known as MEDESS4MS format [108]. The CYCOFOS ocean forecasts have been extensively validated and compared to the parent models, as well as satellite remote SST and in-situ observations [60,158]. The CYCOFOS hydrodynamical forecasts have been used in real oil spill pollution incidents [58,62] and in several oil spill test-cases [59,108,160].

The INGV (Istituto Nazionale di Geofisica e Vulcanologia) provides hydrodynamic data products utilizing the MFS (Mediterranean Forecasting System) model [144,161] and the high-resolution Adriatic Forecasting System (AFS) model [162]. The MFS is based on the NEMO Ocean General Circulation Model (OGCM), which is applied on a model domain with a spatial resolution of 6.5 km. The model spans the whole Mediterranean Sea, and provides daily 10-day forecasts. It has been implemented in many oil spill scenarios, such as the works of Coppini et al. [163] and De Dominicis et al. [19]. The OGCM is linked to a Wave Watch III implementation [164] for the entire Mediterranean Sea, at the same resolution with the hydrodynamics. Furthermore, the Adriatic Forecasting System (AFS) model receives the initial and lateral boundary conditions for temperature, salinity, and velocity from the MFS to produce high resolution (~2.2 km) forecasting outputs for the Adriatic Sea with ECMWF forcing. AFS has been coupled with MEDSLIK-II in oil spill simulations via the De Dominicis et al. [29].

The CNR, Institute for the Marine and Coastal Environment (Naples) (CNR - IAMC) computes hydrodynamics using the Western Mediterranean Model (WMED) [155]. This forecasting system of the marine circulation at sub-regional scale (about 3.5 km) covers the Western Mediterranean, around Sardinia. The model is based on the three- dimensional primitive equation, finite difference hydrodynamic model, named POM [128]. It has been coupled with MEDSLIK-II oil spill model [19,108] and MOTHY model [156].

The IASA provides forecasting ocean products for the Eastern Mediterranean for the next 5 days via the Aegean Levantine Eddy Resolving Model (ALERMO) [165]. This model is a high-resolution implementation of the POM used in the Aegean-Levantine basins, with horizontal resolution of 3.5 km and 25 logarithmic sigma levels in the vertical. The Copernicus Med-MFC (https://marine.copernicus.eu/about/producers/med-mfc) is used to define the one-way nested open boundary conditions [166]. ALERMO has been applied in oil spill modeling as forcing data by Zafirakou-Koulouris et al. [103] and Zodiatis et al. [108,158]. Moreover, the IFREMER (French Research Institute for Exploitation of the Sea) provides oceanographic forecasts, produced by the PREVIMER-MENOR model, which cover the northern part of the Western Mediterranean Sea with 1.2 km horizontal resolution and 60 sigma levels, refined near the surface. The boundary conditions of PREVIMER-MENOR are provided via the Copernicus Med-MFC model. This configuration, which is based on a primitive equation model devoted to regional and coastal modeling, is utilized for both operational and academic reasons [167,168]. The PREVIMER-MENOR model has been coupled to oil spill simulations by De Dominicis et al. [19] and Zodiatis et al. [108].

Table 2. Hydrodynamics models with the corresponding domains and horizontal resolutions.

Hydrodynamics	Provider	Geographic Coverage	Spatial Resolution	Reference
Poseidon Med Model	HCMR	Mediterranean	~ 10 km	[19,108]
Poseidon Aegean Model	HCMR	Aegean Sea	~ 3.5 km	[108]
CYPOM	CYCOFOS	Aegean-Levantine	~ 2 km	[108,160,163]
WMED	CNR IAMC	Western Mediterranean	~ 3.5 km	[19,108,156]
ALERMO	IASA	Eastern Mediterranean	~ 3.5 km	[103,108,158]
MFS	INGV	Mediterranean	~ 6.5 km	[19,108,163]
AFS	INGV	Adriatic Sea	~ 2.2 km	[29,108]
PREVIMER MENOR	IFREMER	North Western Mediterranean	~ 1.2 km	[19,108]
ΜΙΚΕ21	DHI	Regional	-	[42]
CROCO	IRD	Regional	1 km	[68]
NorShelf	Norwegian Meteorological Institute	Norwegian Shelf Sea	2.4 km	[23]
SANIFS	CMCC-OPA	Mediterranean basin	3 km	[56]
SHYFEM	ISMAR-CNR	Regional	4 km, 1km	[104,155]
NEMO	CMEMS	Mediterranean Global Global	(1/24°) ~ 4 km (1/12°) ~ 9 km (1/4°)~27 km	[50,54,55] [51,53,57] [69]
NGOM	NOAA- CSDL	North-eastern and Central GOM	5-6 km	[66]
NCOM	NOAA FNMOC	American Seas and Alaska	3 km	[66,74]
SABGOM	NCSU	GOM	~ 5 km	[66]
IASROMS	NCSU	GOM	~ 2 km	[66]
GLB-HYCOM	NOAA NRL	Global	1/12° (~ 9 km)	[44,70]
GoM-HYCOM	NOAA NRL	GOM	$1/25\text{°}$(~ 4 km)	[66,67]
GoM-HYCOM	NOAA NRL	GOM	1/50° (~ 2 km)	[46,47,49]
Fkeys-HYCOM	NOAA NRL	South Florida coastal, shelf areas and Straits of Florida	1/100° (~ 1 km)	[46,48]

Table 2. Hydrodynamics models with the corresponding domains and horizontal resolutions.

Hydrodynamics	Provider	Geographic Coverage	Spatial Resolution	Reference
Poseidon Med Model	HCMR	Mediterranean	~ 10 km	[19,108]
Poseidon Aegean Model	HCMR	Aegean Sea	~ 3.5 km	[108]
CYPOM	CYCOFOS	Aegean-Levantine	~ 2 km	[108,160,163]
WMED	CNR IAMC	Western Mediterranean	~ 3.5 km	[19,108,156]
ALERMO	IASA	Eastern Mediterranean	~ 3.5 km	[103,108,158]
MFS	INGV	Mediterranean	~ 6.5 km	[19,108,163]
AFS	INGV	Adriatic Sea	~ 2.2 km	[29,108]
PREVIMER MENOR	IFREMER	North Western Mediterranean	~ 1.2 km	[19,108]
ΜΙΚΕ21	DHI	Regional	-	[42]
CROCO	IRD	Regional	1 km	[68]
NorShelf	Norwegian Meteorological Institute	Norwegian Shelf Sea	2.4 km	[23]
SANIFS	CMCC-OPA	Mediterranean basin	3 km	[56]
SHYFEM	ISMAR-CNR	Regional	4 km, 1km	[104,155]
NEMO	CMEMS	Mediterranean Global Global	(1/24°) ~ 4 km (1/12°) ~ 9 km (1/4°)~27 km	[50,54,55] [51,53,57] [69]
NGOM	NOAA- CSDL	North-eastern and Central GOM	5-6 km	[66]
NCOM	NOAA FNMOC	American Seas and Alaska	3 km	[66,74]
SABGOM	NCSU	GOM	~ 5 km	[66]
IASROMS	NCSU	GOM	~ 2 km	[66]
GLB-HYCOM	NOAA NRL	Global	1/12° (~ 9 km)	[44,70]
GoM-HYCOM	NOAA NRL	GOM	$1/25\text{°}$(~ 4 km)	[66,67]
GoM-HYCOM	NOAA NRL	GOM	1/50° (~ 2 km)	[46,47,49]
Fkeys-HYCOM	NOAA NRL	South Florida coastal, shelf areas and Straits of Florida	1/100° (~ 1 km)	[46,48]

3.3. Waves

The ECMWF provides global wave forecasts with a spatial resolution 1/8° using the third generation spectral WAve Model (WAM), [95,169,170]. With 25 frequencies and 24 directions, the WAM model computes the two-dimensional wave distribution. In addition, WAM model has a 1/8° horizontal resolution with outputs every 12-hours. It is the first wave model to solve the complete action density equation, which includes non-linear wave-wave interactions. The WAM model is used operationally in global and regional applications to forecast the sea state. The model may be used for a variety of applications, including ship routing and offshore activities, as well as the validation and interpretation of satellite observations. The output of this wave data product is provided in netCDF format (https://www.ecmwf.int/en/forecasts/datasets/set-ii).

Moreover, the CMEMS, Copernicus Marine System [171] performs operational wave simulations and provides several wave datasets. In this review two datasets are considered: (a) The Global Reanalysis, which is a global wave reanalysis product, describing historic sea states, since 1993. This dataset is based on the WAVERYS model (WAVeReanalYSis) [172], subject to the MFWAM (Météo-France WAve Model) model [173], a third-generation wave model that calculates the directional wave spectrum (i.e., the distribution of sea state energy in frequency and direction) on a 1/5° irregular grid. Average wave quantities derived from this wave spectrum are the significant wave height, the average wave period and the sea surface wave Stokes drift (u and v velocities). These are provided on a regular 1/5° grid, with 3-hour time step. Besides, WAVERYS incorporates oceanic currents from the GLORYS12 physical ocean reanalysis [174] as well as the significant wave height, observed from historical altimetry missions and directional wave spectra from Sentinel 1 SAR, since the beginning of 2017. Detailed information of this dataset is presented by Law-Chune et al. [172]. (b) The main wave product of the Mediterranean Sea Forecasting system (MEDSEA_ANALYSISFORECAST_WAV_006_017) [175], which is made up of hourly wave parameters, with a horizontal resolution of 1/24° that cover the Mediterranean Sea and extend up to 18.125W into the Atlantic Ocean [176]. The wave forecast component (Med-Waves system) is based on the upgraded WAM Cycle 4.6.2. (https://github.com/mywave/WAM). With 24 directional and 32 logarithmically distributed frequency bins, the Med-Waves modeling system resolves the prognostic part of the wave spectrum, and the model solutions are corrected by an optimal interpolation data assimilation scheme of all available along track satellite significant wave height observations. The model employs wave spectra from the GLOBAL ANALYSIS FORECAST WAV 001 027 product (https://data.marine.copernicus.eu/product/GLOBAL_ANALYSISFORECAST_WAV_001_027/description) for Open Boundary Conditions [177,178]. The wave system includes two forecast cycles that provide a Mediterranean wave analysis twice per day and wave forecasts for the next ten days.

The POSEIDON System, established by Hellenic Center for Marine Research (HCMR) [157], provides the sea state forecasting outputs for the Mediterranean and Aegean/Ionian Seas via WAM Cycle 4 model with spatial resolutions of 10 km and 3.5 km, respectively [179]. It is a third-generation wave model capable to compute the spectra of randomly produced short-crested wind waves. The wave forecasting system generates wave forecasts for the next five days, based on hourly analysis and forecasted winds, generated by the POSEIDON weather prediction system.

The Cyprus Coastal Ocean Forecasting and Observing System products (CYCOFOS) [180,181] provides hourly wave forecasting data using the WAM4 model for the Mediterranean and the Black Sea, at 5 km horizontal resolution [19,181], in MEDESS4MS netCDF format to suit the oil spill models.

The UoM provides the wave forecasting data produced by the MALTA Maria WAM model for the central Mediterranean region, with 12.5 km horizontal resolution, using the 3rd generation spectral wave model WAM Cycle 4 [109,182].

The IFREMER provides wave data for the whole Mediterranean, as distributed to the MEDESS-4MS service, through the PREVIMER-MENOR-WW3 model, with 10 km spatial resolution. This wave model is based on the WaveWatch III configuration [183]. Also, it is forced by the atmospheric model of the French Metoffice, the ARPEGE. In parallel, the Puertos del Estado (PdE) (https://www.mitma.gob.es/empresas-fomento/puertos-del-estado ) distributes wave data through the PdE WAM wave forecasting model, with 8 km horizontal resolution, covering the western Mediterranean domain [108].

The Météo-France supports a global forecasting wave system, named MFWAM [184], which is based on the wave model WAM [177]. It is a global model that has 1/12° horizontal resolution and 3-hour instantaneous temporal resolution [83].

Table 3. Wave models and wave products with the corresponding domains and horizontal resolutions.

Wave system	Provider	Geographical area	Spatial Resolution	Data Type	Reference
Poseidon WAM Cycle 4 Med	HCMR	Mediterranean	~ 10 km	Forecast	[108,179]
Poseidon WAM Cycle 4 Aegean	HCMR	Aegean	~ 3.5 km	Forecast	[108,179]
WAM4	CYCOFOS	Mediterranean and Black Sea	~ 5 km	Forecast	[19,108,181]
PdE-WAM	PdE	Western Mediterranean	~ 8 km	Forecast	[108]
PREVIMER-MENOR-WW3	IFREMER	Mediterranean	~ 10 km	Forecast	[19,108]
MALTA/Maria WAMI	UOM	Central Mediterranean	~ 12.5 km	Forecast	[108,109]
WAM Cycle 6, WAM 4.6.2	CMEMS	Mediterranean	1/24° (~ 4.5 km)	Forecast/Reanalysis	[50,176]
MFWAM	Meteo-France	Global	1/12° (~9 km)	Forecast	[83,177]
WAVERYS	CMEMS	Global	1/5° (~ 22 km)	Reanalysis	[51,172]
WAM	ECMWF	Global	0.125° (~ 13 km)	Forecast	[23,46,47,48,49,185]

Table 3. Wave models and wave products with the corresponding domains and horizontal resolutions.

Wave system	Provider	Geographical area	Spatial Resolution	Data Type	Reference
Poseidon WAM Cycle 4 Med	HCMR	Mediterranean	~ 10 km	Forecast	[108,179]
Poseidon WAM Cycle 4 Aegean	HCMR	Aegean	~ 3.5 km	Forecast	[108,179]
WAM4	CYCOFOS	Mediterranean and Black Sea	~ 5 km	Forecast	[19,108,181]
PdE-WAM	PdE	Western Mediterranean	~ 8 km	Forecast	[108]
PREVIMER-MENOR-WW3	IFREMER	Mediterranean	~ 10 km	Forecast	[19,108]
MALTA/Maria WAMI	UOM	Central Mediterranean	~ 12.5 km	Forecast	[108,109]
WAM Cycle 6, WAM 4.6.2	CMEMS	Mediterranean	1/24° (~ 4.5 km)	Forecast/Reanalysis	[50,176]
MFWAM	Meteo-France	Global	1/12° (~9 km)	Forecast	[83,177]
WAVERYS	CMEMS	Global	1/5° (~ 22 km)	Reanalysis	[51,172]
WAM	ECMWF	Global	0.125° (~ 13 km)	Forecast	[23,46,47,48,49,185]

4. Analysis of Selected Oil Spill Models, in terms of Boundary Forcing

Twenty-four oil spill modeling test cases have been selected for the examination of the boundary conditions applied in terms of winds, hydrodynamics and waves. These test cases are summarized in Table 4.

4.1. Wind Fields

All selected oil spill models used operational wind data as input, except for one case, where the wind was constant throughout the duration of the simulation. The majority of oil spill models (16 case studies) utilize wind data products received from the ECMWF (ERA5 or Interim). Moreover, in many test cases, the NOAA wind data products have been (e.g., NCEP-NARR, CFSR, NAM, NO GAPS, WRF, GFS). Apart from the study of Nugroho et al. [68], which uses data from NCOM AMSEAS, the Ribotti et al. [104] use wind data obtained from the SKIRON model, while the Grubesic & Nelson [74] applied a constant wind speed and direction in their oil spill scenario.

The OpenOil model was implemented in seven of the selected test-cases, and the majority of these studies uses as wind forcing data from the ECMWF, with variable resolutions (0.125°, 0.25° and 0.1°), as shown in Table 4. The OpenOil implementation for the Thracian Sea uses wind data from NOAA-GFS, with resolution 0.25°. Regarding the Medslik-II, all applications use wind data from the ECMWF, with 0.125° (~ 13 km) horizontal resolution. Kampouris et al. [55] imposes the ECMWF wind data with $~9kmand~18km$ spatial resolution, while the Ribotti et al. [104] study used wind data from the SKIRON model, with 10 km spatial resolution. Besides, the MEDSLIK applications used wind data from the SKIRON model, with 5 km horizontal resolution. Concerning the GNOME oil spill model, the wind data were obtained from either the ECMWF model, with resolution 1/8°, or from the NOAA-GFS, with a resolution 0.25°. In both oil spill simulation cases that implemented SIMAP and OILMAP DEEP models, the wind data was derived from the NOAA products. Indeed, four different wind reanalysis and forecast products covering the northeastern GoM, obtained from National Oceanic and Atmospheric Administration (NOAA) and the US Navy government websites were used as model inputs (NAM, NARR, CFSR, NOGAPS). Finally, the BLOSOM, STFM, MIKE 21/3 OS FM model applications were mostly coupled to the NCOM AMSEAS model, with 1/36° (~ 3 km) spatial resolution, to the WRF model, with 0.15° resolution, and to the ECMWF data products, with 0.125° resolution.

4.2. Hydrodynamics

All examined oil spill model implementations used hydrodynamic data, to transfer oil as passive pollutant. Most applications (11 test cases) used datasets obtained from CMEMS and NOAA. The high-resolution hydrodynamic HYCOM was applied only for the Gulf of Mexico. The HYCOM-FKEYS model, with 1/100° resolution, was applied in two oil spill simulation studies in the Gulf of Cuba. Two studies utilized the HYCOM-GoM configuration, with 1/25° resolution, and one the Global-HYCOM, with 1/12° resolution. Furthermore, five high-resolution hydrodynamic models were implemented using NorShelf, SANIFS, MIKE 21, CROCO and SHYFEM (Table 4).

Focusing on the OpenOil model, of the 7 studies examined, the high-resolution HYCOM-GoM 1/50° model was used in four cases [46,47,48,49], and the HYCOM-FKEYS with 1/100° resolution, in two test cases [46,48]. Also, the model was coupled with the CMEMS data products, either the Global 1/12C [51], or the Mediterranean, with 1/24° resolution [50]. Moreover, the NorShelf model, with 2.4 km spatial resolution, covering the hydrodynamics of the Norwegian Sea was also used by Röhrs et al. [23]. The examined cases in which MEDSLIK-II oil spill model was implemented, in almost all cases, the hydrodynamic forcing data used were obtained from the CMEMS Global 1/12° configuration [53,57], and the CMEMS Mediterranean 1/24° data product [54,55,58,61]. Liubartseva et al. [56] used the SANIFS high-resolution hydrodynamic model for oil spill simulation studies off the coast of Italy, while Ribotti et al. [104] integrated the SHYFEM high-resolution hydrodynamic model, with 4 km and 1 km spatial resolution. The MEDSLIK applications used hydrodynamic data from CYCOFOS, with 2 km spatial resolution [58,61]. At the same time, GNOME users developed oil spill simulations using variable hydrodynamic forcings: (a) the GLO-CPL (1/4°) [69], (b) the CMEMS-Global (1/12°) [68], (c) the GLB-HYCOM (1/12°)[70], and (d) the high-resolution ocean model CROCO (1km) [68], in which the wind was constant. SIMAP model has been coupled to various hydrodynamic forcing data, with various resolutions, such as the GoM-HYCOM 1/25°, GoM-HYCOM Reanalysis (1/25°), GoM-HYCOM Real-time (1/25°), SABGOM (5 km), IAS ROMS (6 km), NCOM Real Time (3km) and NGOM (5-6 km), by French-McCay et al. [66,67].

Eventually, oil spill models like BLOSOM, STFM and MIKE 21 OS FM used ocean currents from the NCOM-AMSEAS (1/36°) [74], the GLB-HYCOM 1/12° [44] and the high-resolution hydrodynamic model MIKE21 [42].

4.3. Waves

In the majority of oil spill simulations, the impact of waves is not taken into account. As shown in Table 4, nine case studies do not couple with wave models, while only three applications implemented the empirical JONSWAP wave spectrum and Ekman influence. Besides, in six studies, wave data were integrated through the ECMWF and the WAM model outputs, running with 0.125° horizontal resolution [23,46,47,48,49]. In two case studies, wave data were obtained from the forecasting CMEMS products, with WAVERYS model at 1/5° resolution [51] and WAM Cycle 6 model at 1/24° resolution [50]. Based on the above, we may conclude that a high-resolution wave model has never been used, in conjunction with high-resolution hydrodynamic and meteorological models in the operational oil spill modeling works of the latest years.

Focusing specific implementations per case study, it occurs that in all 7 examined OpenOil model cases, wave data were taken into account. The majority of these studies forced the oil spill model with data from the ECMWF and the WAM. In parallel, it is noticed that in two applications, waves were obtained from the CMEMS (Global 1/12° and Mediterranean 1/12°) data products. This means that OpenOil model was never coupled to any wave model with spatial resolution higher than 1/24°. Following Table 4, it is obvious that all eight MEDSLIK-II model cases did not take into account the operational wave data. More specifically, in some studies, the waves were not considered as a mechanism affecting the spillage [53,104], while the empirical JONSWAP wave spectrum and constant Ekman values were used in the works of Liubartseva et al. [54,56], Siqueira et al. [57], and Kampouris et al. [55]. Moreover, concerning the GNOME model, in 3 out of 4 cases, no wave data were used as boundary input. The studies working with SIMAP, OILMAPDEEP, MEDSLIK, BLOSOM, STFM and MIKE OS FM models, did not take into account wave data in their simulations.

4.4. Biogeochemical model

The majority of oil spill test-cases take into account only the main oil weathering processes, such as spreading, evaporation, emulsification, dispersion and beaching. Furthermore, from the 24 applications, only in 4 cases the process of biodegradation was included, in 2 cases dissolution and sedimentation were considered and in one the photo-oxidation process was integrated.

In the OpenOil simulations (5 in 7 cases), the physico-chemical processes integrated are evaporation, emulsification, dispersion, beaching, vertical mixing and resurfacing. In only two test-cases [50,51] the biodegradation process is included. However, the biodegradation in these cases takes into account only temperature as input parameter, so there is a need for further improvement. In addition, in all case-studies of MEDSLIK models, the physico-chemical processes included are evaporation, emulsification, dispersion, spreading and beaching. It is worth mentioning that MEDSLIK-II [18,186] included biogeochemical processes, but none of these simulations was applied to actual or hypothetical scenarios. As for GNOME model, in almost all test-cases (3 of 4 cases), evaporation, emulsification, dispersion, spreading and beaching are obtained, except from one case [68] which did not consider any OWP. The most comprehensive oil spill model in terms of biochemical processes is SIMAP, which takes into account all OWPs, such as evaporation, emulsification, dispersion, spreading, beaching, dissolution, sedimentation, biodegradation and photo-oxidation. However, in both test-cases using SIMAP, no one was coupled with a wave model. Concerning the STFM and MIKE 21 OS FM test simulations, only evaporation, emulsification and dispersion was included, while the BLOSOM simulations involved evaporation, emulsification, dispersion, spreading and beaching, although the model additionally supports the dissolution process.

Table 5. 7-Case studies of OpenOil model (2018-2022) according to Metocean data and initial satellite observations.

Time period	SAR data	Currents	Waves	Wind	Sal/Temp	Reference
22/2/11-4/3/11 3-13/1/2011 5-15/1/2011 7-17/1/2011 (10 days)	No	NorShelf model (2.4 km)	ECMWF (0.125°)	ECMWF (0.125°)	No	[23]
20-27/5/2010 (7-days) 2-10/6/2010 (8 days)	Yes	GoM-Hycom (1/50°)	ECMWF (0.125°)	$\mathrm{ECMWF}$(0.125°)	Only Salinity - GoM-HYCOM (1/50°)	[47]
2011-2016 (6 years)	No	GoM-Hycom (1/50°)	ECMWF (0.125°)	ECMWF (0.125°)	No	[46]
Feb. 2017 (5 days) August 2017 (5 days)	Yes	GoM-Hycom (1/50°)	ECMWF (0.125°)	ECMWF (0.125°)	Only Salinity - GoM- HYCOM (1/50°)	[49]
2010-2017 (7 years)	No	GoM-Hycom (1/50°)	ECMWF (0.125°)	ECMWF (0.125°)	GoM- HYCOM (1/50°)	[48]
25-30/10/2020 (5 days)	No	CMEMS- Mediterranean Forecast (1/24°)	CMEMS- Mediterranean Forecast (1/24°)	$\mathrm{NOAA}-\mathrm{GFS}(0.25^{\mathrm{o}}$)	CMEMS- Mediterranean Forecast (1/24°)	[50]
20-23/07/2014 (3 days) 21-23/08/2014 (2 days)	No	CMEMS-Global Forecast (1/12°)	CMEMS-Global Reanalysis (1/5°)	$\mathrm{ECMWF}-\mathrm{ERA}5\left(0.25^{\mathrm{o}}\right)$	No	[51]

Table 5. 7-Case studies of OpenOil model (2018-2022) according to Metocean data and initial satellite observations.

Time period	SAR data	Currents	Waves	Wind	Sal/Temp	Reference
22/2/11-4/3/11 3-13/1/2011 5-15/1/2011 7-17/1/2011 (10 days)	No	NorShelf model (2.4 km)	ECMWF (0.125°)	ECMWF (0.125°)	No	[23]
20-27/5/2010 (7-days) 2-10/6/2010 (8 days)	Yes	GoM-Hycom (1/50°)	ECMWF (0.125°)	$\mathrm{ECMWF}$(0.125°)	Only Salinity - GoM-HYCOM (1/50°)	[47]
2011-2016 (6 years)	No	GoM-Hycom (1/50°)	ECMWF (0.125°)	ECMWF (0.125°)	No	[46]
Feb. 2017 (5 days) August 2017 (5 days)	Yes	GoM-Hycom (1/50°)	ECMWF (0.125°)	ECMWF (0.125°)	Only Salinity - GoM- HYCOM (1/50°)	[49]
2010-2017 (7 years)	No	GoM-Hycom (1/50°)	ECMWF (0.125°)	ECMWF (0.125°)	GoM- HYCOM (1/50°)	[48]
25-30/10/2020 (5 days)	No	CMEMS- Mediterranean Forecast (1/24°)	CMEMS- Mediterranean Forecast (1/24°)	$\mathrm{NOAA}-\mathrm{GFS}(0.25^{\mathrm{o}}$)	CMEMS- Mediterranean Forecast (1/24°)	[50]
20-23/07/2014 (3 days) 21-23/08/2014 (2 days)	No	CMEMS-Global Forecast (1/12°)	CMEMS-Global Reanalysis (1/5°)	$\mathrm{ECMWF}-\mathrm{ERA}5\left(0.25^{\mathrm{o}}\right)$	No	[51]

5. Conclusions and Recommendations

Nonetheless, numerical models will always have limitations and uncertainties. Waves, for example, have a tendency to travel towards the beach in shallow water close to the shoreline, which could result in the movement of oil onto land. However, not all waves move oil the same way. Waves from distant hurricanes can spread oil much further than ordinary local wind-driven waves. Furthermore, the presence of an oil slick on the water's surface can alter how air flows over it, altering the manner the wind propels the waves. This is because oil slicks change both the rate of evaporation at the water's surface and the force of the wind on the water. Furthermore, oil acts as a physical barrier on the water's surface. With this barrier in place, less water will evaporate in areas with surface oil. The rise in water temperature caused by the decrease in evaporation increases heat flow, which influences surface winds, waves, and oil movement [187]. Oil also reduces friction on the water's surface, resulting in an increase in near-surface wind speed, thus decreasing the wind force exerted on the water's surface. Observations show that surface winds respond to variations in sea surface temperatures. According to modeling studies, this response can happen very quickly, and the ocean then reacts to the change in winds. If responders are to act quickly and effectively in the event of an oil spill, they should have access to relatively accurate models and be aware of their limitations [15]. A trajectory model that is not continuously updated with field measurements, after the initial days of spillage, will be unreliable.

Model constraints can lead to errors and uncertainty in predictions. Constraints on model performance include model resolution, input accuracy, and the ability to accurately describe physical processes such as winds, waves, and currents. When these constraints and uncertainties are removed, the capacity for model prediction will improve. It is critical in operational modeling to understand both the predictions produced by numerical models and any potential inaccuracies [15]. Ambiguities occur in ocean circulation models and oil spill models, so some information may be inaccurate or incomplete. The primary source of uncertainty in ocean circulation modeling and forecasting will vary depending on location and specific mechanisms [19,188]. The main sources of uncertainty in operational oil spill modeling are the drivers, which include the source and amount of oil, as well as wind, wave, and current projections. If these uncertainties are not corrected, then uncertainties are transferred to the oil spill model results. Due to these and other sources of uncertainty, operational oil spill forecasting is only done for short periods of time, typically ranging from 12 hours to two days. Using multiple models allows for output comparison, which can improve the forecast accuracy for weekly time frames. Uncertainty is reduced with each new modeling cycle by combining oil observations with meteorological and oceanic variables [15] and evaluating and correcting model results with satellite observations.

This review highlights that the biggest deficiency in oil spill models at the moment lies in the coupling to high resolution wave models, since in most test cases wave data were not included as forcing or only coarse resolution wave fields were integrated in simulations. Moreover, it is found that some biochemical processes, like biodegradation, sedimentation, dissolution, photo-oxidation, are not included in the operational oil spill scenarios. In the limited models that these OWPs were considered, the algorithms used were fairly simplistic, thus jeopardizing the accuracy of the results.

In order to reduce the model uncertainty, additional steps must be taken in the future, improving the accuracy of the models' representations in terms of chemistry, physics, and ocean dynamics. Thus, operational ocean and oil spill modeling should be improved, creating multiple models and combining them to produce better results. This multi-model oil spill modeling approach, as shown by the MEDESS4MS project [18,108] can provide confidence to the response agencies, especially if they have to act offshore.

Furthermore, uncertainties in ocean, meteorological, and wave modeling should be assessed and communicated to the users of these models. Responders and the general public would benefit from the user-centered design in understanding both the basic information about oil spill trajectories and the corresponding deficiencies and uncertainty levels.