Lagrangian Model PETROMAR-3D to Describe Complex Processes in Marine Oil Spills

Amilcar E. Calzada¹, Iván Delgado¹, Carlos Ramos¹, Frank Pérez¹, Dailín Reyes¹, Dayana Carracedo¹, Alejandro Rodríguez¹, Dayron Chang¹, Javier Cabrales², Alexander Lobaina³

¹Center for Marine Meteorology, Meteorology Institute, Havana, Cuba.
²Center for Atmospheric Physics, Meteorology Institute, Havana, Cuba.
³Air Information System on Meteorology, The Cuban Company of Air Navigation, Havana, Cuba.

DOI: 10.4236/ojms.2021.111002                

Abstract

The version 2.1 of PETROMAR-3D model, created in the Center for Marine Meteorology of the Meteorology Institute of Cuba, is presented. This Lagrangian model has been designed to describe the physical processes of marine oil spills in the face of multiple scenarios of the marine environment. Although it is applicable to any part of the world, it is mainly designed for inter-American seas. The novelty has been to integrate the processes of drift and weathering into a model, with updated methods that incorporate 3D phenomena, a very favorable situation to achieve an operating system in Cuba and the region for the immediate and medium term. Python was chosen as the programming language because it has advanced libraries for numerical modeling, automation work and other useful tools for pre-and post-processing. By means of adapters, an important number of atmospheric, hydrodynamic and wave models have been considered to create the scenarios efficiently. The modular distribution in which the code has been created facilitates its use for other dispersion analysis and biophysical applications. Finally, a set of simple images are presented, aimed at informing decision-makers in order to mitigate the effects of the spill on the environment.

Keywords

Oil Spill Modeling, Lagrangian Model, Inter-American Seas

1. Introduction

The presence of Exclusive Economic Zone (ZEE as abbreviated in Spanish) of the countries that make up the Gulf of Mexico, together with the transit through its territorial waters of millions of barrels of oil daily, has promoted scientific research on the processes involved during a spill. In this regard, one of the efforts made is the creation of the PETROMAR-3D model, whose system is ready to be available in the region of the inter-American seas.

This system has been used in Cuba for the main practical exercises in the face of oil prospecting by some international companies. Additionally, it was used in the event of contingencies such as the Macondo well, an event that not only made profound damages to environment in several ways, but also generated countless important investigations in the region [1] in order to find effective means to mitigate future events like this.

The main antecedents of this work are a national scientific project completed in 2009 and the experience gained from the Gulf Spill in 2010, with own evaluated daily forecasts made for Ministry of Technology and Science functionaries and other national decision-makers. Version 1.0, 1.1 and 1.2 of a drift model, designed and implemented by the authors named PETROMAR [2], are the nearest precedents of the model purposed by the present paper.

The main outputs of the last project, performed by the team and finished in November 2017, are presented in this article. It was initially proposed to update and improve the algorithms and methods of previous versions of the model, with the introduction of the weathering processes and description of the blowout, to obtain more real and objective scenarios.

2. Materials and Methods

The main update made is related with the three-dimensional conception of this version of the model, expressed in its modules that describe the ascent of the pollutant from the depths to the surface and weathering processes. Other improvements were made to the drift module, so the three modules work jointly and in harmony.

2.1. General Model Design. Physical and Mathematical Basis of the Model

In the calculation of drift and weathering of the slick there are, on the one hand, many empirical functions that are still used with well-defined parameters and, on the other hand, numerical models that apply equations of fluid dynamics. Until now, both tendencies have elements of significance and that is why neither of them should be discarded as a possible solution [3]. The following sub-paragraphs will explain the methods used in the three main modules and in the complementary modules, related to pre-processing and post-processing. Any physical process of those considered in the oil slick affects the development of others, an issue that propitiates that the modelers shall be careful in developing their algorithms [4]. Thus, this version 2.1 will describe the physical processes as follows:

- If the spill occurs on the sea surface, the runoff is calculated according to Fay’s equations [2] [5], to then decompose the slick into virtual particles and execute the advection and diffusion processes. The advection is calculated considering the marine currents and the wind. For the diffusion, a Brownian motion is simulated using the random walk method [3]. Taking these processes into account, the movement of partial particles was assumed using three numerical schemes: Euler (E), Rungge-Kutta of second order (RK2) and fourth order (RK4); each one with its features, accuracies and requirements. Interaction with the coastline is analyzed in a first approximation because there is no differentiation of the type of coastline. Despite this, there is a high accuracy in the calculation of the time in which each of the parts making up the pollutant arrives at the coast.

- At the same time, Evaporation processes (with the PC method), vertical dispersion (with Delvigne and Sweeney’s proposal) and emulsification (through Eley’s law of change) are calculated (from the beginning). All this combination of processes will favor, for the first time, that time based parameters such as density, viscosity, thickness and concentration of the oil slick could be calculated.

- If the spill takes place in the sea depths, everything will be executed in a similar way, except for two main aspects: the spreading of the oil on the surface will be the result of the value of the radius of the plume around its axis. The parameters of the slick’s position and dimension are described in the blowout process. Other aspect is that the values of density and viscosity of the pollutant in the surface movement will be linked to the upward movement, introduced by the method proposed by Zheng and Yapa [6].

The system will do the above in four stages:

1) Data selection and interpretation.

2) Creation of the scenarios (framed by dates, zones and supplier models, starting with a configuration file that allows the entry of initial parameters, see Annex).

3) Calculation of the trajectory considering the initial coordinates (X, Y, Z), the volume spilled or the rate of spillage, the type of oil and other aerial or satellite information are also initialized within the configuration file (see Annex).

4) Graphic representations and storage of the numerical data of the results are tasks carried out through the Matplotlib and Basemap libraries.

The Lagrangian method to describe the processes in the slick causes the concentration to be calculated indirectly. This process is done considering the amount of existing virtual particles (with each of its macro-scalar features) in a pre-determined control volume (hereafter CV) around each node of the mesh, as well as its macro-scalar properties that are modified by the environment. Most of the models developed in the first decade of the 21^st century and even during the first years of the second, present similar ideas. Despite this, there is a novel idea to carry out the reconstruction of the oil concentration field from the distribution process of advection-diffusion and transformation processes of the lagrangian particles [7] [8]. At the same time have been some ideas to implement a eulerian approach around Reynolds-averaged Navier-Stokes (RANS) equations, combined with a Lagrangian particle dispersion model [9]. However, the model proposed just will handle some of these ideas in the next versions.

