GISPART: a numerical model to simulate the dispersion of contaminants in the Strait of Gibraltar

Short Communication GISPART: a numerical model to simulate the dispersion of contaminants in the Strait of Gibraltar R. Peria ´n ˜ez* Departamento Fı´sica Aplicada I, E.U. Ingenierı´aTe ´cnica Agrı´cola, Universidad de Sevilla, Ctra. Utrera km 1, 41013-Sevilla, Spain Received 16 January 2004; received in revised form 28 May 2004; accepted 7 June 2004 Abstract GISPART is a three-dimensional particle-tracking code to simulate the dispersion of contaminants in the Strait of Gibraltar. It consists of a hydrodynamic module that is run off-line to determine tidal constants and residuals in the domain. This information is stored in several files that are read by the dispersion module to reconstruct water movements. Several application cases are presented, including instantaneous and continuous releases as well as different wind conditions, to show the information that can be extracted from the model. Ó2004 Elsevier Ltd. All rights reserved. Keywords: Gibraltar Strait; Particle tracking; Numerical modelling; Pollutants Software availability Program name: GISPART (GIbraltar Strait PARticle- Tracking model) Developer: R. Peria ´n ˜ez Contact: Dpto Fı´sica Aplicada 1, EUITA, Ctra Utrera km 1, 41013-Sevilla, Spain. E-mail: rperianez@ us.es Hardware: Compatible PC 200 MHz or higher Programe code: Fortran Cost: Free Availability: http://www.personal.us.es/rperianez or contact R. Peria ´n ˜ez 1. Introduction Recently, there has been an increasing interest in the development of pollutant dispersion models for the marine environment to use for decision making purposes after contaminant spills. In particular, particletracking methods are well suited for problems in which high contamination gradients are involved, since they minimize the effects of numerical diffusion. Also, they can be used to rapidly assess contaminant dispersion if the hydrodynamics are simulated previously off-line. Particle-tracking models have been used to simulate the dispersion of passive tracers (Gomez-Gesteira et al., 1999; Harms et al., 2000), radionuclides (Schonfeld, 1995; Peria ´n ˜ez and Elliott, 2002), oil spills (Proctor et al., 1994a,b) and even contaminated milk (Elliott et al., 2001) in several coastal water environments. The objective of this paper consists of giving a brief description of a new Lagrangian pollutant dispersion model developed for the Strait of Gibraltar, connecting the Mediterranean Sea and the Atlantic Ocean. This is a relevant topic because of the intense shipping activities in the Strait, which include the transport of radioactive material from/to the nuclear fuel reprocessing plants of Sellafield and Cap de la Hague (in UK and France, respectively), the transit of nuclear submarines, as well * Tel.: C34 954486474; fax: C34 954486436. E-mail address: [email protected] 1364-8152/$ - see front matter Ó2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.envsoft.2004.06.008 www.elsevier.com/locate/envsoft Environmental Modelling & Software 20 (2005) 797–802 as massive transport of chemical contaminants. Shipping routes are complex, with intersections of longitudinal routes with some 2000 annual transverse round trips Algeciras–Ceuta, Algeciras–Tanger and Tarifa– Tanger (see Fig. 1). Also, the port of Algeciras is the most important in Spain, with 60.7 Mt of cargo handled in 2003. Fishing activities in the area must also be considered. The area of the Strait of Gibraltar has a high ecological and tourist value, and there are also some important towns. Thus, a contamination release into the Strait as a consequence of an accident (or a deliberate release) can lead to a large ecological and economic impact. 2. Model description The GISPART model consists of two sub-models. First, a hydrodynamic module is run off-line. This provides the tidal constants and residuals that are required to reconstruct water movements in the model domain, and are stored in files that are read by the dispersion module to compute advective transport. Once the hydrodynamic module has been adequately calibrated and all information required by the dispersion computations is stored, it is not necessary to repeat the hydrodynamic calculations. This is an advantage over coupled hydrodynamic and dispersion models. Model equations are summarized in Appendix. 2.1. Hydrodynamic module An important feature of the tidal flow in the Strait is that it can be considered barotropic, as a first approach (Man ˜anes et al., 1998; Tsimplis and Bryden, 2000). As a consequence, 2D depth-averaged models have already been applied to simulate surface tides in the Strait (see for instance Sa ´nchez and Pascual, 1988; Tejedor et al., 1999). Tsimplis et al. (1995) have even used a 2D barotropic model for simulating tides in the whole Mediterranean Sea. The success of these models indicates that the baroclinic component is of secondary importance. The barotropic hydrodynamic equations are solved over the model domain using finite differences. Surface elevations are prescribed from observations along open boundaries and radiation conditions are used to determine the current component that is normal to the open boundary. A quadratic law for bottom friction is applied. Hydrodynamic calculations are carried out separately for the two main tidal constituents, M 2 and S 2 . Thus, spring–neap tidal cycles can be simulated. Once a stable periodic solution is achieved, standard tidal analysis is carried out and residual transport is calculated for each constituent. Tidal constants (amplitudes and phases) for each point in the domain and residual transports for each tidal constituent are stored in files to be read by the dispersion code. Results from the hydrodynamic calculations have been validated through an extensive comparison of tidal amplitudes, phases and current magnitudes and phases with observations for 16 points in the domain. Although not presented here, results are, in general, in good agreement with observations (Peria ´n ˜ez, 2004). 2.2. Dispersion code The dispersion of contaminants is calculated using a particle-tracking method. Essentially, the pollutant discharge is simulated by a number of discrete particles. The path followed by each particle is computed, turbulent diffusion being modelled as a three-dimen- sional random walk process. Decay of particles is also simulated using a Monte Carlo method (Peria ´n ˜ez and Elliott, 2002). This process is relevant, for instance, in the case of radioactive contaminants. The density of particles per water volume unit is finally computed to obtain contaminant concentrations over the Strait at the desired time. Both instantaneous and continuous releases of particles can be simulated. It must be noted that the particle-tracking model is three-dimensional, while the hydrodynamic module provides depth-aver- aged currents. Thus, a current profile is generated from the depth-averaged currents at each location by the dispersion code (see for instance Pugh, 1987). The spatial resolution of the dispersion model is DxZDyZ2500 m and time step is DtZ600 s. The effect of wind is included as usual in particletracking models. Thus, it is assumed that the water Fig. 1. Map of the computational domain showing some important towns (squares). Each unit in the xand yaxes is the grid cell number (thus equal to 2500 m). 798 R. Peria ´n ˜ez / Environmental Modelling & Software 20 (2005) 797–802 surface moves in the direction of wind at a speed equal to 3% of the wind speed 10 m above the sea surface. This current decreases logarithmically to zero at a depth usually taken as 20 m (Elliott, 1986; Pugh, 1987; Proctor et al., 1994a). Date and time of the discharge (and duration in the case of continuous releases) must be specified since the fate of the release will depend on the tidal state when it took place. Thus, the appropriate phase of each tidal constituent at tZ0 must be specified. The values used in this model correspond to the origin of time being January 1, 2003 at 0:15 h Greenwich time. The adsorption of pollutants by suspended and bottom sediments can also be simulated with a particle-tracking model (Peria ´n ˜ez and Elliott, 2002). However, these processes are neglected in the present study since suspended matter concentrations are very low in the Strait, typically 0.1–0.5 mg/L (Leo ´n-Vintro ´et al., 1999). Also, the average depth is 350 m (reaching 900 m in the eastern part) and, as a consequence, interactions of pollutants with bottom sediments can be neglected too. 2.3. Input data The dispersion code automatically reads 6 files that contain the topography of the Strait, amplitudes and phases of the two tidal constituents included in the model (M 2 and S 2 ), and the residual circulation. These files are provided by the hydrodynamic module that has been run off-line, as noted before. The user must also include some information for each simulation, related to the release characteristics. This information is summarized in Table 1. 2.4. Model output The model provides several output files. It gives the coordinates of each particle at several times during the simulation. Twelve snapshots at constant intervals during the simulation are plotted to show the evolution of the pollutant tracers over time. Another file contains a map of the final contaminant concentration over the Strait computed from the density of particles per water volume unit. Optionally, the time evolution of contaminant concentration at desired points over the domain can be obtained. 3. Application cases An example of the type of results that can be obtained from the particle-tracking model is presented in Fig. 2. Conservative particles (i.e. no decay) are considered in the simulations here. Instantaneous discharge of a pollutant was inputted into a grid cell in the area of Camarinal Sill during high water at Tarifa and with no wind. Three thousand particles are used in the simulation, whose tracks are followed during 2 days. The position of each particle at four different times after the release is shown in Fig. 2 (top). The concentration of the pollutant in arbitrary units per m 3 at tZ48 h is also presented in Fig. 2 (center). There is a net transport towards the Mediterranean Sea due to the residual currents, although the patch moves forward and backward following tidal oscillations. This can also be seen in Fig. 2 (bottom), where the time evolution of the number of particles inside an arbitrary grid cell [in this case (15,9)] is shown. The patch moves three times over this point, producing three peaks in the number of particles at 21, 26 and 36 h after the release. The highest peak, 254 particles, is observed 26 h after the release. In this simulation 1.0 !10 6 units of contaminant were released, thus the peak implies a maximum concentration equal to 9.2 !10 5 units/m 3 . For the following peak, at tZ36 h, the concentration is reduced by a factor 5 due to the spreading of the patch. From the center of the patch after 48 h an average velocity of the pollutant of 17 cm/s towards the Mediterranean can be estimated. This number can be compared with the mean speed of the Atlantic inflow equal to 23 cm/s measured by Tsimplis and Bryden (2000). The movement of a patch is obviously influenced by wind conditions. This can be clearly seen in Fig. 3. The same simulation described above has been repeated but with 15 m/s east and west winds, which is common in the Strait. The position of particles 28 h after the release for both simulations is shown in Fig. 3, which can be compared with the 28 h patch in Fig. 2. West winds, directed in the same direction as the residual circulation, produce a faster movement to the eastern part of the Strait, while east winds tend to retain particles into the Strait. Since the particle-tracking model is three-di- mensional, shear diffusion is produced and the patch size increases in the direction of wind. An example of the simulation of a continuous release is presented in Fig. 4. The release occurs at same point and tidal conditions as before (cell (7,9) and high water at Tarifa), and under calm wind. The position of particles 44 h after the release is shown in Fig. 4. This Table 1 Information required by the model to be introduced by the user Release point coordinates Select instantaneous/continuous release option Wind speed Wind direction Release date (day, month, year) Release time UTC, (hours, minutes) Simulation time (days) Magnitude of the release in the corresponding units Contaminant decay constant 799R. Peria ´n ˜ez / Environmental Modelling & Software 20 (2005) 797–802 can be compared with the 44 h patch in Fig. 2. Now there is a plume extending from the release point to the eastern part of the Strait. It is interesting to observe that four patches with larger concentrations of particles are apparent in the plume. They correspond to particles released during slack water, which remain concentrated and move together. Finally, an accident has been simulated at Algeciras harbor. This is a relevant topic, as mentioned in Section 1. The accident consisted of an instantaneous release occurring during high water at Tarifa. This would be the most critical accident, since contaminants are retained into Algeciras Bay. Only in the case of winds blowing from the north would contaminants be flushed out of the bay. The position of particles for calm conditions and a 15 m/s northern wind, 14 days and 2 days after the Fig. 3. Position of particles 28 h after an instantaneous release at grid cell (7,9) with 15 m/s winds from the east and west. The release point and wind directions are indicated. Fig. 4. Position of particles 44 h after the beginning of a continuous release at grid coordinates (7,9). Fig. 2. Up: Dispersion of an instantaneous release. The position of particles at different times after the release is shown. Center: Computed surface concentration in arbitrary units per m 3 48 h after the release. It has been obtained from the density of particles. Down: Time evolution of the number of particles inside grid cell (15,9), whose position is indicated in the central panel. 800 R. Peria ´n ˜ez / Environmental Modelling & Software 20 (2005) 797–802 release respectively, is presented in Fig. 5. It can be seen that for calm wind essentially all the contamination is retained in the bay (this would also be the case for east and west winds). In the case of north winds, contaminants are transported to the axis of the Strait. From this area they are flushed out by longitudinal currents. As an indication, the running time for an instantaneous release (3000 particles are tracked) is approximately 10 s per day of simulation on a 333 MHz Pentium 2 computer. 4. Conclusions A particle-tracking model for simulating pollutant dispersion in the Strait of Gibraltar has been developed. The model solves the depth-averaged hydrodynamic equations for the M 2 and S 2 tidal constituents off-line. Tidal analysis is carried out and tidal constants and residuals are stored in files that are read by the dispersion model. Dispersion is solved using a Lagrangian approach, diffusion and decay being simulated by means of a Monte Carlo method. The dispersion model is three-dimensional, thus standard vertical profiles for the tidal and wind-induced currents have been used to define the variation of currents with depth. Some examples on the dispersion of contaminants have been provided. Generally, the fate of a patch depends on the tidal state when the release was carried out and on wind conditions. East winds oppose the residual current and tend to retain contaminants in the Strait, thus greater contaminant pollution occurs. Under calm conditions and, especially, west winds, the Strait is flushed off rapidly. The most critical accident would take place at, or near, Algeciras harbor. In this case contamination stays in Algeciras Bay, only if wind blows from the north sector a faster flushing of contaminants out of the bay occurs. Appendix The depth-averaged hydrodynamic equations are written as: vz vtCv vxðHuÞCv vyðHvÞZ0ð1Þ vu vtCuvu vxCvvu vyCgvz vxUvCtu rH ZA v2u vx2Cv2u vy2!ð2Þ vv vtCuvv vxCvvv vyCgvz vyCUuCtv rH ZA v2v vx2Cv2v vy2!ð3Þ where uand vare the depth-averaged water velocities along the xand yaxes, Dis the depth of water below the mean sea level, zis the displacement of the water surface above the mean sea level measured upwards, HZDCz is the total water depth, Uis the Coriolis parameter (UZ2wsin b, where wis the earth rotational angular velocity and bis latitude), gis acceleration due to gravity, ris water density and Ais the horizontal eddy viscosity. t u and t v are friction stresses. The surface current is deduced from the depth averaged one: usZmC1 muð4Þ where uis the depth-averaged current provided by the hydrodynamic model and mZ6. The current speed at a level z#below the sea surface is: Fig. 5. Position of particles after an accident at Algeciras harbor, 14 days after the accident in the case of calm conditions and 2 days after the accident in the case of a 15 m/s north wind. The release point and wind directions are indicated. 801R. Peria ´n ˜ez / Environmental Modelling & Software 20 (2005) 797–802 uz#ZusDz# D1=m ð5Þ The wind-induced current at any depth is: uz#Z(u0u) klnz# z0if z#!20 m 0ifz#R20 m ð6Þ where u 0 is the surface wind-induced current, kZ0.4 is the von Karman constant, u* is a friction velocity and z 0 is the sea surface roughness length. The friction velocity can be estimated as u)Z0:0012Wð7Þ for a wide range of conditions, where Wis wind speed 10 m above the sea surface. Of course, these current profiles are solved in the uand vcomponents. Advection is computed solving the following equation for each particle: dr dt Zqð8Þ where ris the position vector of the particle and qis the current vector (due to wind and tide) at the particle position solved in components uand v. The maximum sizes of the horizontal and vertical steps due to turbulent diffusion, D h and D v respectively, are: DhZffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 12KhDt p DvZffiffiffiffiffiffiffiffiffiffiffiffiffi 2KvDt pð9Þ where K h and K v are the horizontal and vertical diffusion coefficients respectively. References Elliott, A.J., 1986. Shear diffusion and the spread of oil in the surface layers of the North Sea. German Journal of Hydrography 39, 113–137. Elliott, A.J., Wilkins, B.T., Mansfield, P., 2001. On the disposal of contaminated milk in coastal waters. Marine Pollution Bulletin 42, 927–934. Gomez-Gesteira, M., Montero, P., Prego, R., Taboada, J.J., Leitao, P., Ruiz-Villareal, M., Neves, R., Perez-Villar, V., 1999. A twodimensional particle-tracking model for pollution dispersion in A Corun ˜a and Vigo Rias (NW Spain). Oceanologica Acta 22, 167–177. Harms, I.H., Karcher, M.J., Dethleff, D., 2000. Modelling Siberian river runoff – implications for contaminant transport in the Artic Ocean. Journal of Marine Systems 27, 95–115. Leo ´n-Vintro ´, L., Mitchell, P.I., Condren, O.M., Downes, A.B., Papucci, C., Delfanti, R., 1999. Vertical and horizontal fluxes of plutonium and americium in the western Mediterranean and the Strait of Gibraltar. The Science of the Total Environment 237, 77–91. Man ˜anes, R., Bruno, M., Alonso, J., Fraguela, B., Tejedor, L., 1998. Non-linear interaction between tidal and subinertial barotropic flows in the Strait of Gibraltar. Oceanologica Acta 21, 33–46. Peria ´n ˜ez, R., 2004. A particle-tracking model for simulating pollutant dispersion in the Strait of Gibraltar. Marine Pollution Bulletin 49, 613–623. Peria ´n ˜ez, R., Elliott, A.J., 2002. A particle tracking method for simulating the dispersion of non-conservative radionuclides in coastal waters. Journal of Environmental Radioactivity 58, 13–33. Proctor, R., Flather, R.A., Elliott, A.J., 1994a. Modelling tides and surface drift in the Arabian Gulf: application to the Gulf oil spill. Continental Shelf Research 14, 531–545. Proctor, R., Elliott, A.J., Flather, R.A., 1994b. Forecast and hindcast simulations of the Braer oil spill. Marine Pollution Bulletin 28, 219–229. Pugh, D.T., 1987. Tides, Surges and Mean Sea Level. Wiley, Chichester. Sa ´nchez, P., Pascual, J.R., 1988. Primeras experiencias en la modelizacio ´n del Estrecho de Gibraltar. In: Almaza ´n, J.L. (Ed.), Oceanografı ´aFı ´sica del Estrecho de Gibraltar. Escuela Superior de Ingenieros de Caminos, Madrid, pp. 251–282 (in Spanish). Schonfeld, W., 1995. Numerical simulation of the dispersion of artificial radionuclides in the English Channel and the North Sea. Journal of Marine Systems 6, 529–544. Tejedor, L., Izquierdo, A., Kagan, B.A., Sein, D.V., 1999. Simulation of the semidiurnal tides in the Strait of Gibraltar. Journal of Geophysical Research 104, 13541–13557. Tsimplis, M.N., Bryden, H.L., 2000. Estimations of the transports through the Strait of Gibraltar. Deep Sea Research 47, 2219–2242. Tsimplis, M.N., Proctor, R., Flather, R.A., 1995. A two dimensional tidal model for the Mediterranean Sea. Journal of Geosphysical Research 100, 16223–16239. 802 R. Peria ´n ˜ez / Environmental Modelling & Software 20 (2005) 797–802