Modelling of hydrocarbons dispersion. Adaptation and improvement of a state of the art model

FINAL DISSERTATION Title Modelling of hydrocarbons dispersion. Adaptation and improvement of a state of the art model Author Paula Martín Laguna Tutor Agustín Sánchez-Arcilla Conejo and Marc Mestres Ridge Department Ingeniería marítima Intensification Modelización Date September 2016 Coastal waters quality is a big concern in our society. It is necessary to assure the parameters of the big ecosystem which is the sea are as invariable as possible so the life in it progresses and does not decrease. For this purpose, it is necessary to own tools in case a spill happens so the authorities can quickly design an action plan. There are several models which simulate the dispersion of hydrocarbons in the sea but few of them do it following both Lagrange and Stokes dispersion and movement theories. MEDSLIK_II is one of them but until now it was only able to process one certain type of input data. This dissertation brings a modification of the software so it is able a widely used type of currents and sea temperature forecast data: ROMS. After understanding the theory behind MEDSLIK_II, a new subroutine has been coded so the software becomes more universal. Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 2 Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 3 Index Index ....................................................................................................................................................... 3 Figures .................................................................................................................................................... 5 Tables ...................................................................................................................................................... 6 1. Introduction ................................................................................................................................... 7 2. Objectives and method ................................................................................................................. 8 3. MEDSLIK_II. Theory background ............................................................................................. 9 3.1. Introduction ............................................................................................................................. 9 3.2. Model equations and state variables ..................................................................................... 10 3.3. MEDSLIK-II tracer grid and solution methodology ............................................................. 14 3.4. Time rate of change of slick state variables .......................................................................... 16 3.5. Time rate of change of particle oil volume state variables ................................................... 18 3.6. Time rate of change of particle positions .............................................................................. 18 3.7. Numerical considerations ...................................................................................................... 22 4. MEDSLIK II. Code Architecture .............................................................................................. 26 4.1. Source Code .......................................................................................................................... 26 4.2. Meteo-oceanographic files .................................................................................................... 27 4.3. Input Data files ...................................................................................................................... 27 4.4. Script files and executables ................................................................................................... 28 4.5. Visualization software........................................................................................................... 28 4.6. Output Data files ................................................................................................................... 28 5. MEDSLIK_II. Input files analysis ............................................................................................. 30 5.1. Introduction ........................................................................................................................... 30 5.2. MEDFF files ......................................................................................................................... 30 5.3. ROMS files ........................................................................................................................... 31 5.4. Comparison between inputs .................................................................................................. 34 6. MEDSLIK_II. Modifications ..................................................................................................... 35 6.1. Introduction ........................................................................................................................... 35 6.2. Extract_II.for ......................................................................................................................... 35 6.3. Bathymetry and coast line data ............................................................................................. 36 6.4. medslik_II.for ........................................................................................................................ 37 7. MEDSLIK_II. Case study .......................................................................................................... 38 8. Conclusions and future modifications ....................................................................................... 40 References ............................................................................................................................................. 41 Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 4 Appendix ............................................................................................................................................... 43 Input parameters in medslik_inputfile.txt ......................................................................................... 44 Input parameters in medslik5.par ...................................................................................................... 47 Subroutine added to Extract_II.for .................................................................................................... 49 Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 5 Figures Figure 1: The advective and diffusive processes which are involved in the transport of contaminants in shallow seas (Noye, 1987) ................................................................................................................................................... 9 Figure 2: Graphic representation of concentration classes (De Dominicis, et al., 2013) .................................... 11 Figure 3: MEDSLIK-II model solution procedure methodology (De Dominicis, et al., 2013) ............................ 15 Figure 4: Weathering processes using Mackay's approach: evaporation (E), dispersion (D) and spreading (S) (De Dominicis, et al., 2013).................................................................................................................................. 17 Figure 5: Ekman transport graphic ...................................................................................................................... 20 Figure 6: MEDSLIK-II Reference Code Architecture (De Dominicis, 2012) ....................................................... 26 Figure 7: Data storage of Medff files and reading direction ............................................................................... 31 Figure 8: Point of measurement of each variable ................................................................................................ 31 Figure 9: ROMS horizontal grid (Anon., 2016) ................................................................................................... 32 Figure 10: ROMS horizontal extracted grid (Anon., 2016) .................................................................................. 32 Figure 11: Example of ROMS vertical grid (Anon., 2015) ................................................................................... 33 Figure 12: Data storage of ROMS files and reading direction ............................................................................ 33 Figure 13: Structure of a file extracted from the main ROMS file ....................................................................... 36 Figure 14: Region for the case study simulation .................................................................................................. 38 Figure 15: Partial view of the files generated by Extract_II.for .......................................................................... 38 Figure 16: Structure of the hourly output file ....................................................................................................... 39 Figure 17: Plot of Algeria test .............................................................................................................................. 39 Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 6 Tables Table 1: Oil spill model state variables (De Dominicis, et al., 2013) .................................................................. 12 Table 2: Meteo-oceanographic model outputs initially run by MEDSLIK-II ....................................................... 27 Table 3: ROMS variables' size ............................................................................................................................. 32 Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 7 1. Introduction The quality of the sea is currently a big concern to our society. Climate change together with uncontrolled spills from factories and cities’ sewages threat everyday our coastal waters. Moreover, it is also important to take into account the oil pipes and overseas oil platforms, which need to follow high standards to ensure they do not compromise the quality of the water. However, although governments are keeping a close eye on it and are continuously developing new strategies and policies, it is still necessary to own tools which help in case of a disaster, such as accidental leakages or spills. Numerical modelling is a powerful tool in these cases. It allows us to work with a huge amount of data in a, relatively, small amount of time. The simulations performed by these models are an important source of information for the actions plans. There are currently several models which have been develop in order to predict how a leak will disperse across the surface and the water column. Special attention is given to hydrocarbons leakages since they are a huge pollution source. The theory behind the simulation is usually either Lagrange, Stokes or a combination of both of them. The main inputs in every case are, on the one hand, sea variables such as currents speed, sea surface temperature and bathymetry; and, on the other, information about the leak such as type of oil and initial location. The sea data comes from different programs which extract it from boys located in the sea and process it. There are many sea data sources around the world but, unluckily, they do not generated files with the same structure. This generates problems when trying to standardize the dispersion of hydrocarbons modelling. This dissertation has its initial point in an already existing free software called MEDSLIK II, developed originally by Robin Lardner and George Zodiatis and later modified by Michaela De Dominicis. This program can simulate oil dispersion across the entire Mediterranean Sea, allowing the user to reduce the sea section to work with by dividing it in different areas. MEDSLIK II combines both Lagrange and Stokes theories when simulating the movement of the oil particles. This document contains a chapter which goes through the theory behind the modelling process. Even though MEDSLIK II is prepared to model any part of the Mediterranean Sea, it cannot be used in certain areas because the sea data locally obtained, ROMS data, does not match the data MEDSLIK II is prepared for. As said before, it is necessary and positive to own the more standardized tools possible so they can be used any time anywhere in the world. Therefore, this dissertation aims to upgrade MEDSLIK II so it can also process the sea data coming from ROMS, which is a free-surface, terrain-following, primitive equations ocean model widely used by the scientific community for a diverse range of applications. Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 14 concentration. Also, it allows defining a solution methodology, since there is no analytical relationship between the oil slick and the particle state variables. 3.3. MEDSLIK-II tracer grid and solution methodology To connect Equation 2 and Equation 3, a discrete oil tracer grid system is defined, xT=(xT, yT), with a uniform but different grid spacing in the zonal and meridional directions, (δxT, δyT). Therefore, the spatially discretized time evolution equations for the structural state variables are 𝑑𝐶𝑆 𝑑𝑡(𝑥𝑇,𝑦𝑇,𝑡)=𝜌 𝛿𝑥𝑇𝛿𝑦𝑇𝑑𝑉𝑆 𝑑𝑡(𝑥𝑇,𝑦𝑇,𝑡) (14) 𝑑𝐶𝐷 𝑑𝑡 (𝑥𝑇,𝑦𝑇,𝑡)=𝜌 𝛿𝑥𝑇𝛿𝑦𝑇𝑑𝑉𝐷 𝑑𝑡 (𝑥𝑇,𝑦𝑇,𝑡) (15) The coastline is designed as a polygonal chain of points connected by segments of different lengths δLi, which finally constitute the coastline segment Li. Following this, Equation 6 becomes 𝐶𝐶(𝐿𝑖,𝑡)=𝜌 𝛿𝐿𝑖·𝑉𝐶(𝐿𝑖,𝑡) (16) By relating the particle state variables to the oil tracer grid, it is possible to write the relationship between structural and particle state variables. That means being able to, for instance, follow the evolution of the different concentrations. The countable groups, IS, ID, of surface and subsurface particles contained in an oil tracer grid cell are defined as 𝐼𝑆(𝑥𝑇,𝑦𝑇,𝑡)= { 𝑛𝑘;𝑥𝑇−𝛿𝑥𝑇 2=𝑥𝑘(𝑡)=𝑥𝑇+𝛿𝑥𝑇 2 𝑦𝑇−𝛿𝑦𝑇 2=𝑦𝑘(𝑡)=𝑦𝑇+𝛿𝑦𝑇 2 𝜎(𝑛𝑘,𝑡)=0 (17) 𝐼𝐷(𝑥𝑇,𝑦𝑇,𝑡)= { 𝑛𝑘;𝑥𝑇−𝛿𝑥𝑇 2=𝑥𝑘(𝑡)=𝑥𝑇+𝛿𝑥𝑇 2 𝑦𝑇−𝛿𝑦𝑇 2=𝑦𝑘(𝑡)=𝑦𝑇+𝛿𝑦𝑇 2 𝜎(𝑛𝑘,𝑡)=1 (18) Taking into account that, the surface and dispersed concentration can be reconstructed as {𝐶𝑆(𝑥𝑇,𝑦𝑇,𝑡)=𝜌 𝛿𝑥𝑇𝛿𝑦𝑇∑𝜌(𝑛𝑘,𝑡) 𝑛𝑘𝜖𝐼𝑆 𝐶𝐷(𝑥𝑇,𝑦𝑇,𝑡)=𝜌 𝛿𝑥𝑇𝛿𝑦𝑇∑𝜌(𝑛𝑘,𝑡) 𝑛𝑘𝜖𝐼𝐷 (19) To calculate the oil concentration on the coast, the set of particles attached to the coastal segment Li, IC, is used: 𝐼𝐶(𝐿𝑖,𝑡)={𝑛𝑘;𝜎(𝑛𝑘,𝑡)=−𝐿𝑖} (20) Then, the concentration of oil on each coastal segment is calculated by Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 15 𝐶𝐷(𝐿𝑖,𝑡)=𝜌 𝛿𝐿𝑖∑𝜌(𝑛𝑘,𝑡) 𝑛𝑘𝜖𝐼𝐷 (21) In order to solve the different concentrations using the oil slick and particle state variable equation, it is necessary to create a sequential solution. As it is represented in Figure 3 MEDSLIK_II sets the initial conditions and solves the transformation processes (evaporation, dispersion, spreading). After, the thin and thick parts are updated, together with the particles volumes. After, the calculations to relocate the particle positions and update the particle status index are performed. Finally, MEDSLIK_II calculates the oil concentrations as described in Equations 19 and 21. Figure 3: MEDSLIK-II model solution procedure methodology (De Dominicis, et al., 2013) One of the most significant approximations is that the oil slick state variables depend only on the slick’s centre geographical position, which is updated after each advection-diffusion time step. The oil spill centre position is defined as follows: 𝑥𝐶(𝑡)=∑𝑥𝑘(𝑡) 𝑁 𝑘=1𝑁 (22) 𝑦𝐶(𝑡)=∑𝑦𝑘(𝑡) 𝑁 𝑘=1𝑁 (23) 3.3.1. Initial conditions The oil release can be instantaneous or continuous. When the leakage lasts several hours or even months, it can happen that the initial volumes spilled have been transported away from the release site when the last ones are released (Liu, et al., 2011). To model that, Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 16 MEDSLIK_II divides the total spill in sub-spills, consisting of a part of the oil released during an interval of time: 𝑁𝑆=𝐷𝐶 𝑇𝐶 (24) where DC(s) is the release duration. For instantaneous releases, the initial oil released is equal to the total oil released VS(xC, t0). However, for a continuous oil spill release, for each interval the volume released is 𝑉𝑆(𝑥𝐶,𝑡0)=𝑅𝐶𝑇𝐶 (25) where Rc is the oil spill rate. During an instantaneous release, N particles are released at the beginning, while in a continuous one, NC particles are released every TC: 𝑁𝐶=𝑁 𝑁𝑆 (26) The initial particle volume can be written as 𝜕(𝑛𝑘,𝑡0)=𝑁𝑆𝑉𝑆(𝑥𝐶,𝑡0) 𝑁 (27) And also the evaporative and non-evaporative oil volume components: 𝜕𝐸(𝑛𝑘,𝑡0)=(1−𝜕𝑁𝐸 100)𝜕(𝑛𝑘,𝑡0) (28) 𝜕𝑁𝐸(𝑛𝑘,𝑡0)=𝜕𝑁𝐸 100𝜕(𝑛𝑘,𝑡0) (29) Where 𝜕𝑁𝐸 is the percentage of the non-evaporative component of the oil that depends on the oil type. The first thin and thick area values are taken from the initial surface amount of oil released using relative thicknesses and F, which is the area ratio. These three values are defined as input. 𝐴𝑇𝑁(𝑡0)=𝐹𝐴𝑇𝐾 (30) 𝐴𝑇𝐾(𝑡0)=𝑉𝑆(𝑥𝐶,𝑡𝑜) 𝑇𝑇𝐾(𝑥𝐶,𝑡0)+𝐹𝑇𝑇𝑁(𝑥𝐶,𝑡0) (31) 3.4. Time rate of change of slick state variables Following Equation 11, the time rate of change of oil volume is written as 𝜕𝑉𝑆 𝜕𝑡 =𝜕𝑉𝑇𝐾 𝜕𝑡 +𝜕𝑉𝑇𝑁 𝜕𝑡 (32) Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 17 Figure 4 shows the three main processes that induce changes of the surface oil volume. They are known as weathering processes: evaporation, dispersion and spreading. The first one acting on the oil slick is evaporation, since initially the main volume is on the surface. In addition, for the first hours, the spill spreads mechanically due to gravitational forces. Figure 4: Weathering processes using Mackay's approach: evaporation (E), dispersion (D) and spreading (S) (De Dominicis, et al., 2013) The three processes are considered separately for the thick and thin slick and only at the slick centre. Therefore, the forecast equations are 𝑑𝑉𝑇𝐾 𝑑𝑡 =𝑑𝑉𝑇𝐾 𝑑𝑡 (𝐸)+𝑑𝑉𝑇𝐾 𝑑𝑡 (𝐷)+𝑑𝑉𝑇𝐾 𝑑𝑡 (𝑆) (33) 𝑑𝑉𝑇𝑁 𝑑𝑡 =𝑑𝑉𝑇𝑁 𝑑𝑡 (𝐸)+𝑑𝑉𝑇𝑁 𝑑𝑡 (𝐷)+𝑑𝑉𝑇𝑁 𝑑𝑡 (𝑆) (34) Each term of these equations is described in detail in Appendices B1, B2 and B4 in De Dominicis, et al (2013) paper, given in terms of modified Mackay fate algorithms for evaporation, dispersion and spreading. Following Mackay’s assumptions, TTN does not change. Hence, ATN is calculated as 𝑑𝐴𝑇𝑁 𝑑𝑡 =1 𝑇𝑇𝑁𝑑𝑉𝑇𝑁 𝑑𝑡 (35) Where On the other hand, for the thick slick 𝑑𝑉𝑇𝐾 𝑑𝑡 =𝑇𝑇𝐾𝑑𝐴𝑇𝐾 𝑑𝑡 +𝐴𝑇𝐾𝑑𝑇𝑇𝐾 𝑑𝑡 (36) The area of the thick slick only changes due to spreading, so 𝑑𝐴𝑇𝐾 𝑑𝑡 =𝑑𝐴𝑇𝐾 𝑑𝑡 (𝑆) (37) VTK is updated using Equation 34 and the thickness changes are calculated by 𝑇𝑇𝐾=𝑉𝑇𝐾 𝐴𝑇𝐾 (38) Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 18 3.5. Time rate of change of particle oil volume state variables After the transformation processes have acted on the oil slick, the particle oil volumes are changed. For all particle status index, the evaporative oil particle volume changes following the empirical relationship 𝜕𝐸(𝑛𝑘,𝑡)=[1−𝜕𝑁𝐸 100−𝑓(𝐸)(𝑥𝐶,𝑡)]𝜕(𝑛𝑘,𝑡) (39) Where f(E) is the fraction of oil evaporated defined as: 𝑓(𝐸)(𝑥𝐶,𝑡)=𝑉𝑇𝐾(𝑥𝐶,𝑡)(𝐸)+𝑉𝑇𝑁(𝑥𝐶,𝑡)(𝐸) 𝑉𝑇𝐾(𝑡0)+𝑉𝑇𝑁(𝑡0) (40) For both surface and dispersed particles, the non-evaporative oil component does not change, while a certain fraction of the non-evaporative oil component of a beached particle can be modified due to adsorption processes occurring on a particular coastal segment, seeping into the sand or forming a layer on rocky shore. For beached particles, the non-evaporative oil component is then 𝜕𝑁𝐸(𝑛𝑘,𝑡)=𝜕𝑁𝐸(𝑛𝑘,𝑡0)0.5𝑡−𝑡0 𝑇𝑆(𝐿𝑖)𝜎(𝑛𝑘,𝑡)=−𝑖 (41) where 𝑡0 is the instant when the particle passes from surface to beached status and vice versa, 𝑇𝑆(𝐿𝑖) is a half-life for seepage or any other mode of permanent attachment to the coasts. It describes the “absorbency” of the shoreline by describing the rate of entrainment of the oil after it has landed at a given shoreline (Shen, et al., 1987). The half-life depends on the coastal type. 3.6. Time rate of change of particle positions Immediately upon entering into a water body, the spilled oil spreads and forms a surface slick which covers a large area of the water surface and can be moved about by the action of winds, waves or currents (Shen, et al., 1987). The evolution of the particle positions is given by nk uncouple Langevin equations: 𝑑𝑥𝑘(𝑡) 𝑑𝑡 =𝐴(𝑥𝑘,𝑡)+𝐵(𝑥𝑘,𝑡)𝜕𝜉(𝑡) (42) where the tensor A represents the deterministic part of the flow field and B the stochastic term together with 𝜉 which is a random factor. If we define the Wiener process 𝑊(𝑡)=∫𝜕(𝜎)𝛿𝜎 𝑡 0 (43) Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 19 and apply the It𝑜 assumption (Tompson & Gelhar, 1990), Equation 42 becomes equivalent to the It𝑜 stochastic differential equation: 𝑑𝑥𝑘(𝑡)=𝐴(𝑥𝑘,𝑡)𝑑𝑡+𝐵(𝑥𝑘,𝑡)𝑑𝑊(𝑡) (44) where dW is a random increment. The Wiener process describes the path of a particle due to Brownian motion modelled by independent random increments dW(t) sampled from a normal distribution with zero mean and second order moment dt. Hence, dW(t) in Equation 44 can be replaced by a vector Z of independent random numbers, normally distributed, i.e. Z𝛻N(0,1), and multiplied by √𝑑𝑡: 𝑑𝑥𝑘(𝑡)=𝐴(𝑥𝑘,𝑡)𝑑𝑡+𝐵(𝑥𝑘,𝑡)𝑍√𝑑𝑡 (45) The unknown tensors 𝐴(𝑥𝑘,𝑡) and 𝐵(𝑥𝑘,𝑡) are written as (Risken, 1984): 𝑑𝑥𝑘(𝑡)=[𝑈(𝑥𝑘,𝑡) 𝑉(𝑥𝑘,𝑡) 𝑊(𝑥𝑘,𝑡)]𝑑𝑡+[√2𝐾𝑥0 0 0√2𝐾𝑦0 0 0 √2𝐾𝑧][𝑍1 𝑍2 𝑍3]√𝑑𝑡 (46) where K are the turbulent diffusivity coefficients. For particles at the surface and dispersed, Equation 46 takes the form 𝑑𝑥𝑘(𝑡)=[𝑈(𝑥𝑘,𝑦𝑘,𝑧𝑘,𝑡) 𝑉(𝑥𝑘,𝑦𝑘,𝑧𝑘,𝑡) 0]𝑑𝑡+[𝑑𝑥′𝑘(𝑡) 𝑑𝑦′𝑘(𝑡) 𝑑𝑧′𝑘(𝑡)] (47) The particle position does not change for particles at the surface. It can only change when the particle becomes dispersed and the horizontal velocity at the vertical position of the particle is used to displace the dispersed particles. The Equation 47 can be expanded in different components according to the particle status index: {𝜎=0𝑑 𝑥𝑘(𝑡)=[𝑈𝐶(𝑥𝑘,𝑦𝑘,0,𝑡)+𝑈𝑊(𝑥𝑘,𝑦𝑘,0,𝑡)+𝑈𝑆(𝑥𝑘,𝑦𝑘,0,𝑡)]𝑑𝑡+𝑑𝑥′𝑘 𝜎=1 𝑑𝑥𝑘(𝑡)=𝑈𝐶(𝑥𝑘,𝑦𝑘,0,𝑡)𝑑𝑡+𝑑𝑥′𝑘 (48) where Uc is the Eulerian current velocity due to non-local wind, Uw is due to local wind and Us is the velocity due to wave-induced currents or Stokes drift. 3.6.1. Current and local wind velocity terms Ocean currents near the surface are generated by the atmospheric forcing. Although it can be subdivided into buoyancy fluxes and wind stresses, the last one is the most important in terms of kinetic energy of the induced motion (Wunsch, 1998). The wind-induced currents are attributable to non-local winds, and are dominated by geostrophic or quasi-geostrophic Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 20 dynamic balances (Pedlosky, 1986), where the pressure gradient force is balanced by the Coriolis effect and creates a flow parallel to the isobars. By definition, geostrophic and quasigeostrophic motions have a timescale of several days and characterize oceanic mesoscale motion, a very important component of the large-scale flow field included in U. These currents dominate below the mixed layer, which its dynamics are typically considered ageostrophic. The dominant time-dependent, wind-induced currents in the surface layer are the Ekman currents due to local winds (Price, et al., 1987). All these components have to be considered in the UC field in Equation 48. Ekman transport is described in Figure 5; surface currents flow at a 45º angle to the wind due to a balance between the Coriolis force and the drags generated by the wind and the water (Mann, 2006). If the ocean is divided vertically into thin layers, the main value of the velocity decreases until it dissipated. The direction also shifts slightly across each subsequent layer (right in the northern hemisphere, left in the southern hemisphere). This is called Ekman spiral (Knauss, 1978). If all flow over the length of the spiral is integrated, the net transportation is at 90º to the right (left) of the surface wind in the northern (southern) hemisphere (Colling, 2001). Figure 5: Ekman transport graphic Therefore, Ekman currents at the surface can be parametrized as a function of wind intensity and angle between winds and currents, i.e. 𝑈𝑊=𝛼[𝑊𝑥𝑐𝑜𝑠𝛽+𝑊𝑦𝑠𝑖𝑛𝛽] (49) 𝑉𝑊=𝛼[−𝑊𝑥𝑠𝑖𝑛𝛽+𝑊𝑦𝑐𝑜𝑠𝛽] (50) where Wx and Wy are the wind zonal and meridional components and α and β are two parameters referred as drift factor and drift angle. Current velocity fields can be obtained from analyses and forecasts produced by highresolution ocean general circulation models (OGCMs). Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 21 3.6.2. Wave current term Waves push pollutants by wave-induced velocities that are known as Stokes drift velocity, US. It has to be added to the currents from OGCMs since most of the ocean models are not coupled with wave models (Röhrs, et al., 2012). Stokes drift is the net displacement of a particle in a fluid due to wave motion, resulting essentially from the fact that the particle moves faster forward when the particle is at the top of the wave circular orbit than it does backward when it is at the bottom of its orbit. Wave-driven transport is often ignored, although it may be the dominant mechanism transporting oil to adjacent beaches and coastal wetlands, whose environment has a big susceptibility to oil spills. It provides a natural mechanism for beaching of surface oil, one of the most damaging outcomes of a coastal oil spill (Sobey & Barker, 1997). Considering the surface, the Stokes drift velocity intensity in the direction of the wave propagation is 𝐷𝑆(𝑧=0)=2∫𝜕𝐾(𝜕)𝑆(𝜕)𝑑𝜕 ∞ 0 (51) where w is angular frequency and S(w) is wave spectrum. To implement Equation 51, it is considered that the direction of the wave propagation is the same as the wind direction. Therefore, 𝑈𝑆=𝐷𝑆𝑐𝑜𝑠𝜗 (52) 𝑉𝑆=𝐷𝑆𝑠𝑖𝑛𝜗 (53) Where 𝜗 is the wind direction. 3.6.3. Turbulent diffusivity terms It is assumed that the particle moving through the fluid gets a random impulse at each time step, due to the action of incoherent turbulent motions and that it has no memory of its previous turbulent displacement: 𝑑𝑥′𝑘(𝑡)=(2𝑟−1)𝛿 (54) Where d is the particle mean path and r is a random real number. The mean square displacement of Equation 54 is |𝑑𝑥′𝑘(𝑡)2|=∫[(2𝑟−1)𝛿]2𝑑𝑟 1 0=13𝛿2 (55) while the mean square displacement of the turbulent terms in Equation 47 is 𝑑𝑥′𝑘(𝑡)2=2𝐾𝑑𝑡. Equalizing both terms the following is obtained: 𝛿2=6𝐾𝑥𝑑𝑡 Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 22 𝛿2=6𝐾𝑦𝑑𝑡 𝛿2=6𝐾𝑧𝑑𝑡 (56) Replacing these terms in Equation 54, it is possible to write the stochastic transport terms: 𝑑𝑥′𝑘(𝑡)=𝑍1√2𝐾𝑥𝑑𝑡=(2𝑟−1)√6𝐾ℎ𝑑𝑡 𝑑𝑦′𝑘(𝑡)=𝑍2√2𝐾𝑦𝑑𝑡=(2𝑟−1)√6𝐾ℎ𝑑𝑡 𝑑𝑧′𝑘(𝑡)=𝑍3√2𝐾𝑧𝑑𝑡=(2𝑟−1)√6𝐾𝑣𝑑𝑡 (57) where Kh and Kv are prescribed turbulent horizontal and vertical diffusivities. 3.7. Numerical considerations In this section, the interpolation method between input fields and the oil tracer grid, to the numerical scheme used to solve Equations 33, 34 and 47, to the model time step and to the oil tracer grid selection is explained. 3.7.1. Interpolation method The environmental variables which are used in this model (atmospheric wind, ocean currents and sea surface temperature) are normally supplied on a different numerical grid than the oil slick centre or particle locations. Therefore, to calculate the advection process, it is necessary to compute the currents and winds at the particle locations. However, for the transformation processes calculation, the sea surface temperature and winds are interpolated at the slick centre. To show these interpolation processes, (𝑥𝐸,𝑦𝐸,𝑧𝐸) is designed as the numerical grid on which the environmental variables are provided by the Eulerian models. A pre-processing procedure is needed to reconstruct the currents in the area between the last water grid node of the oceanographic model and the real coastline. MEDSLIK_II extrapolates the currents in a way that gives a velocity field to the land points: 𝑞𝑥𝐸(𝑖),𝑦𝐸(𝑖)=𝑞𝑥𝐸(𝑖+1),𝑦𝐸(𝑖)+𝑞𝑥𝐸(𝑖−1),𝑦𝐸(𝑖)+𝑞𝑥𝐸(𝑖),𝑦𝐸(𝑖−1)+𝑞𝑥𝐸(𝑖),𝑦𝐸(𝑖+1) 𝑁℘ (58) Then, the winds and currents are computed at the particle position (𝑥𝑘,𝑦𝑘) for a fixed depth 𝑧𝐸 with the following interpolation algorithm: 𝑞1=𝑞𝑥𝐸(𝑖),𝑦𝐸(𝑖)[𝑥𝐸(𝑖+1)−𝑥𝑘] 𝑞2=𝑞𝑥𝐸(𝑖+1),𝑦𝐸(𝑖)[𝑥𝑘−𝑥𝐸(𝑖)] 𝑞3=𝑞𝑥𝐸(𝑖),𝑦𝐸(𝑖+1)[𝑥𝐸(𝑖+1)−𝑥𝑘] Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 23 𝑞4=𝑞𝑥𝐸(𝑖+1),𝑦𝐸(𝑖+1)[𝑥𝑘−𝑥𝐸(𝑖)] 𝑞𝑥𝑘,𝑦𝑘=(𝑞1+𝑞2)[𝑦𝐸(𝑖+1)−𝑦𝑘]+(𝑞3+𝑞4)[𝑦𝑘−𝑦𝐸(𝑖)] ∆𝑥𝐸∆𝑦𝐸 (59) where (xk,yk) is the particle position referenced to the oil tracer grid, (xe(i),ye(i)) and subsequent are the external four grid points. The same algorithm is used to interpolate to the oil slick centre, (𝑥𝐶(𝑡),𝑦𝐶(𝑡) the wind and sea surface temperature. Besides, a vertical interpolation of the currents at the particle positions is also needed: 𝑞𝑥𝑘,𝑦𝑘,𝑧𝑘=1 𝑧𝐸(𝑖)−𝑧𝐸(𝑖+1){𝑞𝑥𝑘,𝑦𝑘,𝑧𝐸(𝑖+1)[𝑧𝐸(𝑖)−𝑧𝑘]+𝑞𝑥𝑘,𝑦𝑘,𝑧𝐸(𝑖)[𝑧𝑘−𝑧𝐸(𝑖+1)]} (60) 3.7.2. Numerical time integration scheme The Lagrangian horizontal particle motion, Equation 41, is solved using a Euler forward scheme as follows 𝑥𝑘(𝑡+∆𝑡)=𝑥𝑘(𝑡)+𝑈(𝑥𝑘,𝑡)∆𝑡+∆𝑥′𝑘(𝑡) (61) where xk(t) represents the particle position at the current time step. Equations 33 and 34 are also solved the same way as before but with a different time step, socalled weathering time step, 𝑉𝑇𝐾(𝑡+𝛿𝑡)=𝑉𝑇𝐾(𝑡)+𝑑𝑉𝑇𝐾 𝑑𝑡 𝛿𝑡 (62) 𝑉𝑇𝑁(𝑡+𝛿𝑡)=𝑉𝑇𝑁(𝑡)+𝑑𝑉𝑇𝑁 𝑑𝑡 𝛿𝑡 (63) The two processes above described are not solved with the same speed: transformation processes are faster than advection-diffusion processes. Since transformation equations are stiff, to integrate them it is necessary that the time step is shorter, a fraction, of the Lagrangian time step (Butenschön, et al., 2012). 3.7.3. Particle status updates After the particles have moved for a Lagrangian time step, their oil volume and status are updated. A particle becomes dispersed if the following probability function 𝑃(𝐷)(𝑡)=𝑓(𝐷)(𝑥𝐶,𝑡)−𝑓(𝐷)(𝑥𝐶,𝑡−∆𝑡) (1−𝑓(𝐷)(𝑥𝐶,𝑡−∆𝑡)) (64) Is greater than a random number defined between 0 and 1: 𝑟<𝑃(𝐷)(𝑡)=𝛻𝜎(𝑛𝑘,𝑡)=1 (65) Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 30 5. MEDSLIK_II. Input files analysis 5.1. Introduction All the information related to the currents, winds, temperature and bathymetry is stored in netCDF format files. The Network Common Data Form (netCDF) software was developed by Unidata, a National Science Foundation-sponsored program empowering U.S. universities which has extended its influence around the world. It works as an I/O library, which can be called by C, FORTRAN, C++ and other languages. The library stores and retrieves data in self-describing, machineindependent datasets in the form of arrays. An array is an n-dimensional rectangular structure containing terms which all have the same data type. The type of data encompasses from single-point observations and time series to satellite or radar images. Array values are accessed directly, without knowing how data is stored. Extra information such as units is stored with the data (Rew, et al., 2016). The netCDF libraries allow simultaneous access to multiple netCDF datasets which are identified by ID numbers in addition to ordinary file names. Each dataset contains dimensions, variables and attributes. The dimensions have a name and a length, which is a positive integer. The variables have a name and an ID number. They represent an array of values of the same type. Each variable has a name, a data type, and associated attributes. It is important to identify these aspects so they can be properly called from the main program. MEDSLIK II is prepared to work with netCDF files in a certain way. Here relays the aim of this dissertation, to upgrade MEDSLIK II so it can read an extensively used nectCDF type of file: ROMS. The following section describe ROMS structure, although first, the original input files are reviewed so it is easier to identify conflict points. 5.2. MEDFF files The Medff system provides with one file for each variable: temperature, u-currents and vcurrents. One file contains hourly the information of one day. However, this system does not provide the data from 00:00 to 23:59; instead, it is given from 13:00 to 12:59 of the following day. Each variable is calculated in the middle of the cell, so the string where they are stored has the same size: imx x jmx x kmx x ktmx (longitude length, latitude length, depth levels, measurement hours). Important is to notice that the depth levels are measured from the surface down, levels that are constant at any (lon , lat) position. There is a variable which stores the depth values in meters. Figure 7 shows a diagram of how the data is stored. Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 31 Figure 7: Data storage of Medff files and reading direction Figure 8: Point of measurement of each variable The file also contains information relatively the longitude and the latitude of each cell centre, although the program does not use it because it works with already set up areas. 5.3. ROMS files ROMS is a free-surface, terrain-following, primitive equations ocean model widely used by the scientific community for a diverse range of applications. The data resolution is quite convenient, since the user can adjust both vertical and horizontal resolution. A unique ROMS file contains all three variables: temperature, u-currents and v-currents. Also, when extracting the file, the user can ask for the duration of the forecast. In other words, instead of having one file per day, there is only one file containing all the information. Since the user selects when downloading the data which period of time they want the information about, once they have the file, the only way of figuring out the date is throw an internal calculation based on the file’s variable scrum_time. This variable refers to the time since initialization since 1st January at 00:00. Special attention needs to be given to where ROMS calculates each variable. Scalar ones are at the cell’s centre (temperature) while vector ones are calculated at the contour (u-currents and v-currents). Figure 9 shows the grid and the location of each variable. Even though currents (both components) are defined in all the cells contour, when extracting data, one can usually just find the interior range of parameters, as it is shown in Figure 10. That results in not having any current data in the external contour. Furthermore, the dimensions of the three variables (temperature, u-currents and v-currents) are not the same, as it is shown in Table 3. Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 32 Figure 9: ROMS horizontal grid (Anon., 2016) Figure 10: ROMS horizontal extracted grid (Anon., 2016) Variable Size Temperature imx x jmx x kmx x ktmx u-currents imx-1 x jmx x kmx x ktmx v-currents imx x jmx-1 x kmx x ktmx Table 3: ROMS variables' size Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 33 Regarding depths, ROMS has its own reference system. It does not have, as Medff, horizontal levels at determined depths, which means that at shallow waters it stores null data at a certain point. Instead, it only takes into consideration the actual water levels by stretching curves over the ocean bottom. Figure 11 shows an example of these curves. Figure 11: Example of ROMS vertical grid (Anon., 2015) There is a simple equation which transforms from ROMS reference system to z-system in case it is needed to know to which z-depth a certain data belongs. In Section 6.2 this transformation is reviewed. ROMS data starts at the bottom and goes up to the ocean surface. Figure 12 shows a diagram of how the data is stored (the same type of diagram as in Medff files has been used, although the depth levels do not correspond since they are not horizontal). Figure 12: Data storage of ROMS files and reading direction Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 34 ROMS files contain many other variables regarding the currents and temperature. For this dissertation, special relevance is given to bathymetry at RHO-points variable and to the coefficients needed to transform S-coordinates to Z-coordinates. 5.4. Comparison between inputs Once the structure of both input types is clear, the main differences are identified. The target is to create a subroutine that extracts the data from the ROMS file and transforms it so it is equal to the one MEDSLIK II is programmed to work with. Overall, the main differences are:  Medff files store data of one day counting from 13:00 to 12:59 (including then data from two different calendar days). Therefore, there will be several files, one for each day. ROMS, however, contains all the information of every day in a unique file, starting at 00:00 to 23:59 of the last day of forecast.  MEDSLIK_II identifies the files by its file name which has the following structure: YYMMDDHH. ROMS files have a standard name, “ocean_his”. Therefore, when extracting the data, the files need to be renamed so medslik_II.for can read them.  The original MEDSLIK_II is programmed in order to get information from three different files, one for each current component and one for the temperature. However, ROMS includes these three variables in a unique file.  In the original program, one can find the bathymetry information already stored in a library in .bath files. However, the ROMS file includes the bathymetry of the area to model. Therefore, it is necessary to create the .bath file from the information in there. This will allow the user to perform the model in any area, not only in the ones predefined.  The three variables of interest, both components of the currents and temperature, are taken from the centre of the grid cell. However, ROMS presents the currents on the border of the cell and the temperature in the centre (see Figure 10). Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 35 6. MEDSLIK_II. Modifications 6.1. Introduction MEDSLIK_II was programmed to run at specific areas of the Mediterranean Sea. This dissertation brings a modification which makes MEDLIK_II universal; it can be used anywhere if the user has de required input files. In this section, the main additions and modifications to the code are listed and reviewed. One has to consider that minor changes are not mentioned, although in the appendix of this paper the entire code is attached with the modifications and additions marked. Most of the changes affect the source code, specially Extract_II.for. This code processes the input files (in this case only one) based on the data of the spill and extracts the information needed for the simulation. One of the targets of this dissertation is to create a clean and clear code so it can be read and modified easily by other users. Therefore, the new code follows the structure of the old one so it does not affect its harmony. 6.2. Extract_II.for One of the ROMS highlights is its data resolution. Since the area covered by the data is rather small (compared with original MEDSLIK_II data, which covers the entire Mediterranean), it is also the area where MEDSLIK_II performs the modelling. Therefore, it is not necessary, as it is with the other type of files, to select a subarea from the original area. This results in having a much bigger grid (in terms of cells number), but with better resolution. First step is to interpolate the u and v-currents strings so the values are located in the cells centre. Since only two values (from each border of the cell) are available, the interpolation is lineal Since temperature values are already stored at the centre, after the interpolation is done, the three variables available to the next step, which is the vertical interpolation. As it has been described before, ROMS data has a terrain-following vertical grid, so at shallow waters it does not store null numbers (as it would with a standard vertical grid). MEDSLIK_II works with the data at certain depths (10, 30 and 120 m) in all cells, so it is necessary to perform a transformation of the vertical grid to know the depth location of each cell: 𝑧(𝑥,𝑦,𝜎,𝑡)=𝜁(𝑥,𝑦,𝑡)+[𝜁(𝑥,𝑦,𝑡)+ℎ(𝑥,𝑦)]𝑆(𝑥,𝑦,𝜎) (75) 𝑆(𝑥,𝑦,𝜎)=ℎ𝑐𝜎+ℎ(𝑥,𝑦)𝐶(𝜎) ℎ𝑐+ℎ(𝑥,𝑦) (76) Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 36 Where S(x,y,σ) is a nonlinear vertical transformation functional, ζ(x,y,t) is the time-varying free-surface, h(x,y) is the unperturbed water column thickness (a.k.a. bathymetry), σ is a fractional vertical stretching coordinate ranging from -1 to 0, C(σ) is a nondimensional, monotonic, vertical stretching function ranging from -1 to 0 and hc is a positive thickness controlling the stretching. It is important to notice that, 𝑆(𝑥,𝑦,𝜎)={0 𝑖𝑓 𝜎=0, 𝐶(𝜎)=0 𝑎𝑡 𝑓𝑟𝑒𝑒 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 1 𝑖𝑓 𝜎=−1, 𝐶(𝜎)=−1 𝑎𝑡 𝑡ℎ𝑒 𝑜𝑐𝑒𝑎𝑛 𝑏𝑜𝑡𝑡𝑜𝑚 [78] As it is shown by the transformation above, ROMS files have a changing vertical grid throw time. This, again, gives even more precision to the modelling, since the data (u and v currents) is exactly measured at the exact depth. Once the depth in meters string is available (z), several loops are performed to define the value of both u and v string that is located at the desired depth. Finally, all the information which MEDSLIK_II requires to run is written in a file. The is a file for each time step. Figure 13: Structure of a file extracted from the main ROMS file 6.3. Bathymetry and coast line data A part from the u and v-currents and temperature, MEDSLIK_II also needs the bathymetry of the modelling area, its coast line and the type of coast for each coast line segment. MEDSLIK_II has a library where the user can find bathymetry, coast line and type of coast files for each area for which it was programmed to run for. Therefore, it is necessary that the user manually introduce these files for their case since ROMS files do not include this information (only the bathymetry which is extracted in Extract_II.for and stored in a file ready to be read by MEDSLIK_II). These files need to have the following structure: - Bathymetry (*.bath): Heading Min Longitude Max Longitude Min Latitude Max Latitude No. ‘x-cells’ Nº ‘y-cells’ Cells bathymetry starting from the higher left one and moving to the left. - Coast line (*.map) Total no. of contours (coast lines) Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 37 No of-points at the coast line 0 (island) or 1 (open coast) lon lat ... - Type of coast (*.cst): For each coordinate of a coast point, the user must definí a type of coast according to this code: 1 sand beach 2 sand and gravel beach 3 cobble beach 4 rocky shore 5 seawall; concrete, wharf, etc 6 exposed headland 7 sheltered sand or gravel beach 8 sheltered rocky shore 9 sheltered marsh or mud flats These three files have to be place at the MEDSLIK_II/EXE/data folder. For each simulation, these files need to be created accordingly. 6.4. medslik_II.for medslik_II.for contains the simulation code based on the theory reviewed in section X. Minor modifications need to be applied since one of the improvement targets is to upgrade the code following its original structure: the ROMS file is processed so the information obtain has the same structure as the one medslik_II.for is programmed for. Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 38 7. MEDSLIK_II. Case study To check the new code, a ROMS file containing data from Vigo has been used. The main effort has been generating through the routine Extract_II.for forecast data files which the same structure are the previously admitted by the program. A test case of the original program (which modelled a region of Algeria) has also been used to compare different aspects. In this case, the ROMS file contains information of a 176x153x10x169 grid. This means a 7days forecast of the region showed in Figure 14: Region for the case study simulation In this case, the grid will contain information regarding land areas. This increases unnecessarily the size of the files, resulting in a higher computational time if these points weren’t considered. However, the geography and size of the area make impossible to reduce the grid because the leakage can disperse from its original point (e.g. Vigo’s harbour located inside the ria) to any part of the current grid. The files are correctly extracted and stored in the forecast folder for the subsequent simulation (a total of 192 files). Figure 15: Partial view of the files generated by Extract_II.for Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 39 The bathymetry and coast line where already available. The coast line includes both Cíes islands although only part of them are considered in the simulation (see Figure 14). Finally, medslik_II.for processes all input data so files containing oil concentration are created for each hour after the spill. Figure 16: Structure of the hourly output file In this case, the graphic output is not valid due to the resolution. Further modifications need to be added to adapt ROMS resolution to the plot’s. However, in Figure 17 one of the plots obtained in the case study is shown. This plot is useful to know the direction of the spill’s movement. However, it is not really useful regarding numerical data. Figure 17: Plot of Algeria test Overall, except the plot section, the model is considered as validated. Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 46 for the observation is read from a file named initial.txt in the main model directory ($HOME/MEDSLIK_II_1.01/EXE). This file has been converted by the ReadSatData.py from a ‘.gml’ or ‘.xml’ file derived from the satellite image. Satellite data is received from EMSA in a file of type ‘*.gml’ which can be visualized on Google Earth. Such files may contain data for several oil slicks, with the boundary of each digitized to a fine resolution, especially in the case of a gml file. If SAT_DATA=YES is chosen, the following parameters have to be filled (otherwise leave them blank): namefileGML has to be filled with the name of the GML file (e.g. namefileGML=ASA_WSM_1PNACS20080806_095116_000000612071_00022_336 43_0001.N1.00000_Oil.gml). The file has to be saved in the directory $HOME/MEDSLIK_1.01/EXE. N_OS has to be filled with the number of the oil slick in the file to be simulated (e.g. N_OS=1) Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 47 Input parameters in medslik5.par The following will provide a brief description of the significance of the parameters in the medslik5.par. These are the parameters that a user is most likely to want to change from their default values. Stokes drift correction. Choosing 01 allows the model to use the wave ‐ induced velocity (Stokes drift) calculated using an empirical formulation. Choosing 00 Stokes drift velocity will not be applied. The default value is 01. Wind correction (Drift Factor). The drift speed of the slick is equal to this factor multiplied by the wind speed. The default value is 0. Wind correction (Drift Angle). The wind - ‐ driven drift of the slick occurs at this angle to the right of the wind direction. The default value is 0.0 degrees, which causes the slick to move directly downwind. Variable Drift Angle. Choosing 01 allows the model to use a drift angle that decreases as the wind speed increases; a wind speed at which the drift angle is reduced by 50% must then be entered in the line below. The default is that such a reduction is not made. Reduction of Forecast Wind Speed. When using forecast water circulation in a simulation, the forecast water velocities already include the effect of the wind forces on the water surface. It may thus be considered appropriate in some cases to reduce the wind speed used in the drift formula by a fraction of the winds used in the forecast. This can be done writing 01 in effective wind speed and entering in the line below the reduction fraction (between 0 and 1). The default is that such a reduction is not made. Smagorinsky Scheme. Choosing 01 allows the horizontal diffusivity to be computed from the water currents using the Smagorinsky scheme. The default value is 0. Horizontal Diffusivity. Enter the diffusivity that determines the horizontal diffusive spreading of the slick. The default value is 2.0 m2/s. A larger value will cause the slick to spread faster. Vertical Diffusivities. The model allows the use of two values of vertical diffusivity: a larger value in the top well - ‐ mixed layer and a smaller value below the mixed layer. The defaults are 0.01 and 0.0001 m2/s respectively. Depth of Mixed Layer. The default for this depth is 30 m. Number of parcels. Diffusion and dispersion are modelled using a Monte Carlo algorithm, representing the oil by a large number of particles which are then given appropriate random displacements, and the resulting cloud of parcels used to estimate the concentration of the oil in the slick. MEDSLIK - ‐ II uses 10,000 parcels as the default, but the user can increase this up to 300,000. Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 48 Depths of Forecast Currents. The current fields are extracted from the Nectdf and are given in the text files in $HOME/MEDSLIK_II_1.01/EXE/fcst_data at fixed depths. Whatever depths are used in the text files (they can only be changed by modifying the Extract_II.for code), the same depths must be entered in the parameters form. In the MEDSLIK - ‐ II_1.01 version the current fields are given at 3 depths (10 m, 30 m and 120 m) plus the surface. Selection of currents for convection of slick. In the MEDSLIK - ‐ II_1.01 version it is possible to choose between 4 options: 00 for surface currents, 01 for 10 m deep currents, 02 for 30 m deep currents and 03 for 120 m deep currents. The default strategy is to use the surface velocities for advection of the slick and to ignore the wind drift altogether (by setting the drift factor equal to zero). Number of Time Steps per Hour. The default value is 2. The default time step for computation of the convection and diffusion is 30 minutes. (A shorter step is used for the fate processes.) In general, this is adequate, but for a continuous spill with strong winds and/or currents, the 30 - ‐ minute step may cause the computed slick to appear as a number of discrete slicks that do not merge for several hours. This lack of reality can be reduced by using a shorter time step. Dimension of the Array used for representing the slick. Each time the results are printed in one of the output files, the particles are aggregated inside ‘pixels’, the size of which was entered in the medslik_inputfile.txt (grid_size). These pixels form an array the dimension of which is set by this parameter. The default array dimension is 2000×2000 and the maximum is 4000×4000. Choosing this maximum causes the run program to slow very significantly. On the other hand, too small a dimension causes the displayed slick to have straight, barrier - ‐ like, edges. If this occurs, you must either choose a larger array dimension or a larger pixel size. Modelling of hydrocarbons’ dispersion in the sea Paula Martín Laguna 49 Subroutine added to Extract_II.for Archivo: /home/paula/MEDSLIK_II_1.01/EXE/source/Extract_II.for Página 1 de 12 1c----------------------------------------------------------------------------------- 2c MEDSLIK-II_1.01 3c oil spill fate and transport model 4c----------------------------------------------------------------------------------- 5c Extract_II.for 6c This routine reads winds and currents from 7c meteo-oceanogrpahic model output (NetCDF files) 8c----------------------------------------------------------------------------------- 9c Copyright (C) <2012> 10 c This program was originally written 11 c by Robin Lardner and George Zodiatis. 12 c Subsequent additions and modifications 13 c have been made by Michela De Dominicis. 14 c This version contains additions and modifications 15 c by Paula Martín Laguna so the program is able to process 16 c as inputs ROMS files. 17 c---------------------------------------------------------------------------------- 18 c The development of the MEDSLIK-II model is supported by a formal agreement 19 c Memorandum of Agreement for the Operation and Continued Development of MEDSLIK-II 20 c signed by the following institutions: 21 c INGV - Istituto Nazionale di Geofisica e Vulcanologia 22 c OC-UCY - Oceanography Center at the University of Cyprus 23 c CNR-IAMC - Consiglio Nazionale delle Ricerche – Istituto per 24 c lo Studio dell’Ambiente Marino Costiero 25 c CMCC - Centro Euro-Mediterraneo sui Cambiamenti Climatici 26 c 27 c This program is free software: you can redistribute it and/or modify 28 c it under the terms of the GNU General Public License as published by 29 c the Free Software Foundation, either version 3 of the License, or 30 c any later version. 31 c 32 c This program is distributed in the hope that it will be useful, 33 c but WITHOUT ANY WARRANTY; without even the implied warranty of 34 c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 35 c GNU General Public License for more details. 36 c You should have received a copy of the GNU General Public License 37 c along with this program. If not, see < http://www.gnu.org/licenses/ >. 38 c----------------------------------------------------------------------------------- 39 40 character regn*4, indate(30)*8,indate_wind(30)*8, fc_dir*120, 41 & filename*9 42 common regn, alon1, alon2, alat1, alat2, numfiles, indate, 43 & numfiles_wind,indate_wind,iviod, icurrents,fc_dir, 44 & filename 45 integer len_dir 46 c-------------------------------------------------------------------- 47 c read subregion limits & file dates. Adjust limits to lie on OPA grid 48 c-------------------------------------------------------------------- 49 call getarg(1,fc_dir) 50 len_dir=120 51 do while(fc_dir(len_dir:len_dir).eq.' ') 52 len_dir=len_dir-1 53 enddo 54 55 56 open(1,file='medslik.tmp') 57 read(1,*) regn, icurrents, iwind 58 print*, icurrents 59 read(1,*) alon1,alon2 60 read(1,*) alat1,alat2 61 read(1,*) numfiles 62 if(icurrents.ne.100) then 63 do n=1,numfiles+1 64 read(1,'(a8)') indate(n) 65 enddo 66 read(1,*) numfiles_wind 67 do n=1,numfiles_wind 68 read(1,'(a8)') indate_wind(n) 69 enddo 70 else 71 read(1,'(a9)') filename 72 print*, filename Archivo: /home/paula/MEDSLIK_II_1.01/EXE/source/Extract_II.for Página 2 de 12 73 endif 74 read(1,*) iviod 75 close(1) 76 77 open(99,file='Extract.log') 78 if(icurrents.eq.10) call ExtractOPA(fc_dir,len_dir) 79 if(icurrents.eq.11) call ExtractADRI24(fc_dir,len_dir) 80 if(icurrents.eq.12) call ExtractSICI24(fc_dir,len_dir) 81 if(icurrents.eq.13) call ExtractTYRR24(fc_dir,len_dir) 82 83 if(icurrents.eq.70) call ExtractOPA_1hr(fc_dir,len_dir) 84 if(icurrents.eq.71) call ExtractSICI(fc_dir,len_dir) 85 if(icurrents.eq.72) call ExtractADRI(fc_dir,len_dir) 86 if(icurrents.eq.73) call ExtractTYRR(fc_dir,len_dir) 87 if(icurrents.eq.74) call ExtractRELO(fc_dir,len_dir) 88 if(icurrents.eq.75) call ExtractWESTMED(fc_dir,len_dir) 89 90 if(icurrents.eq.100) call ExtractROMS_1hr(fc_dir,len_dir, 91 & filename) 92 93 if(iwind.eq.6) call ExtractECMWF(fc_dir,len_dir) 94 if(iwind.eq.25) call ExtractECMWF25(fc_dir,len_dir) 95 stop 96 end 97 98 c********************************************************************** 99 c Extract medslik files from ROMS data 1hr 100 c********************************************************************** 101 subroutine ExtractROMS_1hr(fc_dir,len_dir,filename) 102 103 parameter(ktmx=169, imx=176, jmx=153, kmx=10, imx_U=175, 104 & jmx_V= 152, theta_s=6, thetha_b=0.9, hc=1) 105 106 real sc_r(10), cs_r(10) 107 real fmis !netcdf 108 parameter(fmis=0.) !netcdf 109 integer start(4), start_lonlat (2), start_mask(2), start_h(2), 110 & start_st(2), count(4), count_lonlat(2), count_mask(2), 111 & count_h(2), count_st(2) !netcdf 112 integer id, idlon, idlat, idU, idV, idT, idmask_rho, idh, idst !netcdf 113 integer Status !netcdf 114 integer msk(imx,jmx), mask(imx,jmx), st(ktmx,1) 115 real oplon(imx), oplat(jmx), 116 & ts(imx,jmx,kmx), u(imx,jmx,kmx),v(imx,jmx,kmx), 117 & ts_tmp(imx,jmx,kmx), u_tmp(imx,jmx,kmx), 118 & v_tmp(imx,jmx,kmx), lon_tmp(imx), lat_tmp(jmx), 119 & ts_169(imx,jmx,kmx,ktmx), u_169(imx,jmx,kmx,ktmx), 120 & v_169(imx,jmx,kmx,ktmx), lon(imx,jmx), lat(imx,jmx), 121 & u_in(imx_U,jmx,kmx,ktmx), v_in(imx,jmx_V,kmx,ktmx), 122 & z(imx,jmx,kmx), h(imx,jmx) 123 character Startd*6, filename*9, prdate*16, outfile*40,infile*120, 124 & heads*150, empty*80, regn*4, ora*2, ore*2, fc_dir*120, 125 & dia*2, hour*2, dayc*2, monthc*2 126 logical ex 127 integer t, o, kount, nore, len_dir, Startdi, Startdd, d, s, dur, 128 & month, day, dayS 129 common regn, alon1, alon2, alat1, alat2, numfiles, indate, 130 & iviod, icurrents 131 data udef /9999./, rhoa /1.19/ 132 133 c---------------------------------------------------------------------- 134 c Read ROMS file 135 c---------------------------------------------------------------------- 136 c Get to the directory 137 138 Status = 0 139 infile=fc_dir(1:len_dir)//'/fcst_data/ROMS/' 140 & //filename(1:9)//'.nc' 141 142 len_file=120 143 do while(infile(len_file:len_file).eq.' ') 144 len_file=len_file-1 Archivo: /home/paula/MEDSLIK_II_1.01/EXE/source/Extract_II.for Página 3 de 12 145 enddo 146 147 c Open netCDF file 148 149 Status = nf_open(infile(1:len_file),nf_nowrite,id) 150 print *,id 151 call handle_err(Status) 152 153 c---------------------------------------------------------------------- 154 c ROMS horizontal grid 155 c---------------------------------------------------------------------- 156 c Check of the variables from ROMS file 157 158 Status = nf_inq_varid (id, 'lon_rho', idlon) 159 call handle_err(Status) 160 Status = nf_inq_varid (id, 'lat_rho', idlat) 161 call handle_err(Status) 162 Status = nf_inq_varid (id, 'u', idU) 163 call handle_err(Status) 164 Status = nf_inq_varid (id, 'v', idV) 165 call handle_err(Status) 166 Status = nf_inq_varid (id, 'temp', idT) 167 call handle_err(Status) 168 Status = nf_inq_varid (id, 'scrum_time', idst) 169 call handle_err(Status) 170 Status = nf_inq_varid (id, 'mask_rho', idmask_rho) 171 call handle_err(Status) 172 Status = nf_inq_varid (id, 'h', idh) 173 call handle_err(Status) 174 175 print*,idh 176 c Since longitudes and latitudes at RHO points are avaliable, the grid 177 c is already built. It's just necessary to extract the strings holding 178 c the data. 179 180 c Longitudes string (at RHO points) 181 182 start_lonlat(1) = 1 183 start_lonlat(2) = 1 184 185 count_lonlat(1) = imx 186 count_lonlat(2) = 1 187 188 Status = nf_get_vara_real (id, idlon, start_lonlat, count_lonlat, 189 & lon_tmp) 190 call handle_err(Status) 191 oplon=lon_tmp 192 c Latitudes string (at RHO points) 193 194 start_lonlat(1) = 1 195 start_lonlat(2) = 1 196 197 count_lonlat(1) = 1 198 count_lonlat(2) = jmx 199 200 Status = nf_get_vara_real (id, idlat, start_lonlat, count_lonlat, 201 & lat_tmp) 202 call handle_err(Status) 203 oplat=lat_tmp 204 205 c The spill will affect the entire grid due to the size of the area. 206 c Therefore, it is easy to update the coordinates of the RHO points 207 c and number of total cells affected by the spill. 208 209 alon1=oplon(1) 210 alon2=oplon(imx) 211 alat1=oplat(1) 212 alat2=oplat(jmx) 213 214 i_first = 1 215 i_last = imx 216 j_first = 1 Archivo: /home/paula/MEDSLIK_II_1.01/EXE/source/Extract_II.for Página 4 de 12 217 j_last = jmx 218 219 imax = imx 220 jmax = jmx 221 222 c Write information in the file Extract.log 223 224 write(99,*) 'i-limits = ',i_first,i_last,imax 225 write(99,*) 'j-limits = ',j_first,j_last,jmax 226 write(99,*) 'lon-limits = ',alon1,alon2,imx 227 write(99,*) 'lat-limits = ',alat1,alat2,jmx 228 229 c---------------------------------------------------------------------- 230 c read ROMS data file 231 c---------------------------------------------------------------------- 232 c All the information regarding currents and temperature is stored in 233 c the same netcdf file, which has been already open to extract the 234 c longitudes and latitudes. In this section, a part from extracting 235 c the relevant strings, a interpolation is performed, since the data 236 c corresponding to the currents has been measured at the borders of the 237 c grid and MEDSLIK_II works with all the data at the centre of it. 238 c A linear interpolation is performed, taking into acount that the 239 c borders of the grid do not have information, therefore no 240 c interpolation is possible there. More details can be found at the 241 c report. 242 c---------------------------------------------------------------------- 243 c Extract U 244 245 start(1) = 1 246 start(2) = 1 247 start(3) = 1 248 start(4) = 1 249 250 count(1) = imx-1 251 count(2) = jmx 252 count(3) = kmx 253 count(4) = 1 254 255 do t = 1,ktmx 256 start(4)=t 257 Status = nf_get_vara_real (id, idU, start, count, u_tmp) 258 call handle_err(Status) 259 u_in(1:imx_U,1:jmx,1:kmx,t) = u_tmp(1:imx_U,1:jmx,1:kmx) 260 enddo 261 262 c Linear interpolation. Because of the data avaliable, 263 c first and last columns can't be interpolated 264 265 u_169(1,1:jmx,1:kmx,1:ktmx)=u_in(1,1:jmx,1:kmx,1:ktmx) 266 u_169(imx,1:jmx,1:kmx,1:ktmx)=u_in(imx_U,1:jmx,1:kmx,1:ktmx) 267 268 c Interpolation for the rest 269 270 do l=2,imx_U 271 do m=1,jmx 272 do n=1,kmx 273 do o=1,ktmx 274 u_169(l,m,n,o)=u_in(l-1,m,n,o)+ 275 & (u_in(l,m,n,o)-u_in(l-1,m,n,o))/2 276 enddo 277 enddo 278 enddo 279 enddo 280 281 c Extract V 282 283 start(1) = 1 284 start(2) = 1 285 start(3) = 1 286 start(4) = 1 287 288 count(1) = imx Archivo: /home/paula/MEDSLIK_II_1.01/EXE/source/Extract_II.for Página 5 de 12 289 count(2) = jmx-1 290 count(3) = kmx 291 count(4) = 1 292 293 do t = 1,ktmx 294 start(4)=t 295 Status = nf_get_vara_real (id, idV, start, count, v_tmp) 296 call handle_err(Status) 297 v_in(1:imx,1:jmx_V,1:kmx,t) = v_tmp(1:imx,1:jmx_V,1:kmx) 298 enddo 299 300 c Linear interpolation. Because of the data avaliable, 301 c first and last rows can't be interpolated 302 303 v_169(1:imx,1,1:kmx,1:ktmx)=v_in(1:imx,1,1:kmx,1:ktmx) 304 v_169(1:imx,jmx,1:kmx,1:ktmx)=v_in(1:imx,jmx_V,1:kmx,1:ktmx) 305 306 c Interpolation for the rest 307 308 do l=1,imx 309 do m=2,jmx-1 310 do n=1,kmx 311 do o=1,ktmx 312 v_169(l,m,n,o)=v_in(l,m-1,n,o)+ 313 & (v_in(l,m,n,o)-v_in(l,m-1,n,o))/2 314 enddo 315 enddo 316 enddo 317 enddo 318 319 c Extract T 320 321 start(1) = 1 322 start(2) = 1 323 start(3) = 1 324 start(4) = 1 325 326 count(1) = imx 327 count(2) = jmx 328 count(3) = kmx 329 count(4) = 1 330 331 do t = 1,ktmx 332 start(4)=t 333 Status = nf_get_vara_real (id, idT, start, count, ts_tmp) 334 call handle_err(Status) 335 ts_169(1:imx,1:jmx,1:kmx,t) = ts_tmp(1:imx,1:jmx,1:kmx) 336 enddo 337 338 339 c There is no need to interpolate temperature since this variable 340 c is avaliable at them RHO points. 341 342 c---------------------------------------------------------------------- 343 c Obtain starting and ending date of the data 344 c---------------------------------------------------------------------- 345 c MEDSLIK_II input files are hour files. These files' names have the 346 c following structure: nameYYMMDDHH, where name is a 4 digits 347 c description of the region where the data comes from. 348 c The main target of this work is to adapt the ROMS files so MEDSLIK_II 349 c can process them. Therefore, it is necessary to create these hour 350 c files and give them a name that follows the structure showed above. 351 c ROMS data, unlike MFS data, start at 00:00. The netcdf file contains 352 c a string (scrum_time) which stores the time since initialization. 353 c This variable will be used to obtain: inital date, final date and 354 c number of days. 355 c---------------------------------------------------------------------- 356 c Extract scrum_time (st) 357 358 start_st(1) = 1 359 start_st(2) = 1 360 Archivo: /home/paula/MEDSLIK_II_1.01/EXE/source/Extract_II.for Página 6 de 12 361 count_st(1) = ktmx 362 count_st(2) = 1 363 364 Status = nf_get_vara_int (id, idst, start_st, count_st, st) 365 call handle_err(Status) 366 367 c Starting month and day of the data 368 369 im=st(1,1)/86400 370 371 im1=im-31 372 im2=im-31-28 373 im3=im-31-28-31 374 im4=im-31-28-31-30 375 im5=im-31-28-31-30-31 376 im6=im-31-28-31-30-31-30 377 im7=im-31-28-31-30-31-30-31 378 im8=im-31-28-31-30-31-30-31-31 379 im9=im-31-28-31-30-31-30-31-31-30 380 im10=im-31-28-31-30-31-30-31-31-30-31 381 im11=im-31-28-31-30-31-30-31-31-30-31-30 382 im12=im-31-28-31-30-31-30-31-31-30-31-30-31 383 384 if (im1.le.0) then 385 month=1 386 day=im1+31 387 else if (im2.le.0) then 388 month=2 389 day=im2+28 390 else if (im3.le.0) then 391 month=3 392 day=im3+31 393 else if (im4.le.0) then 394 month=4 395 day=im4+30 396 else if (im5.le.0) then 397 month=5 398 day=im5+31 399 else if (im6.le.0) then 400 month=6 401 day=im6+30 402 else if (im7.le.0) then 403 month=7 404 day=im7+31 405 else if (im8.le.0) then 406 month=8 407 day=im8+31 408 else if (im9.le.0) then 409 month=9 410 day=im9+30 411 else if (im10.le.0) then 412 month=10 413 day=im10+31 414 else if (im11.le.0) then 415 month=11 416 day=im11+30 417 else if (im12.le.0) then 418 month=12 419 day=im12+31 420 else 421 print *, 'There is an error with the scrum time string' 422 endif 423 424 c Variable with the starting day of the data (YYMMDD) 425 426 write (monthc,'(i2)') month 427 write (dayc,'(i2)') day 428 429 if (month .lt. 10) then 430 if (day .lt. 10) then 431 Startd='150'//monthc(2:2)//'0'//dayc(2:2) 432 else