2.1.1. Upward Oil Movement. Description of the Blowout

This model is developed based on the Lagrangian concept proposed by Lee and Cheung [10]. The considerations taken into account for implementing the model (Figure 1) are the following:

1) The plume cross-section is circular and perpendicular to the trajectory.

2) The control volume has the shape of an inclined cone.

3) The effects of oil viscosity are not included in the calculations due to the turbulent plume development.

4) Forced penetration of ambient fluid into the plume occurs from its front.

5) The model is of the quasi-stationary type.

1) Governing plume equations

In this work, a Lagrangian element is defined as a control volume moving through the central plume line. Element moves with the speed V of its central line. Thickness and mass are defined by $h=V\text{d}t$ and $m=\rho \pi b^{2}h$ , been $\rho$ and b the density and radius of the volume respectively. The following governing equations are applied to the control volume.

a) Conservation of mass equation for a plume

$\frac{\text{d}m}{\text{d}t}={\rho }_a{Q}_e-{\displaystyle {\sum }_i^n\left(\frac{\text{d}{m}_i}{\text{d}t}-\frac{\text{d}{m}_d}{\text{d}t}\right)}$ (1)

Figure 1. Schematic diagram of the upward movement: (a) graphical representation of the control volumes in the plume over the upward movement of the oil from horizontal perspective, and (b) distributions of the oil and gas drops in the vertical section.

In the Equation (1), ${\rho }_a$ is the density of the ambient fluid, $Q_e=Q_f+Q_s$ the volume flow entering the control volume, due to forced penetration $Q_f$ and induced shear penetration $Q_s$ . The second term is the summation of the differences of the loss of mass due to dissolution and turbulent diffusion.

In practice, it is convenient to define the term $Q_e$ as a result of a constant of proportionality E multiplied by an environmental function.

$Q_e=E\cdot f\left({\rho }_a,{\rho }_p,v_a,v_p\right)$ (2)

This constant is empirical and must be adjusted through experiments. The aforementioned function f in (2) is proportional to ambient density ( ${\rho }_a$ ), CV density ( ${\rho }_p$ ), ambient speed ( $v_a$ ) and CV speed ( $v_p$ ).

b) Momentum conservation

$\frac{\text{d}\left(mV\right)}{\text{d}t}=V_a\frac{\text{d}m}{\text{d}t}+m\frac{\Delta \rho }{\rho }gk-2\rho bh{C}_D{\left(\left|V\right|-{V^{\prime }}_a\right)}^2\frac{\left|V\right|}{\left|{V^{\prime }}_a\right|}$ (3)

In the above Equation (3), $V_a$ is the average vector of the flow velocity over the exposed plume surface while ${V^{\prime }}_a$ is a projection of $V_a$ in the direction of $V$ and k is a unitary vector in the vertical direction. On the other hand, the deficiency of density between oil and the ambient fluid ${\Delta }_p$ is presented as the difference ${\Delta }_p=p_a-p$ and $C_D$ is a drag coefficient.

The first term on the right in the equation represents the momentum of the mass penetrating the plume, the last two terms represent the net forces acting on a control volume.

c) Conservation of scalar magnitudes

$\frac{\text{d}\left(mI\right)}{\text{d}t}=I_a\frac{\text{d}m}{\text{d}t}-2a\pi \rho bhk\frac{I-I_a}{b}$ (4)

In this case, I is a symbol representing any of the scalar magnitudes that characterize the plume, for heat $I=c_{p}T$ , for salinity $I=S$ and for oil mass concentration $I=C$ . Diffusive constants, such as heat and salinity diffusivity or oil concentration, are represented by k.

In the case of concentration C, an additional term is necessary to take into account the dissolution rate

${\displaystyle {\sum }_i^n\frac{\text{d}{m}_i}{\text{d}t}}$ . (5)

This continuity equation states that scalar changes in the control volume are equal to the scalar magnitude that mainly penetrates through the mass entering the control volume, and the loss of this magnitude due to turbulent diffusion.

d) Equation of state

The set of governing equations needed to model the behavior of a submerged plume can be completed with an equation of state. Its purpose is to describe the changes in density due to salinity, temperature and concentration taking the general form

$\rho =\rho \left(T,S,C\right)$ . (6)

Exact shape of the above equation depends on the type of fluid in the plume, for example, the density of hot water released into seawater depends on the temperature and salinity of the mixed water.

Model selected includes the main processes affecting plume development and correctly represents the physical behavior of the oil that ascends in a broad spectrum of initial and environmental conditions. Because of its formulation, it ceases to be valid once the vertical velocity is zero or very close to it, known as the terminal level (which may be found below the surface). In these cases, it is necessary to take into account other processes that affect oil, such as horizontal spreading.

e) Implementation details

Two classes were implemented to programmatically encapsulate the methods used in this model. One for the numerical processing of the variables of interest on the control volume; the second was implemented for the visualization of the results providing a wide variety of graphical outputs to facilitate the analysis and verification of the results.

The model was implemented in the Python programming language version 2.7, using the integrated development environment Spyder 3.0 (Scientific Python Development Environment). Python Seawater module [11] was used for seawater pressure and density calculations and was implemented following the standards established by UNESCO in 1983 [12].

In order to plot the results, a Python module was implemented where MATPLOTLIB library [13] was used.

2.1.2. Drift Module

The drift processes are governed by the equations

$\frac{\text{d}{x}_i}{\text{d}t}=u_a\left(x_i,t\right)+u_d\left(x_i,t\right)$ (7)

and

$\frac{\text{d}{y}_i}{\text{d}t}=v_a\left(y_i,t\right)+v_d\left(y_i,t\right),$ (8)

considering the motion of a particle i in each axis. The $u_a$ and $v_a$ are the components of the speed of the particle, consequence of the advection and in the other hand, $u_d$ and $v_d$ are those referring to diffusion. Thus, putting equations (6) and (7) in perspective, you would obtain that the advection part could be modelled like

$u_a\left(x_i,t\right)=c_{cor}{u}_{cor}\left(x_i,t\right)+c_{vto}{u}_{vto}\left(x_i,t\right)$ (9)

and

$v_a\left(y_i,t\right)=c_{cor}{v}_{cor}\left(y_i,t\right)+c_{vto}{v}_{vto}\left(y_i,t\right),$ (10)

in which $u_{cor}$ and $v_{cor}$ are the current components and $u_{vto}$ and $v_{vto}$ the wind ones. Coefficients $c_{cor}$ and $c_{vto}$ represent the interaction between the current and the wind with the oil slick. They were selected, taking into account the values obtained in tropical laboratories exercises according to Betancourt [14].

