Modelling the environmental behaviour of pollutants in Algeciras Bay (south Spain)

Modelling the environmental behaviour of pollutants in Algeciras Bay (south Spain) R. Periáñez ⇑ Dpto. Física Aplicada I, ETSIA, Universidad de Sevilla, Ctra. Utrera km 1, 41013 Sevilla, Spain article info Keywords: Numerical modelling Hydrodynamics Sediments Dispersion Heavy metals Algeciras Bay abstract An environmental study of Algeciras Bay is carried out through numerical modelling. First, a 2D barotropic model is applied to calculate tides and mean circulation. Results of this model are used by a sediment transport model which provides suspended matter concentrations and sedimentation rates in the Bay. It includes three particle classes. An effective diffusion coefficient has been calibrated simulating temperature distribution inside the Bay. An additional validation is obtained from an independent nitrate dispersion simulation. Then heavy metal dispersion patterns are investigated using a model which includes water–sediment metal interactions and uses the outputs of the hydrodynamic and sediment transport models. The metal transport model has been applied to simulate the dispersion of Zn, Cu and Ni. Results from the hydrodynamic, sediment and metal transport models have been compared with measurements. Model results also indicate that transport inside the Bay is relatively weak. Numerical experiments have been carried out to determine flushing times for conservative and non-conservative pollutants. Flushing time is about 20 days for a conservative tracer, and this value is mainly due to the M 2 residual current. Tides are not effective in removing pollutants. Ó2011 Elsevier Ltd. All rights reserved. 1. Introduction Algeciras Bay (south Spain) is located in Gibraltar Strait, which connects the Atlantic Ocean and the Alboran Sea, in the Mediterranean (Fig. 1 top). Large amounts of pollutants are discharged into the Bay: urban waste from some towns and industrial wastes from one of the largest industrial areas in Andalusia, with a significant number of petrochemical plants and refineries, in addition to the production of steel, paper and power (thermal power plants). The Palmones and Guadarranque rivers empty into the Bay (Fig. 1 bottom), and the Bay is also polluted by heavy maritime traffic (Algeciras port is one of the largest in Spain). In previous work (Morillo et al., 2008), it has been stated that water in the Bay has a high turnover because of its proximity to the Strait of Gibraltar, where the Mediterranean Sea and the Atlantic Ocean meet, and the strong currents that predominate in the area. However, some preliminary calculations with a low resolution model of the Strait of Gibraltar (Periáñez, 2004b) indicated that contaminants released into the Bay tend to remain there. Only north winds are efficient for a rapid cleaning of the Bay. The objective of this work consists of developing a high resolution numerical model of Algeciras Bay which may be used to understand the dynamics of pollutants in this environment, and to carry out a quantitative study of flushing times. In particular, we are interested in assessing how water circulation affects such flushing times. The model consists of a hydrodynamic module, which provides water circulation, a sediment transport model which provides suspended matter concentrations and sedimentation rates, and the pollutant dispersion model. The dispersion model may be applied to non-conservative pollutants (do not remain in solution). Uptake/release reactions between the dissolved and solid phases (suspended matter and bed sediments) are described in a dynamic way, using kinetic transfer coefficients. Results from the hydrodynamic and sediment transport models have been compared with measurements in the Bay. The pollutant dispersion model has been applied to simulate the behaviour of a conservative substance (nitrates) and to heavy metals. In the case of metals, computed concentrations in water and bed sediments are compared with measurements. Model results indicate that, in general, the dynamics of pollutants in the Bay has been correctly described. Then, numerical experiments are carried out to estimate flushing times for both conservative and non-conservative substances. The effects of the different tidal constituents and residuals on flushing times have been investigated. It is worth commenting that this is the first study on environmental conditions, pollutant dynamics and flushing carried out for Algeciras Bay with numerical modelling. Consequently this study completes previous field work carried out in the area (Morillo et al., 2008; Sánchez-Moyano et al., 2002; Guerra-García et al., 2009; Morales-Caselles et al., 2007; Morillo and Usero, 2008; Estacio et al., 1997). 0025-326X/$ - see front matter Ó2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.marpolbul.2011.11.030 ⇑ Tel.: +34 954486474; fax: +34 954486436. E-mail address: [email protected] Marine Pollution Bulletin 64 (2012) 221–232 Contents lists available at SciVerse ScienceDirect Marine Pollution Bulletin journal homepage: www.elsevier.com/locate/marpolbul The three submodels (hydrodynamics, sediments, pollutants) are described separately in the next section, together with some indications on the numerical solution of equations. Later results are presented and discussed, again separately for each submodel. A study on flushing times closes the paper. 2. Model description 2.1. Hydrodynamic model An important feature of the tidal flow in the Strait is that it can be considered, as a first approach, as barotropic. Indeed, 93% of the current velocity variance in the semidiurnal band has a barotropic character in the Strait (Mañanes et al., 1998). As a consequence, 2D depth-averaged models have already been applied to simulate tides in the Strait (Tejedor et al., 1999; Periáñez and Pascual- Granged, 2008). Tsimplis et al. (1995) have even used a 2D barotropic model for simulating tides in the whole Mediterranean Sea. Other authors have also stated that using a 2D model for simulating tides is a reasonable approach (Dyke, 2001; Yanagi, 1999). Inside Algeciras Bay, it has been found that salinity and temperature are essentially homogeneous both horizontally and vertically (University of Cadiz, 2007). Thus, the use of a 2D depth averaged barotropic model is justified. It is based on the following equations: @f @tþ@ @xðHuÞþ @ @yðH v Þ¼0ð1Þ @u @tþu@u @xþ v @u @yþg@f @xX v þ s u q w H¼A@ 2 u @x 2 þ@ 2 u @y 2 !ð2Þ @ v @tþu@ v @xþ v @ v @yþg@f @yþXuþ s v q w H¼A@ 2 v @x 2 þ@ 2 v @y 2 !ð3Þ where uand vare the depth averaged water velocities along the x and yaxis, his the depth of water below the mean sea level, fis the displacement of the water surface above the mean sea level measured upwards, H¼hþfis the total water depth, Xis the Coriolis parameter ðX¼2wsin k, where wis the Earth rotational angular velocity and kis latitude), gis acceleration due to gravity, q w is water density and Ais the horizontal eddy viscosity. s u and s v are friction stresses which have been written in terms of a quadratic law: s u ¼k q w uffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 þ v 2 pð4Þ s v ¼k q w v ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 þ v 2 p where kis the bed friction coefficient. The solution of these equations provides the water currents at each point in the model domain and for each time step. Currents are treated through standard tidal analysis (Pugh, 1987, Chapter 4) and tidal constants are stored in files that will be read by the dispersion codes to calculate the advective transport. The model includes the two main tidal constituents, M 2 and S 2 . Thus, the hydrodynamic equations are solved for each constituent and tidal analysis is also carried out for each constituent separately. A residual transport cannot be produced by the pure harmonic currents given by the tidal analysis, thus tidal residuals have been calculated as well, as explained below. Diurnal constituents are not included since most of the variance of current velocities is given by the M 2 signal alone, S 2 being the second important constituent. Therefore they can be used to characterize a very significant fraction of tides in the area (Mañanes et al., 1998). Once a stable periodic solution of hydrodynamic equations is achieved, tidal analysis is carried out to determine tidal constants that are used by the sediment and dispersion codes. Tidal residual transports are also calculated. This is done for the M 2 and S 2 tides separately. Tidal residuals for each constituent are calculated from the equation: ~ q r ¼hH~ q t i hHið5Þ which corresponds to the Eulerian residual transport velocity (Delhez, 1996). In this equation hiis the time averaging operator, ~ q r is the tidal residual and ~ q t is the instantaneous tidal current. In addition, the mean current due to geostrophic flow is calculated from an additional run of the hydrodynamic model. In this run constant values of the sea surface elevation along the open boundaries of the domain are specified (Sannino et al., 2004; Periáñez and Pascual-Granged, 2008). The total mean current is obtained adding the tidal residuals of the M 2 and S 2 constituents plus the mean geostrophic flow. Current ~ qat any time and position in the Strait is obtained adding this total mean current and the instantaneous tidal currents at the corresponding point. This is the water velocity value used in the dispersion calculations. This procedure is the same used in Periáñez and Pascual-Granged (2008). 2.2. Sediment transport The transport of sediments is described by a 2D advection– diffusion equation to which some terms are added. These are external sources of particles, terms describing particle deposition on the seabed and erosion from the bed to the water column. Fig. 1. Up: general localization of Algeciras Bay (inside the red circumference). Down: topography (m) of the Bay. Each unit in the axis corresponds to 100 m (grid cell spacing). |indicate urban waste release points and }industrial waste release points. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 222 R. Periáñez / Marine Pollution Bulletin 64 (2012) 221–232 The formulation of these processes is based upon standard formulae. Thus, the erodability constant is used for the erosion term. Particle deposition is described using the settling velocity, which is obtained from Stoke’s law. Critical erosion and deposition stresses are applied as usually (Periáñez, 2005a,b; Liu et al., 2002a; Lumborg and Windelin, 2003; Cancino and Neves, 1999). Three particle classes are considered in the model, each one characterized by a mean size, to obtain a better representation of suspended particles dynamics and sedimentation in the Bay. The equation for suspended sediment transport is: @ðHm k Þ @tþ@ðuHm k Þ @xþ@ð v Hm k Þ @y¼@ @xHK h @m k @x  þ@ @yHK h @m k @y  þðER k DEP k ÞþS k ð6Þ where k¼1;2;3 denote each sediment class, m k is the suspended matter concentration, K h is an effective horizontal diffusion coefficient (as explained in Section 3.2), S k is the external particle source and ER k and DEP k are the erosion and deposition terms, respectively. The deposition term is written in the following form: DEP k ¼w k m k 1j~ s j s cd  ð7Þ where ~ s is the total bottom friction stress (tides plus mean flow) and s cd is a critical deposition stress above which no deposition occurs since particles are kept in suspension by turbulence. The settling velocity of particles is obtained from Stokes’s law as commented above: w k ¼ q  q w q w gD 2 k 18 m ð8Þ where qand D k are suspended particle density and diameter, respectively, and mis the kinematic viscosity of water. The erosion rate is written in term of the erodability constant: ER k ¼Ef k j~ s j s ce 1  ð9Þ where Eis the erodability constant, f k gives the fraction of kclass particles in the bed sediment and s ce is a critical erosion stress below which no erosion occurs. The model can also calculate sedimentation rates (SR) as the balance between the deposition and erosion terms. 2.3. Pollutant transport Non-conservative pollutants are those which do not remain dissolved in the water column, but have a certain affinity to be fixed to particles. If the pollutant is introduced in the water column, it will be fixed to settling suspended particles and their deposition on the sea bed will contaminate the bottom sediment. Of course there are also advection/diffusion processes in water and direct adsorption of pollutants on the seabed. The exchanges between the dissolved and solid phases may be described in terms of kinetic transfer coefficients. Thus, assuming that adsorption/release reactions are governed by a single reversible reaction, a coefficient k 1 governs the transfer from the liquid to the solid phase and a coefficient k 2 governs the inverse process. Also, the migration of contaminants to the deep sediment may be included. Thus, contaminants deposited on the sediment surface will be buried by particle deposition and will migrate below the mixed sediment layer which directly interacts with the dissolved phase (see for instance the discussion in Monte et al., 2006). This effect may be easily treated as a decay process with constant k burial given by Periáñez (2009): k burial ¼SR q s Lð10Þ where Lis the sediment mixing depth (the distance to which the dissolved phase penetrates the sediment) and q s is the sediment bulk density (dry mass divided by wet volume). Nevertheless, this process is not significant in Algeciras Bay due to the low SR values and the typical simulation times (several months). Sediment bulk density is related to particle density through porosity p: q s ¼qð1pÞ. The adsorption process is a surface phenomenon that depends on the surface of particles per water volume unit into the grid cell. This quantity has been denoted as the exchange surface (Periáñez, 2003a, 2004a, 2008, 2009). Thus in general: k 1;k ¼ v ðS m;k þS s;k Þ¼k spm 1;k þk sed 1;k ð11Þ where S m;k and S s;k are the exchange surfaces for suspended matter and bottom sediments, respectively (dimensions ½L 1 ) for the corresponding particle class k. v is a parameter with the dimensions of a velocity. It is denoted as the exchange velocity (Periáñez, 2004a, 2008, 2009). Assuming spherical particles, the exchange surfaces are written as (see references cited above): S m;k ¼3m k q R k ð12Þ and S s;k ¼3Lf k ð1pÞ/ R k Hð13Þ where R k is particle radius, pis sediment porosity and /is a correction factor that takes into account that part of the sediment particle surface may be hidden by other sediment particles. This formulation has been successfully used in all modelling works cited above. Real particles are not spheres, but with this approach it is possible to obtain an analytical expression for the exchange surface (Duursma and Carroll, 1996). The equation that gives the temporal evolution of pollutant concentration in the dissolved phase, C d , is: @ðHC d Þ @tþ@ðuHC d Þ @xþ@ð v HC d Þ @y¼@ @xHK h @C d @x  þ@ @yHK h @C d @y  H X k k spm 1;k þk sed 1;k  C d þk 2 H X k m k C s;k þk 2 L q s /X k f k A s;k ð14Þ where A s;k and C s;k are concentrations in bed sediment and suspended matter class k, respectively. The temporal evolution of pollutant concentration in suspended particles is given, for each class k, by: @ðHC s mÞ k @tþ@ðuHC s mÞ k @xþ@ð v HC s mÞ k @y ¼@ @xHK h @ðC s mÞ k @x  þ@ @yHK h @ðC s mÞ k @y  þk spm 1;k C d H k 2 m k C s;k HþSED k ð15Þ where SED k expresses the pollutant exchange between suspended particles and the bed sediment resulting from erosion/deposition: R. Periáñez / Marine Pollution Bulletin 64 (2012) 221–232 223 SED k ¼SR k C s;k SR k >0 SR k A s;k SR k <0 ð16Þ The equation for the temporal evolution of concentration in the bed sediment is: @A s;k @t¼k sed 1;k C d H L q s f k k 2 A s;k /þSED k ð17Þ where now the exchange due to erosion/deposition of suspended particles is written as SED k ¼ SR k C s;k L q s f k SR k >0 SR k A s;k L q s f k SR k <0 8 < :ð18Þ The total concentration of pollutants in the sediment, A tot ,is computed from: A tot ¼X k f k A s;k ð19Þ and in the case of suspended matter it is: C total s ¼P k m k C s;k P k m k ð20Þ 2.4. Numerical solution All equations are solved using explicit finite difference schemes (Kowalick and Murty, 1993) on a grid with resolution D x¼ D y¼100 m. Second order accuracy schemes are used for advective and diffusive terms. Time step in the hydrodynamic equations, limited by the CFL condition, is D t¼0:5 s. Since the sediment transport and pollutant dispersion models are run off-line, time step could be increased to 30 s. Open boundary conditions consist of, in hydrodynamic calculations, specifying sea surface elevations compiled from observations and from a larger model covering the complete Strait of Gibraltar (Tejedor et al., 1999; Periáñez and Caravaca, 2010). The suspended matter model (without the erosion term) is started from a sea bottom containing no sediments. Then the accumulation of particles of each class is calculated to have a first estimation of the distribution of sediment particle sizes over the model domain. Next, the suspended matter model is started (with erosion) from the estimated particle class distribution until a Fig. 2. Maps of computed tide amplitude and tidal current amplitude for both constituents included in the model. 0 50 100 150 200 250 0 50 100 150 200 250 X Y • Pta Carnero • Algeciras • Gibraltar • Pta Europa • Interior Fig. 3. Map indicating points where results from the hydrodynamic calculations have been compared with observations. Table 1 Observed, index obs, and computed, index comp, amplitudes (A, m) and phases (g, deg) of tidal elevations at several locations indicated in Fig. 3. Station M 2 S 2 A obs g obs A comp g comp A obs g obs A comp g comp Pta Carnero 0.311 47.5 0.318 46.3 0.115 71.0 0.119 71.2 Algeciras 0.310 48.0 0.291 45.5 0.111 73.9 0.112 72.4 Gibraltar 0.298 46.0 0.291 45.5 0.107 72.0 0.112 72.4 224 R. Periáñez / Marine Pollution Bulletin 64 (2012) 221–232 steady state is reached. This way a self-consistent distribution of different particle classes, f k on the sea bed over the model domain can be obtained. This procedure was successfully applied in Periáñez (2005a). Suspended matter concentrations are defined along open boundaries for solving the sediment transport model (Dirichlet boundary condition) according to measured values in the Strait of Gibraltar (León-Vintró et al., 1999). Once currents, suspended matter concentrations and sedimentation rates have been calculated and stored in files, the pollutant transport model may be run. In this case a no-gradient open boundary condition has been used. 3. Results Results from the hydrodynamic, sediment transport and pollutant dispersion models are presented separately in the following subsections. Next, results from the flushing numerical experiments are discussed. 3.1. Hydrodynamics As an example, computed tide amplitudes and tidal current amplitudes for M 2 and S 2 constituents are presented in Fig. 2.It may be seen that there is a significant reduction in the tide amplitude from the Strait into Algeciras Bay and towards the Mediterranean. Tidal currents are also weak inside the Bay, specially in the case of S 2 . Indeed, negligible currents are obtained for this constituent in the inner part of the Bay. A comparison between computed and measured (Tsimplis et al., 1995; Tejedor et al., 1999) tide amplitudes and phases at some points indicated in Fig. 3 may be seen in Table 1. Current amplitudes are also compared with available observations (Universidad de Cadiz, 2007) in Table 2. It may be seen that, generally speaking, the model is giving a realistic representation of tides in the Bay. A map showing the computed M 2 residual transport velocity may be seen in Fig. 4. Currents are of the order of 10 cm/s in the Strait of Gibraltar, but extremely weak inside the Bay. It is interesting to notice the cyclonic gyre which is apparent at the mouth of Algeciras Bay. 3.2. Sediment transport Winds and waves are not explicitly included in the hydrodynamic model. It is evident that they will produce an enhanced mixing inside the Bay (although waves inside the Bay are small, predominantly <0.5 m height; University of Cadiz, 2004). Thus an effective diffusion coefficient K h has been defined, which would include mixing induced by wind and waves. To select its value, the water temperature distribution has been simulated. Water temperature Tis described by an advection/diffusion equation in the form: @T @tþu@T @xþ v @T @y¼@ @xK h @T @x  þ@ @yK h @T @y  ð21Þ The model is started from an uniform water temperature over all the domain, but temperature is defined at the release points of the thermal power plants which are cooled with water from the Bay (these are points 4 and 6 in Fig. 5). The effects of these heat emissions are restricted to the shallow waters along the coastline. Thus, we assume that they do not affect the general circulation of the Bay and the barotropic approach is still valid. Indeed, it was commented in Section 2.1 that measurements of salinity and temperature showed generally homogeneous distributions (horizontally and vertically). These measurements were made in the main water body of the Bay and do not reflect the heat emissions from power plants. The model is run until a steady temperature distribution is obtained. The diffusion coefficient K h is changed by trial and error until the computed distribution agrees with measurements. When this occurs we would have the value of the effective diffusion coefficient which implicitly includes the average mixing induced by wind and waves. Such value resulted 10 m 2 =s, and a comparison between measured (Guerra-García et al., 2009) and computed water temperatures may be seen in Fig. 5. This approach consisting of not including winds explicitly in the simulation has been successfully used in previous modelling works (Periáñez, 2002, 2003a,b, 2005a, 2008, 2009). Although it may be a too simple approach, it can give a realistic view of the mean transport processes in the Bay, given the generally good agreement between measured and computed suspended sediment concentrations and metal distributions that will be shown below. However, the sensitivity of suspended sediment and pollutant transport in Algeciras Bay to changing winds is an interesting problem that has to be addressed in the next future, since dispersion processes related to meteorological conditions may become apparent. Moreover, erosion events may be induced by wind waves in the shallower areas. Nevertheless, wave-induced erosion is not included in the model since calculations are carried out under calm conditions. Once the effective diffusion coefficient has been defined, sediment transport simulations may be carried out. Three particle classes have been used. Mean particle radius defining each class are 0.1, 1.0 and 15 l m. The smallest class represents the background suspended sediment in the water column (as in Liu et al., 2002b). The following classes correspond to clay and medium silt, similarly to, for instance, Liu et al. (2002b),Jiang et al. (2000) and Holt and Table 2 Observed, index obs, and computed, index comp, magnitudes (cm/s) of uand v tidal currents at some locations indicated in Fig. 3. Station M 2 S 2 u obs v obs u comp v comp u obs v obs u comp v comp Pta Europa 24.3 19.3 28.3 19.6 7.26 8.06 8.50 4.20 Pta Carnero 3.33 13.5 1.40 12.9 1.73 1.65 1.80 6.70 Interior 2.88 0.55 0.10 0.30 0.53 0.40 0.10 0.10 0 50 100 150 200 250 0 50 100 150 200 X Y 10 cm/s Fig. 4. Computed M 2 residual. Only one of each 100 computed vector is drawn for clarity. R. Periáñez / Marine Pollution Bulletin 64 (2012) 221–232 225 James (1999). Settling velocities are then obtained from Stoke’s law, as mentioned above. Other parameters are required to run the sediment transport module. The critical deposition stress typically ranges between 0.04 and 0:1N=m 2 (Tattersall et al., 2003), while the critical erosion stress ranges 0:1—1:5N=m 2 (Tattersall et al., 2003). In the present application intermediate values of 0.06 and 1:0N=m 2 have been taken for the critical deposition and erosion stresses, respectively, as in the models for southern Spain described in Periáñez (2008, 2009). The erodability constant is fixed as E¼1:6 10 3 kg=m 2 s. This parameter typically varies between 2 10 4 and 3 10 3 kg=m 2 s(Tattersall et al., 2003). A standard value q ¼2600 kg=m 3 has been used for particle density. Sediment discharges by Palmones and Guadarranque rivers, as well as from waste release have been defined by trial and error, until the model showed a reasonable agreement with observations. Of course, these magnitudes present seasonal (and at a shorter scale as well) variations. However, these are not quantified and our objective has consisted, simply, of having a correct estimate of mean suspended matter concentrations over the domain to be able to simulate the dispersion of reactive pollutants. It must be considered that available sediment concentrations and SR to be compared with model output are mean (annual) averages (Sánchez-Moyano et al., 2002), which justifies our approach. The sediment model is started from an uniform suspended matter concentration of 3 g=m 3 (total) and integrated until a steady state is achieved (see also Section 2.4). The computed distribution of suspended sediments over the domain may be seen in Fig. 6 and a comparison between computed and measured (Sánchez-Moyano et al., 2002) suspended matter concentrations and sedimentation rates is presented in Fig. 7. The distribution in Fig. 6 shows an increase in suspended matter towards the interior of the Bay due to inputs from rivers and urban wastes, as stated in Sánchez-Moy- ano et al. (2002). The magnitude of suspended sediment concentration is in general correctly predicted by the model, except along the southwest shore of the Bay (points 1 and 2), where suspended matter concentration is underestimated. Sedimentation rates are, in general, underestimated. This is not surprising given the approximations made. Although constant suspended matter concentrations have been defined in the rivers, there will be seasonal variations depending on pluviometry, for instance. Heavy rain episodes will enhance mean sedimentation rates, and this is not described in the model. Sediment bed load transport is neglected. Finally, atmospheric deposition events of particles coming from the Sahara Desert have not been considered in the model since they cannot be easily quantified (as in Periáñez, 2008). 3.3. Pollutant transport The dispersion of metals Zn, Cu and Ni has been simulated since there are field measurements which can be used to test the model results. These three metals were studied in a previous modelling work of the Gulf of Cadiz (Periáñez, 2009), at the west side of Gibraltar Strait. Thus, the same parameters as in such model have 0 100 200 0 50 100 150 200 250 X Y • 1 • 2 • 3 • 4 • 5 • 6 • 7 • 8 1 2 3 4 5 6 7 8 14 16 18 20 22 Sampling point T (°C) Fig. 5. Computed and measured water temperature along the Bay’shoreline. Measurement points are indicated on the top panel. Fig. 6. Computed distribution of total suspended matter concentrations ðg=m 3 ) over the domain. A logarithmic scale is used to appreciate details more clearly. Numbers indicate points were suspended matter concentrations and sedimentation rates have been compared with measurements. 226 R. Periáñez / Marine Pollution Bulletin 64 (2012) 221–232 been used in the present work. Essentially, we require the exchange velocity v and the kinetic coefficient k 2 for each metal. Also, a mean porosity, sediment mixing depth Land correction factor / have to be specified. The exchange velocity is deduced from the equilibrium partition coefficient k d as described in detail in Periáñez (2009). Values for the different parameters are summarized in Table 3, and detailed justification may be seen in Periáñez (2009), thus it is not repeated here. Sensitivity analysis of model formulation to these parameters may also be seen in such reference. Nevertheless, a first numerical experiment has been carried out applying the model to simulate the dispersion of a conservative contaminant (thus v ¼k 2 ¼0). Releases from Acerinox plant have been simulated. It is a stainless steel production plant located between Palmones and Guadarranque rivers. Wastes released to the Bay contain significant amounts of nitrates (University of Cadiz, 2007). The model is started from zero concentrations over the whole domain and nitrate concentration in the release point from Acerinox has been fixed as 11.5 mg/l, according to University of Cadiz (2007). A simulation over 40 days is carried out. In this time steady state is reached. Model results are compared with the observed nitrate distribution in September 2006, when winds in the area were weak. It must be commented that concentration in the release point, 11.5 mg/l, is the measured one in this date. A \comparison between the observed and calculated distributions may be seen in Fig. 8. Measurements indicate that the 1 mg/l isoline is within a 0.5 miles radius circumference centered in the industrial plant. Model results are in agreement with these observations. The nitrate patch has an almost circular shape, due to the 1 2 3 4 5 6 7 8 9 10 0 5 10 15 20 25 Sampling point SPM (mg/l) 1 2 3 4 5 6 7 8 9 10 0 10 20 30 40 Sampling point SR (kg/m2month) Fig. 7. Comparison between measured (points) and computed (lines) suspended matter concentrations and sedimentation rates. Sampling points are indicatedinFig. 6. Measured values correspond to annual mean values (one sample per month). Table 3 Summary of model parameters. Details and justification may be seen in Periáñez (2009). Parameter description Value sediment mixing depth L¼0:1m particle density q¼2600 kg=m 3 correction factor /¼0:1 sediment porosity p¼0:5 desorption kinetic coefficient k 2 ¼1:16 10 5 s 1 Zn distribution coefficient k d ¼7:010 4 Ni distribution coefficient k d ¼2:010 4 Cu distribution coefficient k d ¼2:610 4 Fig. 8. Top: measured nitrate concentrations. The 1 mg/l concentration isoline is shown, as well as the 0.5 miles circumference centered in Acerinox. Bottom: computed concentrations. Lines are the 5 (closer to the source) and 1 mg/l isolines. R. Periáñez / Marine Pollution Bulletin 64 (2012) 221–232 227 very weak currents in this part of the Bay. Also, the dilution degree, in close agreement with the observed one, confirms the correct calibration of the effective diffusion coefficient previously carried out through the temperature simulations. To simulate reactive metal dispersion the model is again started from zero concentrations in water, suspended matter and bed sediments and integrated until steady concentrations are obtained over the domain in the three phases. In practice, simulations over some 100 days are required, longer than in the case of nitrates because of water–sediment interactions and the presence of multiple sources. Metal concentrations in water and bed sediments are compared with measurements in the Bay (Morillo et al., 2008; Morales-Caselles et al., 2007). A mean metal discharge, in dissolved form, is introduced from industrial and urban waste sources. For each metal, the same magnitude is used for all urban sources and for all industrial sources. These average inputs were selected by calibration. Of course, this is a very rough approach: inputs are not the same for all sources and will present temporal variability. However, it allows obtaining an estimation of average metal concentrations in the Bay and it may also be observed if the metal partition between water and sediment is correctly reproduced by the model. Using a different source magnitude for each release point, obtained from calibration, may look tricky since such magnitudes could be fitted to improve model agreement with observations. Thus it was decided to use the same magnitude for sources of the same nature (industrial or urban), although this will clearly make model/measurement comparisons worse. It must finally be taken into account that sediments integrate contamination over a long temporal scale, thus short period fluctuations are not relevant. This approach to the source term has been used in previous modelling works (Periáñez, 2002; Periáñez et al., 2005). Points where metals have been measured in the dissolved phase and bed sediments are shown in Fig. 9. A comparison between calculated and measured concentrations may be seen in Figs. 10–12 for Zn, Cu and Ni, respectively. The same pattern is obtained in the two phases for the three metals. This is a logical result since the bed sediment is more contaminated as concentration in water over it increases. However, it is worth mentioning that, in spite of the procedure used to define the source term, realistic concentration levels are in general simultaneously obtained for the dissolved phase and the sediment, and for the three metals. This indicates that metal partition between the dissolved and solid phases has been correctly described. Metal distribution over the Bay may be seen, in the case of Zn as an example, in Figs. 13–15, where the computed distributions in water, suspended matter and bed sediments are presented. Concentrations decrease quickly, even inside the Bay, specially in the case of the dissolved phase. This confirms the weak transport inside the Bay which was observed when simulating nitrates. Metals 0 50 100 150 200 250 0 50 100 150 200 250 X Y • 1 • 2 • 3• 4 • 5 • 6 • 7 • 8 Fig. 9. Map indicating points where metal concentrations in water and sediments have been compared with measurements. 1 2 3 4 5 6 7 8 0 5 10 15 20 25 Sampling point Cdis (µ g/l) 1 2 3 4 5 6 7 8 0 200 400 600 800 Sampling point Csed (mg/kg) Fig. 10. Computed (lines) and measured (points) Zn concentrations along the Bay’shoreline in water (top) and bed sediments (bottom). Measurement points are indicated on Fig. 9. 228 R. Periáñez / Marine Pollution Bulletin 64 (2012) 221–232 which leave the Bay are mainly directed to the east (Mediterranean Sea) due to the prevailing residual currents in the Strait (Periáñez and Caravaca, 2010). However, tidal mixing in the Strait makes a metal fraction to move to the west as well. High metal concentrations in suspended matter ð 10 2 mg=kgÞare due to the low particle concentrations (Fig. 6), but a very significant portion of the metal content, per water volume, is in dissolved form. Bed sediments are contaminated as water and suspended matter containing metals travel above them. A smooth distribution is obtained outside the Bay (Fig. 15) since sediments integrate long-scale processes as commented above. Model results, however, have to be interpreted with care (as commented before for suspended matter transport) since plumes of suspended particles and hence contaminants are sensitive to changes in wind speed and direction as well as they depend on pluviometry. Moreover, other meteorological conditions, such as atmospheric pressure differences between the Atlantic and Mediterranean, can induce flow variations through the Strait of Gibraltar. These facts make the comparison of computed metal concentrations in the water column with the corresponding measurements specially difficult (as already discussed by Dyke, 2001). However, bed sediments integrate all this variability and 1 2 3 4 5 6 7 8 0 1 2 3 4 5 Sampling point Cdis (µ g/l) 1 2 3 4 5 6 7 8 0 20 40 60 80 100 Sampling point Csed (mg/kg) Fig. 11. Computed (lines) and measured (points) Cu concentrations along the Bay’shoreline in water (top) and bed sediments (bottom). Measurement points are indicated on Fig. 9. 1 2 3 4 5 6 7 8 0 0.5 1 1.5 2 2.5 3 Sampling point Cdis (µ g/l) 1 2 3 4 5 6 7 8 0 10 20 30 40 50 60 70 Sampling point Csed (mg/kg) Fig. 12. Computed (lines) and measured (points) Ni concentrations along the Bay’shoreline in water (top) and bed sediments (bottom). Measurement points are indicated on Fig. 9. R. Periáñez / Marine Pollution Bulletin 64 (2012) 221–232 229