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

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


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

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.

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 A. E. Calzada et al. Open Journal of Marine Science 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.
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.  In practice, it is convenient to define the term e Q as a result of a constant of proportionality E multiplied by an environmental function.
This constant is empirical and must be adjusted through experiments. The In the above Equation (3), a V is the average vector of the flow velocity over the exposed plume surface while a ′ V is a projection of a V 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 p ∆ is presented as the difference p a p p ∆ = − and D C 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. 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.

Drift Module
The drift processes are governed by the equations and considering the motion of a particle i in each axis. The a u and a v are the components of the speed of the particle, consequence of the advection and in the other hand, d u and d v are those referring to diffusion. Thus, putting equations (6) and (7) in perspective, you would obtain that the advection part could be modelled like 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 In these formulations, h D represents the horizontal diffusion coefficient, choosing 10 5 m 2 /s for the calculations [15]. R 1 and R 2 are independent random variables within the [ ] 0,1 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 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.

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 d d First five PCs of equal volume are created and their boiling points are determined by Vapor pressure of each PC is based on the Antoine equation [21] such that ( ) 2 2,1 0 2, 2, The relative molar volume i V 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).
Reference mass transfer coefficient In addition, the exposure evaporation matrix is constructed consisting of four columns: exposure evaporation, volatility factor Z, evaporated fraction 2 F and vapor pressure.
Volatility factor is a measure of how volatile the PC is and is defined by 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.
The incremental change in evaporation due to exposure also increases with the index of each PC as follows.
( ) The differential equations describing evaporation are then solved in these preset values of evaporation due to exposure.
A numerical approximation in terms of lnV is used: 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: • Evaporative fraction: The last one is calculated from the amount of hydrocarbon that evaporates over time, using the previously created exposure evaporation so The integration time was carried out as follows: 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 where disp C is an experimentally determined parameter, bw f is the fraction of breaking wave per wave period, and disp V 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 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].
where 0 N and 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 where the interfacial parameter s k is sensitive to wave energy; S and max S 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 w d such that

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.
where 1 c and 2 c are empirical constants, ρ is the hydrocarbon density, w ρ the water density, ref ρ the reference hydrocarbon density, Y the emulsion fraction, T the water temperature, ref T the reference temperature and evap f 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 τ stress is proportional to the velocity gradient.

ln
The proportionality constant T C 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.
The following relationship given by Mackay et al. [18] is consistent with this mixing rule: The empirical proportionality constant 5 B C = is used for this research.
Mooney's equation [28] is widely used to describe the viscosity of emulsions.
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: where ( )

7) Module implementation
The libraries "matplotlib", "basemap" and "cmocean" were used, mainly in the preparation of the new graphic outputs.

Data-Providing Models. Bathymetry and Work Zones
Weather data used in the creation of the scenarios come from the following atmospheric models. 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 Open Journal of Marine Science 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.

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.

Description of the Blowout
Images describing the upward movement are illustrated in Figure 2 below. The 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.

Drift and Weathering Modules
Once the calculations of each weathering process are attributed to the virtual Open Journal of Marine Science  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.

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