Furthermore, if we incorporate a random walk method, such as the Brownian motion, into these formulations, the position of the particles at any time would be now described by the stochastic differential equations

$\begin{array}{c}{x}_i^{n+1}=x_i^n+\left(c_{cor}{u}_{cor}\left(x_i^n,t\right)+c_{vto}{u}_{vto}\left(x_i^n,t\right)\right)\text{Δ}t\\ +\sqrt{12D_h\text{Δ}t}{\left[R_1\right]}_0^1\cos \left(2\pi {\left[R_2\right]}_0^1\right)\end{array}$ (11)

$\begin{array}{c}{y}_i^{n+1}=y_i^n+\left(c_{cor}{v}_{cor}\left(y_i^n,t\right)+c_{vto}{v}_{vto}\left(y_i^n,t\right)\right)\text{Δ}t\\ +\sqrt{12D_h\text{Δ}t}{\left[R_1\right]}_0^1\sin \left(2\pi {\left[R_2\right]}_0^1\right)\end{array}$ (12)

In these formulations, $D_h$ represents the horizontal diffusion coefficient, choosing 10⁵ m²/s for the calculations [15]. R₁ and R₂ are independent random variables within the $\left[0,1\right]$ interval.

The spills described by the system could be instantaneous or continuous (considering the type of source), on the surface or subsurface (taking into account the depth) and in a fixed location.

1) Implementation details

This module has many details from its previous versions. The programming language selected for this version to re-elaborate the codes of the algorithms was Python v2.7.6. This language covers the majority of the entire code and serves as a coordinating element between different scientific libraries used, mostly written in FORTRAN and C, languages designed for high performance. Numpy, Scipy, Matplotlib, Trapezoid Map Finder, Basemap, Cmocean and Seawater are just examples of them.

Precisely, the ability of Python to interact with other languages to optimally solve specific problems, is one of the main arguments that support its choice as the main language. This feature makes the drift module have an efficient implementation.

On the other hand, the fact that Python and the mentioned scientific libraries are supported by multiple platforms, makes this implementation highly subjected to be implemented in them. The current version was implemented for Linux operating system (Ubuntu 14.04 or superior to be more exact) to facilitate the automation processes.

2.1.3. Mass Transfer or Weathering Module

Three mass transferring processes will be explained below. Two of them (evaporation and vertical dispersion) reduce the initial volume of the pollutant; the first transfers it to the atmosphere and the second to the water column underlying the spill. The third, emulsification, occurs when to water column does not enter any more oil molecules (saturation state) resulting in the mixture. Therefore, the latter increases the initial volume and makes the modeling more complex by giving rise to a new substance with very complex characteristics. Details are given below.

1) Evaporation

Evaporation is the most important mass removal process in an oil spill. It is controlled by the mass transfer coefficient, wind speed, oil diffusivity (Schmidt number) and vapor pressure; although other authors also include water temperature.

The authors of this paper choose the pseudo-component (hereafter PC) model to describe this process. Therefore, the hydrocarbon mixture is divided into several fractions from boiling point data. Gros et al., Kotzakoulakis and George, and Spanoudaki et al. [16] [17] and [18] mention that there are many PC models for measuring the process, but all, in general, are modifications of the models previously developed by Mackay and Matsugu [19], and Mackay et al. [20]. Each PC

is constructed from the first-order approximation of $\frac{\text{d}T}{\text{d}f}$ using empirical correlation based on its API.

$T_0=457.16-3.3447\text{API}$ (13)

$\frac{\text{d}T}{\text{d}f}=1356.7-247.36\ln \text{API}$ (14)

First five PCs of equal volume are created and their boiling points are determined by

$B{P}_i=T_0+\frac{\text{d}T}{\text{d}f}\frac{i-\frac{1}{2}}{5}$ . (15)

Vapor pressure of each PC is based on the Antoine equation [21] such that

$\ln \frac{P_i}{p^0}=\frac{\text{Δ}{S}_i{\left(B{P}_i-C_{2,1}\right)}^2}{\text{Δ}{Z}_{b}RB{P}_i}\left[\frac{1}{B{P}_i-C_{2,i}}-\frac{1}{T-C_{2,i}}\right]$ (15)

in which

$C_{2,i}=0.19B{P}_i-18$ (16)

and

$\text{Δ}{S}_i=8.75+1.987\log B{P}_i.$ (17)

The relative molar volume ${\bar{V}}_i$ of each PC is obtained from a correlation between the molar volume and the boiling point for a series of normal alkanes (C3 to C20).

${\bar{V}}_i=7.0e^{-5}-\left(2.102e^{-7}B{P}_i\right)+\left(9.494e^{-7}{\left(B{P}_i\right)}^2\right)$ (18)

Effective molecular weight of each PC is obtained in the same way.

$M{W}_i=0.04132-\left(1.985e^{-4}B{P}_i\right)+\left(9.494e^{-7}{\left(B{P}_i\right)}^2\right)$ (19)

Reference mass transfer coefficient $K_i^0$ is defined at a wind speed of 1 m/s, and a pool diameter of 1 m. A single Schmidt number value is used for all components and is based on the hydrocarbon moles weighted average.

$K_i^0=0.0048\left(1.3676{\left(\frac{0.018}{M{W}_{avg}}\right)}^{\frac{1}{2}}\right)$ (20)

$M{W}_{avg}=\frac{{\displaystyle {\sum }_{i=0}^{ncp}\frac{M{W}_i{V}_i\left(0\right)}{V}}}{{\displaystyle {\sum }_{i=0}^{ncp}\frac{V_i\left(0\right)}{V}}}$ (21)

In addition, the exposure evaporation matrix is constructed consisting of four columns: exposure evaporation, volatility factor Z, evaporated fraction $F_2$ and vapor pressure.

Volatility factor is a measure of how volatile the PC is and is defined by

$Z\left(t\right)={\displaystyle {\sum }_{i=0}^{ncp}\frac{P_i{V}_t}{RTN\left(t\right)}}$ . (22)

The evaporation exposure values are selected in such a way that the separation between them increases with the index of each PC. The minimum is set to zero, and the maximum was obtained using the maximum time and maximum wind speed.

${\varphi }_{\max }=t_{\max }{K}^0{U}_{\max }^{\frac{7}{9}}$ (23)

The incremental change in evaporation due to exposure also increases with the index of each PC as follows.

$\text{Δ}\varphi \left(j\right)=\frac{2{\varphi }_{\max }j}{j_{\max }^2}$ (24)

The differential equations describing evaporation are then solved in these preset values of evaporation due to exposure.

$\frac{\text{d}{V}_i}{\text{d}\varphi }=-\bar{V}\left(\varphi \right)\left(\frac{P_i}{RT}\right)V_i\left(\varphi \right)$ (25)

A numerical approximation in terms of $\ln V$ is used:

$V_i\left({\varphi }_{i+1}\right)=V_i\left({\varphi }_i\right){\text{e}}^{-\bar{V}\left({\varphi }_i\right)\frac{P_i}{RT}\left({\varphi }_{i+1}-{\varphi }_i\right)}.$ (26)

Once the PC volumes are at certain discrete values of evaporation due to exposure, the volatility factor, the evaporated fraction and the vapor pressure of each PC are easily calculated by the following ratios.

● Volatility factor:

$Z\left(\varphi \right)={\displaystyle {\sum }_{i=0}^{ncp}\frac{P_i{V}_i\left(\varphi \right)}{RTN\left(\varphi \right)}}$ . (27)

● Evaporative fraction:

$F_z\left(\varphi \right)=\frac{{\displaystyle {\sum }_{i=0}^{ncp}{V}_i}}{V_0}.$ (28)

● Vapor pressure is the sum of the partial pressure of the individual components such that

$P_T=\frac{{\displaystyle {\sum }_{i=0}^{ncp}{n}_i{P}_i}}{{\displaystyle {\sum }_{i=0}^{ncp}{n}_i}}.$ (29)

The last one is calculated from the amount of hydrocarbon that evaporates over time, using the previously created exposure evaporation so

$\frac{\text{d}V}{V}\text{d}t=-\left\{\frac{K^0\left(1-Y\left(t\right)\right)U^{\frac{7}{9}}}{h\left(t\right)}\right\}Z\left(\varphi \right).$ (30)

The integration time was carried out as follows:

$V\left(t_{n+1}\right)=V\left(t_n\right){\text{e}}^{\left[-\left\{\frac{K^0\left(1-Y\left(t_n\right)\right)U^{\frac{J}{9}}}{h\left(t_n\right)}\right\}Z\left(\varphi \left(t_n\right)\right)\left(t_{n+1}-t_n\right)\right]}.$ (31)

2) Vertical dispersion

Natural dispersion of oil in a spill is understood as the process of fragmentation of oil into droplets due to the breaking of the waves on the slick. According to Delvigne and Sweeney [22], and used recently by Li et al. [23], it is possible to estimate the number, size and distribution of droplets each time they break. The entry of oil into the water column is given by

$Q=C_{disp}{D}_e^{0.57}{f}_{bw}{V}_{disp},$ (32)

where $C_{disp}$ is an experimentally determined parameter, $f_{bw}$ is the fraction of breaking wave per wave period, and $V_{disp}$ is the volume of oil introduced into the column per unit volume of water.

The energy dissipation of the wave per unit area is defined by

$D_e=0.0034{\rho }_{w}g{H}_{rms}^2$ , (33)

in which $H_{rms}^2$ is the mean root of the wave height in meters, assuming the significant spectral height of wave $h_0$ as $H_{rms}=0.707H_0$ .

$V_{disp}$ is obtained by integrating the product of the droplet volume and the frequency distribution of the droplets on the oil volume. In practice, the integration is carried out between the minimum and maximum size that the droplet could reach, considering the experimental analyses.

$V_{disp}\propto {\displaystyle {\int }_{d_{\min }}^{d_{\max }}N\left(\delta \right){\delta }^3\text{d}\delta }$ (34)

Larger droplets, whose imposed maximum allowable size is 70 microns, will refloat in less time than necessary for the stain to pass through the area covered by the dispersed oil, thus reintegrating into the oil slick (forming the tail). Droplets 70 microns or less in size remain in suspension for a time determined by the buoyancy and depth of penetration into the water column. The minimum droplet size is generally considered to be 5 microns, from which the droplets are considered to be precipitated to the seabed [22].

$N\left(\delta \right)$ is the number of droplets of oil per unit volume of water per droplet size measured by

$N\left(\delta \right)=N_0{\left(\frac{{\delta }_0}{\delta }\right)}^{\frac{2}{3}}$ (35)

where $N_0$ and ${\delta }_0$ are experimental reference values.

3) Emulsification

Wave energy spontaneously mixes hydrocarbon with seawater, forming a “water-in-oil” emulsion. In a series of laboratory experiments with crude oils [24] it was found that the rate of emulsion formation is best described by a law of first-order change in the interfacial area. Quiao et al. [25] put it into practice in a series of modern experiment to obtain asphaltenes. The formulation to describe how water content and droplet size distribution modify the process with respect to time is

$\frac{\text{d}S}{\text{d}t}=k_s\left(1-\frac{S}{S_{\max }}\right),$ (36)

where the interfacial parameter $k_s$ is sensitive to wave energy; S and $S_{\max }$ are the oil-water interfacial area and the maximum interfacial area, respectively. The water fraction contained in the oil slick is related to the interfacial area and to the mean diameter of the water droplets $d_w$ such that

$Y=\frac{S{d}_w}{6+S{d}_w}$ . (37)

4) Change of physical-chemical properties

Evaporation and emulsification increase the density and viscosity of the slick, even when the oil is freshly spilled. Most hydrocarbons and their by-products are more viscous than water. Emulsified hydrocarbons may have a viscosity of several orders of magnitude greater than the corresponding fresh hydrocarbon and often behave like a non-Newtonian fluid. Density will also increase due to the weather, but not to the point of losing the buoyancy of the oil.

5) Density calculation

The two circumstances in which slick sinking has been recorded are when the slick mixes with the mineral sediment in the water column or when it burns in-situ, generating a residue that may be of high density. The weathering algorithms used will not identify any hydrocarbon that is denser than the surrounding water, as these will simply not be modeled.

A common way to present the density of a liquid is to compare it with the density of fresh water, giving a specific gravity of a fluid in relation to the water. However, in the oil industry, the common system of units is the API. The two references are related by the formula:

$G_e=\frac{141.5}{131.5+\text{API}}$ . (38)

Three factors affecting the slick density are considered: water temperature, degree of emulsification, and evaporated hydrocarbon fraction. The formula used that estimates this parameter in comparison with the density of fresh oil is

$\rho =Y{\rho }_w+\left(1-Y\right){\rho }_{ref}\left[1-c_1\left(T-T_{ref}\right)\right]\left(1+c_2{f}_{evap}\right),$ (39)

where $c_1$ and $c_2$ are empirical constants, $\rho$ is the hydrocarbon density, ${\rho }_w$ the water density, ${\rho }_{ref}$ the reference hydrocarbon density, Y the emulsion fraction, T the water temperature, $T_{ref}$ the reference temperature and $f_{evap}$ the evaporated hydrocarbon fraction.

6) Calculation of viscosity

Viscosity is the property that determines how freely a fluid flows and may be physically regarded as the internal friction of the fluid. Viscosity is significantly affected by temperature and increases with the time the oil has been outdoors. For a Newtonian fluid it is defined by the relationship between stress and shear rate. For a dimensional flow in the X direction, the $\tau$ stress is proportional to the velocity gradient.

$\tau =-\eta \frac{\text{d}{v}_x}{\text{d}y}$ (40)

Taking in account the mentioned temperature sensibility of viscosity, Andrade correlation is generally valid for small temperature variations [26]

$\ln \frac{{\eta }_2}{{\eta }_1}=C_T\left(\frac{1}{T_2}-\frac{1}{T_1}\right)$ (41)

The proportionality constant $C_T$ is taken with the value 5000 K.

As oil evaporates, the relative concentrations of the various components change. Perry’s rule [27] shows that the viscosity of a non-polar mixture may be derived from the viscosities of the components.

$\ln {\eta }_{\max }={\displaystyle {\sum }_{i=0}^{ncp}{w}_i\ln {\eta }_i}$ (42)

The following relationship given by Mackay et al. [18] is consistent with this mixing rule:

$\ln \frac{\eta }{{\eta }_0}=C_B{F}_z$ (43)

The empirical proportionality constant $C_B=5$ is used for this research.

Mooney’s equation [28] is widely used to describe the viscosity of emulsions.

$\frac{\eta }{{\eta }_0}={\text{e}}^{\left(\frac{2.5Y}{1-C_{Y}Y}\right)}$ (44)

Nguyen [29] uses it efficiently at the east of Vietnam Sea. This has been modified so that the relative viscosity may be adjusted for the size of the water droplets in the oil:

$C_{\eta }\left(t\right)=\frac{d_{\min }}{Y_{\max }d\left(t\right)}\left(1-\frac{2.5Y}{\ln R_{\eta }}\right)$ (45)

where $d\left(t\right)$ is the diameter of droplets based on time, $d_{\min }$ the minimum diameter of droplets, $F_z\left(t\right)$ the fraction of evaporated oil based on time, $R_n$ the maximum emulsion radius, T the temperature (K), v the fluid velocity, $Y\left(t\right)$ the fraction of water content based on time, $Y_{\max }$ the fraction of maximum water content, $n_{\max }$ the maximum dynamic viscosity and $n\left(t\right)$ the dynamic viscosity based on time.

7) Module implementation

The libraries “matplotlib”, “basemap” and “cmocean” were used, mainly in the preparation of the new graphic outputs.

2.2. Data-Providing Models. Bathymetry and Work Zones

Weather data used in the creation of the scenarios come from the following atmospheric models.

● WRF: Weather Research and Forecasting [30], implemented in the Center for Atmospheric Physics (CFA), with a spatial resolution of 18 km and the second domain with 6 km. Both domains are centered in 22 degrees north latitude and 80 degrees west longitude, and were used in both 30 vertical levels.

● MM5: Mesoscale Model 5 [31], implemented in the Atmosphere Physical Centre (CFA in Spanish) as a research tool. The work team used the domain that has a temporal resolution of 3 hours and spatial resolution of 27 km.

● NAM-CARIBE: North American Mesoscale-Caribbean Forecast System, downloaded from Internet. The domain considers the entire Greater Caribbean and presents a spatial resolution of 12 km and a time of 3 hours [32]. The download site is: http://nomads.ncep.noaa.gov/pub/data.

Among the models and projects that provide the variables of marine circulation, those mentioned below were chosen.

● ROMS: Regional Oceanic Modeling System [33] implemented in the CFA. It is performed in a domain with 3 km of spatial resolution, with 12 vertical levels and centered in Cuba.

● HYCOM: Hybrid Coordinate Ocean Model, downloaded from Internet. The horizontal grid has an average space between nodes of 7 km and 32 vertical layers [34]. The website http://www.hycom.org/dataserver/ provides access to the output of the global ocean prediction system in near real time.

● NCOM: Navy Coastal Ocean Model, downloaded from the Internet, is the U.S. Navy’s worldwide operating ocean model that includes several domains within which the Gulf of Mexico and the Caribbean Sea [35] appear (identified by “ncom_relo_amseas”). The download website used is: http://nomads.ncep.noaa.gov/pub/data.

● CW: COSTWATCH, downloaded from the Internet) this project does not yields predictive scenarios, just diagnostic. However, its values have been corrected for their drifter service, improving the quality of the data [36]. Access to the data was obtained from the following website: http://www.aoml.noaa.gov/.

In order to consider the parameters of the waves, the outputs of the following models were incorporated.

● WW III: WAVEWATCH III [37], implemented in the CFA, operates to a depth of 7.5 m, so it is a deep-sea wave generation model.

● SWAN: Simulating Waves Nearshores [38], run in the CFA, has a spatial resolution of 1 Km and best describes the wave transformation.

The bathymetric information for the Gulf of Mexico, Caribbean Sea and Near Atlantic, referring to bathymetric and coastline data, was obtained by integrating two sources. In first place The GEBCO project [39] and Nautical charts available at the Centre for Marine Meteorology (CMM) in the second one.

Characteristics of the basemap were conformed creating the regions and other components in the following way: scale 1:250,000, WGS84 projection, latitude range from 7.20 degrees to 33 degrees latitude north and longitude range from 61 degrees to 101 degrees longitude west. The model has twelve pre-conceived zones related to the Interamerican Seas, spatially designed by finite elements [2]. Either way, the user may design a new zone anywhere in the world and be perfectly incorporated into the system.

3. Results and Discussion

Simple ways to illustrate the graphical outputs results is one of the advantages of PETROMAR-3D model. The purpose is to give as accurate as possible forecasts to facilitate decision-makers quick and effective analyses and responses.

3.1. Description of the Blowout

Images describing the upward movement are illustrated in Figure 2 below. The first of these depicts the evolution of a plume in a three-dimensional projection and its latitudinal and longitudinal deviation in transit to the surface. Similarly, the concentration of oil is considered by means of different tonalities in the circles, which correspond to the values presented in the color palette. The others describe the vertical velocity and density deficiency as a function of time, and the radius of the plume as it ascends.

(a) (b) (c) (d)

Figure 2. Some model outputs for E = 0.01 of the CV: (a) Tree-dimensional tracking of a CV; (b) Vertical velocity during upward movement; (c) Density deficiency over time and (d) Relationship between radius of the CV and depth.

Outputs were obtained taking the value of the constant E equal to 0.01, as the best fit for the northwestern area of Cuban territorial waters.

Horizontal displacement is one of the most important variables for those responsible for taking action in the event of disasters. Its variations for different values of parameter E (Figure 3) are appreciable. A similar dispersion of the results of the zonal and southern displacements may be observed, while the percentage variations with respect to the variations of E are between 5% and 9%, greater in the case of the southern component of the movement. This may be explained by the weak southern component of the currents in the study area, which results in the inherent errors of the model and those of truncation having a greater relative weight in the results.

In dotted lines are shown the results of the ASCENSOMAR model, proposed by Pérez [40]. The curves obtained are substantially different, expected result since ASCENSOMAR is a simpler model that does not take into account the exchanges of mass and energy with the marine environment; however, for the proposed case study the total displacement is within the range of results obtained by the new model implemented. In particular, for the zonal displacement a maximum and minimum deviation of 35% and 1.80% respectively was calculated. In the southern displacements, these were 76% and 7.33% respectively.

Additionally, the new curves reflect, to a better extent, the behavior of the particle in its ascent, mainly the fact of having the concavity downwards (the slope to the curve decreases with height) which is theoretically more accepted. The displacement in the near spill is very limited because the exit speed is imposed and then gradually reduced to a limit speed.

3.2. Drift and Weathering Modules

Once the calculations of each weathering process are attributed to the virtual

(a)(b)

Figure 3. Comparison of the horizontal displacement for different values of E. (a) Zonal displacement; (b) Southern displacement.

particles of the Lagrangian method, the oil concentration values will be obtained at all times considered in the scenario as a sample of the existing intermodular relationship.

Numerical and graphical outputs are yielded every one hour giving information of state and position of all particles conforming the modelled oil slick which are been transported. The first ones are contained in a file with the variables of each run: initial time, final time and time of greatest arrival of the pollutant to a coastal sector, as well as location of virtual particles with their macroscale variables every one hour. This kind of output becomes very convenient because it can be used as input for other analysis and post-processing tools like geographic information systems to obtain, for example, the affected coastal sector in kilometers. In addition, the model can suggest the type of optimal response technique to be used over time to mitigate the environmental impact of the spill.

Graphical outputs are released to illustrate not only all the numerical releases but also other useful information like cartographical maps and meshes configured on PETROMAR-3D.

Figures below are results of a 5000 bb spill from an average oil of 30 API. They are very important for users because of their integration of useful information in a simple way. This feature allows them to take quick actions with respect to the incident.

Figure 4 illustrates four model outputs for first 72 hours after the oil spill. This scenario was taken on September 10, 2004 at 01 UTC. Initial coordinates are represented by an orange star.

Figure 5 illustrates the overlap of the oil slick drift at every time. In Figure 6 is shown a superposition of the wind and marine currents fields that are also incorporated. It is one of the benefits presented by the software so user visually notices the influences of each field.

Figure 4. Oil concentration forecast for (a) 1 hour, (b) 24 hours, (c) 48 hours and (d) 72 hours of the spill.

Figure 5. Graphical output with the final drift of the slick.

Figure 6. Graphical output with the slick and the fields of current and wind.

4. Conclusions and Recommendations

The version 2.1 of PETROMAR-3D model has a module that simulates a blowout in the sea depths. The enhancement of this process includes the characteristics of the type of oil, marine environment and the spill itself. With the introduction of weathering, the spill is better represented and the rest of the processes are simulated in a more realistic way. With these conditions it is possible to represent the concentration of the slick based on time, favoring a better analysis of the scenario.

Given the choice of the mentioned theory to describe the physical processes, the mathematical methods and the numerical schemes, the proposed model is considered to have advantages to be used in the simulation of spills in the inter-American seas. So the novelty of the work is directed to three fundamental aspects. First, the design of the region’s areas has been refined with great geographic details and the outputs of the main existing operating models have been considered. The drift, weathering and blowout processes have been unified in the same model, achieving similarity with others of advanced performance. Finally, the graphical and numerical outputs have been adjusted to the demands of decision-makers present in the region and there is still a list of suggestions that authors will consider.

It is recommended to include the influence of waves on surface drift through Stokes’s law [41], and convenient to incorporate a new procedure when the vertical speed becomes zero below sea surface, during the blowout. To prepare another article to evaluate the numerical precision and system optimization of the model, as well as present a proposal to give users different kinds of ways to access PETROMAR-3D software, it is also very appropriate.

Annex

– Configuration file of PETROMAR-3D

### AUTHOR: CENTRE FOR MARINE METEOROLOGY ###

General:

model_data: /home/petromar/PetroData

output_data: /home/petromar/PetrOut

#############################

### SCENARIO CREATION ###

### INPUT VARIABLES ###

#############################

Scenario:

name_scenario: Macondo_z16_ROMS_Py

initial_date: 22/04/2010

number_days: 30

### DESIGNED ZONES ###

###Golfo_zon; zona01; zona02; zona03; zona04; zona05; zona06; zona07; zona08

###zona09; golfo_zona01; golfo_zona01_06; zona01_02_03; zona01_02_03_06_11; zona01_06; zona03_04; zona04_05_06

###zona07_08; zona08_09

zona_s: zona01_06

### MODELS TO USE ###

### MARINE CURRENT MODELS: Curr ###

### hycom1; hycom2; coastwatch; roms; ncom

curr: roms

### ATMOSPHERE WIND MODELS: wd ###

### mm5; gfs; wrf

wd: 'wrf'

### WAVE MODELS: wv ###

### ww3; swan

wv: swan

### TYPE OF RUNS (tr): FORWARD (1) or BACKWARD (−1)

tr: 1

##############################

### TRAJECTORY CALCULATION ###

### INPUT VARIABLES ###

##############################

Model:

name_scenario: Macondo_z16_ROMS_PY

name_track: t1_rk4

### Runge-Kutta numerical schemes ###

### rk1; rk2; rk4

numeric_schema: rk4

### longitude(geog. degrees): x

### latitude(geog. degrees): y

x: −83.49

y: 23.31

### Kind of the spill (ks): “instantaneous” or “continuous”

ks: instantaneous

initial_time: 10

### Volume magnitude (mag_vol): “m3”, “liters”, “barrels”

mag_vol: barrels

Volume: 500

### Density (Ro): “density”, “specific gravity: ge”,

### “American Petroleum Institute: api”

Ro: api

Value: 20

Oilvis (cSt): 259.0

Oilvisreftem (C): 25.0

#Sea water density (kg/m³): Roa

Roa: 1027

#Water cinematic viscosity (m²/s): visc

visc: 140.267

#Superficial tension coeficient: coef_tension

coef_tension: 27

###############################

### Subsurface spills ###

###############################

salinity_p: 27

temperature_p: 35

depth_p: 100.0

density_p: 893

density15_p: 892

radio_v: 0.4016

### Release speed of blowout ###

v_p(m/s): 1

u_p(m/s): 1

w_p(m/s): 2.1

final_time: 700

constante_E: 0.5

######################## PETROMAR-3D END #######################

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Liu, Y., MacFadyen, A., Ji, Z.G. and Weisberg, R.H. (2013) Monitoring and Modeling the Deepwater Horizon Oil Spill: A Record Breaking Enterprise. John Wiley & Sons, Hoboken, 195.
[2]	Calzada, A.E., Acuna, F.E.P., Perdomo, D.R. and Taylor, R.C. (2015) Modelación de los derrames de petróleo mediante el empleo de PETROMAR. Revista Cubana de Meteorología, 21, 57-69. http://rcm.insmet.cu/index.php/rcm/article/view/413/473
[3]	Spaulding, M.L. (2017) State of the Art Review and Future Directions in Oil Spill Modeling. Marine Pollution Bulletin, 115, 7-19. https://doi.org/10.1016/j.marpolbul.2017.01.001
[4]	Arkhipov, B.V. and Shapochkin, D.A. (2019) Modeling Oil Spills in the Sea. Mathematical Models and Computer Simulations, 11, 107-120. https://doi.org/10.1134/S2070048219010046
[5]	Fay, J.A. (1971) Physical Processes in the Spread of Oil on a Water Surface. International Oil Spill Conference, 1, 463-467. https://doi.org/10.7901/2169-3358-1971-1-463
[6]	Zheng, L., Yapa, P.D. and Chen, F. (2003) A Model for Simulating Deepwater Oil and Gas Blowouts Part I: Theory and Model Formulation. Journal of Hydraulic Research, 41, 339-351. https://doi.org/10.1080/00221680309499980
[7]	De Dominicis, M., Pinardi, N., Zodiatis, G. and Lardner, R. (2013) MEDSLIK-II, a Lagrangian Marine Surface Oil Spill Model for Short-Term Forecasting. Part 1: Theory. Geoscientific Model Development, 6, 1851-1869. https://doi.org/10.5194/gmd-6-1851-2013
[8]	Ribotti, A., Antognarelli, F., Cucco, A., Falcieri, M.F., Fazioli, L., Ferrarin, C., Olita, A., Oliva, G., Pes, A. and Quattrocchi, G. (2019) An Operational Marine Oil Spill Forecasting Tool for the Management of Emergencies in the Italian Seas. Journal of Marine Science and Engineering, 7, 1. https://doi.org/10.3390/jmse7010001
[9]	Golshan, R., Boufadel, M.C., Rodriguez, V.A., Geng, X., Gao, F., King, Th., Robinson, B. and Tejada-Martínez, A.E. (2018) Oil Droplet Transport under Non-Breaking Waves: An Eulerian RANS Approach Combined with a Lagrangian Particle Dispersion Model. Journal of Marine Sciences and Engineering, 6, 1-16. https://doi.org/10.3390/jmse6010007
[10]	Lee, J.H. and Cheung, V. (1990) Generalized Lagrangian Model for Buoyant Jets in Current. Journal of Environmental Engineering, 116, 1085-1106. https://doi.org/10.1061/(ASCE)0733-9372(1990)116:6(1085)
[11]	Fernandez, F. (2018) Seawater: Seawater Library for Python. https://github.com/pyoceans/python-seawater
[12]	Fofonoff, N.P. and Millard, R.C. (1983) Algorithms for the Computation of Fundamental Properties of Seawater, Ocean Best Practices. https://www.oceanbestpractices.net/handle/11329/109
[13]	Hunter, J.D. (2018) User’s Guide—Matplotlib 3.1.1 Documentation. https://matplotlib.org/3.1.1/users/index.html
[14]	Betancourt, F.O. (2005) Modelado de la evolución de derrames de petróleo en zonas litorales. PhD Dissertation, Universidad Nacional Autónoma de México, México.
[15]	Gutiérrez, A.R. (2012) Dispersión y caoticidad de partículas pasivas en las aguas oceánicas de la región occidental de Cuba. PhD Dissertation, Instituto de Oceanología, Cuba. https://www.oceandocs.org/handle/1834/4847
[16]	Gros, J., Dissanayake, A.L., Daniels, M.M., Barker, C.H., Lehr, W. and Socolofsky, S.A. (2018) Oil Spill Modeling in Deep Waters: Estimation of Pseudo-Component Properties for Cubic Equations of State from Distillation Data. Marine Pollution Bulletin, 137, 627-637. https://doi.org/10.1016/j.marpolbul.2018.10.047
[17]	Kotzakoulakis, K. and George, S.C. (2018) Predicting the Weathering of Fuel and Oil Spills: A Diffusion-Limited Evaporation Model. Chemosphere, 190, 442-453. https://doi.org/10.1016/j.chemosphere.2017.09.142
[18]	Spanoudaki, K., Kampanis, N., Kalogerakis, N., Zodiatis, G. and Kozyrakis, G. (2018) Modelling of Oil Spills from Deep Sea Releases. EGU General Assembly Conference Abstracts, Vienna, 8-13 April 2018, 16401.
[19]	Mackay, D. and Matsugu, R.S. (1973) Evaporation Rates of Liquid Hydrocarbon Spills on Land and Water. The Canadian Journal of Chemical Engineering, 51, 434-439. https://doi.org/10.1002/cjce.5450510407
[20]	Mackay, D., Shiu, W.Y., Hossain, K., Stiver, W., and McCurdy, D. (1982) Development and Calibration of an Oil Spill Behavior Model. Toronto Univ. (Ontario) Dept. of Chemical Engineering and Applied Chemistry, Toronto. https://www.osti.gov/biblio/5464992
[21]	Bi, Y.C., et al. (2015) A Vapor-Liquid Equilibrium Model for Petroleum Based on Two-Parameter Antoine Equation. Computers and Applied Chemistry, 32, 324-328. http://en.cnki.com.cn/Article_en/CJFDTotal-JSYH201503014.htm
[22]	Delvigne, G.A.L. and Sweeney, C. (1988) Natural Dispersion of Oil. Oil and Chemical Pollution, 4, 281-310. https://doi.org/10.1016/S0269-8579(88)80003-0
[23]	Li, Z., Spaulding, M.L. and French-McCay, D. (2017) An Algorithm for Modeling Entrainment and Naturally and Chemically Dispersed Oil Droplet Size Distribution under Surface Breaking Wave Conditions. Marine Pollution Bulletin, 119, 145-152. https://doi.org/10.1016/j.marpolbul.2017.03.048
[24]	Eley, D., Hey, M. and Symonds, J. (1988) Emulsions of Water in Asphaltene-Containing Oils 1. Droplet Size Distribution and Emulsification Rates. Colloids and Surfaces, 32, 87-101. https://doi.org/10.1016/0166-6622(88)80006-4
[25]	Qiao, P., Harbottle, D., Tchoukov, P., Wang, X. and Xu, Z. (2017) Asphaltene Subfractions Responsible for Stabilizing Water-in-Crude Oil Emulsions. Part 3. Effect of Solvent Aromaticity. Energy Fuels, 31, 9179-9187. https://doi.org/10.1021/acs.energyfuels.7b01387
[26]	Green, D.W. and Perry, R.H. (2008) Perry’s Chemical Engineers’ Handbook. Eighth Edition, McGraw-Hill, New York.
[27]	Perry, R.L., Leung , P.T., Zwissler, C., Bouchard, R., Hervey, R., Fiorentino, L., Sharma, N., Martin, K., Howden, S., Kirkpatrick, B. and Kim, H.S. (2016) Collaborative Observational and Operational Oceanography to Improve 2015 and 2016 Loop Current Forecasting in the Gulf of Mexico. OCEANS 2016 MTS/IEEE Monterey, Monterey, 19-23 September 2016, 1-8. https://doi.org/10.1109/OCEANS.2016.7761275 https://ieeexplore.ieee.org/document/7761275
[28]	Mooney, M. (1951) The Viscosity of a Concentrated Suspension of Spherical Particles. Journal of Colloid Science, 6, 162-170. https://doi.org/10.1016/0095-8522(51)90036-0
[29]	Nguyen, T. (2017) Research and Development for Oil Spill Simulation Backward in Time at East Vietnam Sea. Journal of Petroleum & Environmental Biotechnology, 8, 344. https://doi.org/10.4172/2157-7463.1000344
[30]	Skamarock, W.C., Klemp, J.B., Dudhia, J., Gill, D.O., Barker, D.M., Wang, W. and Powers, J.G. (2008) A Description of the Advanced Research WRF Version 3. NCAR Technical Note-475+ STR.
[31]	Zou, X., Huan, W. and Xiao, Q. (1998) A User’s Guide to the MM5 Adjoint Modeling System. Mesoscale and Microscale Meteorology Division, National Center for Atmospheric Research, Boulder. https://opensky.ucar.edu/islandora/object/technotes%3A343/datastream/PDF/view
[32]	Nam (2017) North America Mesoescale Forecasting System. National Center for Environmental Information (NCEI) Formerly Known as National Climatic Data Center (NCDC). https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/north-american-mesoscale-forecast-system-nam
[33]	Moore, A.M., Arango, H.G., Broquet, G., Powell, B.S., Weaver, A.T. and Zavala-Garay, J. (2011) The Regional Ocean Modeling System (ROMS) 4-Dimensional Variational Data Assimilation Systems: Part I System Overview and Formulation. Progress in Oceanography, 91, 34-49. https://doi.org/10.1016/j.pocean.2011.05.004
[34]	Chassignet, E.P., Hurlburt, H.E., Smedstad, O.M., Halliwell, G.R., Hogan, P.J., Wallcraft, A.J., Baraille, R. and Bleck, R. (2007) The HYCOM (HYbrid Coordinate Ocean Model) Data Assimilative System. Journal of Marine Systems, 65, 60-83. https://doi.org/10.1016/j.jmarsys.2005.09.016
[35]	Martin, P.J., Barron, Ch.N., Smedstad, L.F., Campbell, T.J., Wallcraft, A.J., Rhodes, R.C., Rowley, C., Townsend, T.L. and Carroll, S.N. (2009) Navy Coastal Ocean Model (NCOM) Version 4.0 (User’s Manual). Naval Research Lab Stennis Space Center Ms Oceanography Div. https://doi.org/10.21236/ADA494652 https://www.ncdc.noaa.gov/data-access/model-data/model-datasets/navoceano-ncom-reg
[36]	Lumkin, R., Grodsky, S.A., Centurioni, L., Rio, M.-H., Carton, J.A. and Lee, D. (2013) Removing Spurios Low-Frecquency Variability in Drifter Velocities. Journal of Atmospheric and Oceanic Technology, 30, 353-360. https://doi.org/10.1175/JTECH-D-12-00139.1
[37]	Tolman, H.L. (2002) User Manual and System Documentation of WAVEWATCH III Version 2.22. National Oceanic and Atmospheric Administration, Technical Note 222. https://polar.ncep.noaa.gov/mmab/papers/tn276/MMAB_276.pdf
[38]	The SWAN Team (2020) SWAN User Manual. SWAN Cycle III Version 41.31A. Delft University of Technology, Delft, 143 p. http://swanmodel.sourceforge.net/online_doc/swanuse/swanuse.html
[39]	GEBCO (2017) Gridded Bathymetry Data (General Bathymetric Chart of the Oceans). British Oceanographic Data Centre, Liverpool. https://www.bodc.ac.uk/data/hosted_data_systems/gebco_gridded_bathymetry_data
[40]	Pérez, F. (2013) Simulación del ascenso del petróleo derramado en el lecho marino. VII Congreso de Meteorología SOMETCUBA, La Habana.
[41]	Robert, L. and Wettre, C. (2012) Comparison of Operational Oil Spill Trajectory Forecast with Surface Drifter Trajectories in the Barens Sea. Journal Geology and Geosciences, 1, 1-8.