Environmental modelling in the Gulf of Cadiz: Heavy metal distributions in water and sediments

Environmental modelling in the Gulf of Cadiz: Heavy metal distributions in water and sediments R. Periáñez ⁎ Dpto. Física Aplicada I, E.U. Ingeniería Técnica Agrícola, Universidad de Sevilla. Ctra. Utrera km 1, 41013-Sevilla, Spain abstract Keywords: Numerical modelling Hydrodynamics Sediment transport Heavy metals Gulf of Cadiz The Gulf of Cadiz (GoC) connects the Atlantic Ocean and the Mediterranean Sea. An environmental study of the GoC is carried out through numerical modelling. First, a 3D baroclinic model is used to obtain the residual circulation and a 2D barotropic model is applied to calculate tides. The results of these models are used by a 3D sediment transport model which provides suspended matter concentrations and sedimentation rates in the GoC. Then heavy metal dispersion patterns are investigated using a 3D 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 introduced into the GoC from three rivers draining the Iberian Pyrite Belt, in the southern Iberian Peninsula. Results from the hydrodynamic, sediment and metal transport models have been compared with measurements in the GoC. In particular, the contamination of sediments collected along the southern coast of Spain is well reproduced by the model. Metal plumes reach the Strait of Gibraltar, thus the three rivers constitute a source of pollutants into the Mediterranean Sea. 1. Introduction The Gulf of Cadiz (GoC) is the sub-basin of the Atlantic Ocean which is nearest to the Strait Of Gibraltar, connecting the Atlantic Ocean and the Mediterranean Sea. Its northern, southern and eastern boundaries are, respectively, the Atlantic coasts of Spain, Morocco and the Strait of Gibraltar. The western boundary is usually defined by the 9°W meridian (Fig. 1). Recently, the GoC has been the subject of oceanographic studies dealing with surface and deep circulation, aimed at understanding the mechanisms of water exchange between the Atlantic Ocean and the Mediterranean Sea; as well as the behaviour of the dense plume of Mediterranean water (Ambar and Howe, 1979; García-Lafuente et al., 2006; Criado-Aldeanueva et al., 2006; Machín et al., 2006; García- Lafuente and Ruíz, 2007). The distribution of suspended matter and sediment transport in the GoC has also been investigated (Gonzalez et al., 2007; Freitas and Abrantes, 2002; Lobo et al., 2004; Cravo et al., 2006; Palanques et al., 1986–1987). The freshwater inputs of rivers discharging in the GoC are relatively small. However, the Guadalquivir, Guadiana and Odiel–Tinto rivers, in the southern Iberian Peninsula (Fig.1), present strongly enhanced heavy metal concentrations since they drain the Iberian Pyrite Belt (Sainz and Ruiz, 2006), one of the most important mining areas in the south of Europe. Mineral resources have been extracted in the last 5000 years during two main periods: the Roman age and the last two centuries. During the last period, intensive exploitation has led to a relevant environmental impact, with vast surfaces covered with mining residues and subjected to erosion (Sainz and Ruiz, 2006). The GoC is responsible for 5–10% of fish and shell-fish catches of Spain and Portugal (Beckers et al., 2007), holding important living resources of commercial and ecological interest. Consequently, it is relevant to study and understand the geochemistry and dispersion patterns of heavy metals in the GoC system, since this will help assessing the potential influence of metals on ecosystem functioning. Thus, several papers concerning the distribution of metals in the GoC have been published in the last years (Sainz and Ruiz, 2006; Morillo et al., 2004; Elbaz-Poulichet et al., 2001; Fernández- Caliani et al.,1997). Indeed, studies on the distribution of trace metals in coastal waters are frequently published in recent times (see for instance Li et al., 2007; Suntornvongsagula et al., 2007; Cuong et al., 2008; Chen and Jiao, 2008; Valdés et al., 2008; Marín-Guirao et al., 2008). The objective of this work consists of studying the dynamics of heavy metals in the northern GoC by means of numerical modelling. Models have been widely applied to simulate contaminant dispersion since they may provide insights on the main environmental processes governing such dispersion and, consequently, may help to describe and characterize the environment (Scott, 2003). In particular, models have been widely applied to heavy metal (Tappin et al.,1997; Wu et al., 2005; Prandle et al., 1996; Liu et al., 1998) and radioactive element (Harms, 1997; Cetina et al., 2000; Periáñez, 2003; Monte et al., 2006) dispersion in coastal waters. ⁎Fax: +34 954486436. E-mail address: [email protected]. Fig. 1. General localization of the study area and topography of the GoC (depths in m). The localization of Guadiana, Odiel–Tinto and Guadalquivir rivers is also shown. Although some interesting modelling works describing water circulation off Iberia and Morocco coasts have been published (Johnson and Stevens, 2000; Batteen et al., 2000), these models have a relatively low resolution, not providing detailed information about the GoC basin circulation features. In the first case, spatial resolution is 10 min in both longitude and latitude (5 times the grid cell in the present model). In the second reference it is 10 km in longitude and latitude, about three times larger than in this work. Other modelling works are specifically devoted to the study of Mediterranean water spreading (Jungclaus and Mellor, 2000; Johnson et al., 2002; Serra et al., 2005). The excellent paper by Peliz et al. (2007) presents three nested model domains aimed at reproducing known features of the Azores current and of circulation inside the GoC. Published models describing trace metal dispersion in the GoC consider metals as conservative tracers, no interacting with sediments and without any other sources and sinks (Elbaz-Poulichet et al., 2001; Beckers et al., 2007). The first authors use a model to estimate the dilution of a conservative tracer released by the Odiel–Tinto rivers. Beckers et al. (2007) apply a numerical model to reproduce observed metal (again considering metals as conservative tracers) concentrations in surface waters of the northern GoC. Analysis of model results showed that sources/sinks of metals due to interactions with sediments (adsorption/desorption reactions as well as erosion and deposition processes) were apparent. Models which try to reproduce measured levels of metals in bed sediments of the GoC have not been published yet, to the author's knowledge. The model described in this paper consists of three sub-models: firstly, a hydrodynamic module which provides currents over the domain. Two hydrodynamic models are used. A 2D barotropic model is applied to calculate tides and 3D baroclinic model is used to obtain the residual (mean) circulation. Tidal currents must be calculated since they may increase the bed stress and hence enhance sediment resuspension and affect deposition of particles as well. Indeed, the second sub-model is a sediment transport model which provides suspended matter concentrations and sedimentation rates over the domain. The third sub-model is the metal transport module, which includes advection/diffusion plus uptake/release reactions of metals between the dissolved and solid (suspended matter in the water column and bed sediments) phases. The three modules, as well as the numerical techniques used to solve the involved equations are described in Section 2. Next, model results are discussed separately for water circulation, sediment transport and metal distributions. Some sensitivity analysis are finally described. 2. Model description 2.1. Hydrodynamic models As commented above, tides are required to calculate bed stress over the domain, since it will affect sedimentation rates in the shallower areas and/or where stronger tidal currents exist. A 2D depth-averaged model has been applied to calculate surface tides. This is a reasonable approach that has already been successfully used in the Strait of Gibraltar (Tejedor et al.,1999; Periáñez and Pascual-Granged, 2008), the Alborán Sea (Periáñez, 2008) and even in the complete Mediterranean Sea (Tsimplis et al.,1995).As has been shown, it is safe to neglect density differences in tidal computations (Dyke, 2001; Yanagi, 1999). Equations are standard and may be seen for instance in Periáñez (2005a). 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)andtidal constants are stored in files that will be read by the sediment transport model. The barotropic 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. The long-term circulation is obtained from a full 3D, primitive equation, baroclinic hydrodynamic model. It is based upon the hydrostatic and Boussinesq approximations on a βplane. The model includes equations for salinity and temperature evolution and water density is calculated from them using a standard state equation (Batteen et al., 2000). The one-equation turbulence model described in Periáñez (2005b) has been used to calculate the vertical eddy viscosity. Details on the 3D equations may be seen, for instance, in Kowalick and Murty (1993). A summary of the main equations involved in both hydrodynamic models is presented in Appendix A. 2.2. Sediment transport The transport of sediments is described by a 3D advection–diffusion equation to which some terms are added. These are external sources of particles (river supply), terms describing particle deposition on the seabed and erosion from the bed to the water column, and vertical settling. The formulation of these processes is based upon standard formulae. Thus, the erodability constant is used for the erosion term. Particle settling and deposition are described using the settling velocity, which is obtained from Stoke's law. Critical erosion and deposition stresses are used as usual. Details on the mathematical formulation may be seen elsewhere (Periáñez, 2005b; Liu et al., 2002a; Lumborg and Windelin, 2003; Cancino and Neves, 1999). Finally, it is also possible to calculate net sedimentation rates (SR) as the balance between the deposition and erosion terms. A summary of the equations may be seen in Appendix B. 2.3. Metal 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 (suspended matter and bed sediments). Thus, three phases are included in the model: dissolved, suspended particles and bed sediments. Only fine sediments are considered (particles with diameter b62.5 µm) since metals are predominantly fixed to these (Eisma, 1993). 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 metals to the deep sediment is included. Thus, metals 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. This effect may be easily treated as a decay process with constant λ burial given by (Periáñez, 2008): λburial =SR ρsLð1Þ where Lis the sediment mixing depth (the distance to which the dissolved phase penetrates the sediment), ρ s is the sediment bulk density (dry mass divided by wet volume) and SR is the sedimentation rate provided by the sediment transport sub-model. It is known that adsorption depends on the surface of particles per water volume unit at each point and time. This quantity has been denoted as the exchange surface (Periáñez, 1999, 2000, 2002). Thus, the kinetic coefficient k 1 is written as: k1=χSm+Ssδb ðÞ=km 1+ks 1ð2Þ where S m and S s are the exchange surfaces for suspended matter and bottom sediments respectively (dimensions [L] −1 )andχis a parameter with the dimensions of a velocity denoted as the exchange velocity (Periáñez, 1999, 2000, 2002). The delta function is introduced to take into account that only the deepest water layer interacts with the bed sediment. Thus δ b =1 for the deepest layer and δ b =0 elsewhere. Assuming spherical particles, the exchange surfaces are written as (see references cited above): Sm=3m ρRð3Þ and Ss=3Lf1−pðÞ Rψð4Þ where Rand ρare particle radius and density respectively, mis the suspended matter concentration, fis the fraction of fine particles in the sediment, 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. Finally, ψis the thickness of the deepest water layer. 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 kinetic coefficient k 2 is considered to be constant. The equation that gives the time evolution of metal concentration in the dissolved phase, C d , is: @Cd @t+u@Cd @x+v@Cd @y+w@Cd @z=A@2Cd @x2+@2Cd @y2! +@ @zK@Cd @z  −k1Cd +k2mCs+δbk2 LρsfAs ψ ð5Þ where C s and A s are, respectively, the concentrations of metals in suspended matter and bottom sediments. u,vand ware water velocities along the x,yand zaxis and Aand Kare, respectively, the horizontal and vertical diffusion coefficients. The equation that gives the time evolution of metal concentration in suspended matter is: @mCs ðÞ @t+u@mCs ðÞ @x+v@mCs ðÞ @y+w−ws ðÞ @mCs ðÞ @z=A@2mCs ðÞ @x2+@2mCs ðÞ @y2! +@ @zK@mCs ðÞ @z  +km 1Cd−k2mCs−δbSED ð6Þ where w s is the particle settling velocity and SED is the deposition of metals from the deepest water layer to the sediment evaluated according to: SED =SRCsbðÞ ψ:ð7Þ Note that (b) means that the corresponding magnitude is evaluated at the deepest water layer. The equation for the temporal evolution of metal concentration in the bottom sediment mixed layer is: @As @t=ks 1 CdbðÞψ Lρsf −k2As−λburialAs+SED ð8Þ where the deposition is now calculated as: SED =SRCsbðÞ ρsLf :ð9Þ The total metal content, A p , in the sediment below the mixed depth is given by the following equation: @Ap @t=λburialLρsfAs:ð10Þ 2.4. Numerical solution and parameters The topography of the GoC was obtained from the NOAA (US National Ocean and Atmosphere Administration) GEODAS database, available on-line, with a resolution of 2 min in both longitude and latitude. The grid consists of 64×136 cells, extending from −8° to −5.9° in longitude and from 33° to 37.5° in latitude. Fifty vertical levels are considered, with increasing thickness from the sea surface to the bottom. All the equations were solved using explicit finite difference schemes. In particular, the scheme described in Periáñez (2005a) is used in the 2D barotropic hydrodynamic model. In the 3D baroclinic model the scheme in Kowalick and Murty (1993) is applied. The second order accuracy MSOU (Monotonic Second Order Upstream) advection scheme (Vested et al.,1996) is applied for all advective terms and a second order scheme is used for diffusion as well (Kowalick and Murty, 1993). Time-step in the different sub-models is selected with care to assure that stability conditions are satisfied. A time-splitting procedure was used to treat the water–sediment metal interactions. This is essentially the same as in hydrodynamic models, when time-splitting is used to solve the external and internal modes. The procedure is required because a very restrictive stability condition is introduced by such interaction terms (Periáñez, 2005a). Essentially, it consists of using a smaller time step to solve water–sediment metal interactions. Thus, for each advection/diffusion time step, a number of smaller steps are needed to solve metal partition. Some open boundary conditions must be provided. In the barotropic model, tidal surface elevations and phases are specified along the open boundaries of the computational domain from Schwiderski (1980a,b). A radiation condition is applied to the water velocity component that is normal to the open boundary (Jensen, 1998). The equations are integrated, starting from rest, until stable oscillations are obtained. Then tidal analysis is carried out and tidal constants are stored in files. In the case of the 3D baroclinic model, water velocities are specified along the Strait of Gibraltar from the average inflow (into the Mediterranean) of Atlantic Water and deep outflow (from the Mediterranean to the Atlantic) of Mediterranean water. Water temperature and salinity along the open boundaries are specified from the NOAA WOA05 (World Ocean Atlas 2005), available on-line. Summer values are used since sampling campaigns were carried out in this season. The 3D model is again started from rest and a wind stress of 0.045 N/m 2 , directed to the south, is applied. It corresponds to the summer characteristic value (Machín et al., 2006). Residual currents are stored in files that will be later read by the sediment and metal transport sub-models. A summary of parameters involved in the models is given in Table 1. Suspended particles (diameter b63 µm) are characterized by a mean size that is considered to be representative of suspended matter in this environment. The particle size controls, through Stoke's law, the settling velocity of particles and, as a consequence, affects the sedimentation rate as well. It also affects the adsorption of metals from the dissolved phase (Eqs. (3) and (4)). The value used for particle diameter was 8 µm. Freitas and Abrantes (2002) have found that particles with diameter b10 µm are dominant in all water masses in the Gulf of Cadiz. Moreover, this is the same value selected for a model of the Alborán Sea (Periáñez, 2008) after a sensitivity analysis (the GoC reflects the Alborán Sea since surface waters in the first flow into the second and Mediterranean waters flow, in the deep layer, from the Alborán Sea into the GoC). Finally, the selected value is also consistent with Liu et al. (2002b), who used a radius of 2.5 µm to describe fine silts in their model. Of course, it would be convenient to have a detailed particle size spectra and then using several particle classes (Periáñez, 2005b), but it is not available in current literature. The critical deposition stress typically ranges between 0.04 and 0.1 N/m 2 (Tattersall et al., 2003), while the critical erosion stress ranges 0.1–1.5 N/m 2 (Tattersall et al., 2003). In the present application intermediate values of 0.06 and 1.0 N/m 2 have been taken for the critical deposition and erosion stresses respectively, as in the Alborán Sea model described in Periáñez (2008). 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). The fraction fof fine particles in bed sediments must be defined over the model domain. This information has been compiled from Hernández- Molina et al. (2006) and Gonzalez et al. (2007). Schematically, sand and sandy mud dominate to a depth of 30 m. The middle shelf, to a depth of 100 m, is characterized by an extensive mud belt. Finally, outer sediments below 100 m are dominated by mixtures of sand and clay. Suspended matter concentrations in the three main rivers (Fig. 1) have been defined. Of course, these magnitudes present seasonal (and at a shorter scale as well) variations. However, representative mean values have been used for the present application. In the Odiel–Tinto estuary, the mean value measured in Periáñez et al. (1996) has been used: 30 g/m 3 . The mean value measured by Gómez-Parra et al. (2000), 120 g/m 3 , was selected for the Gualdalquivir estuary. Finally, in the Guadiana estuary, a mean value of 30 g/m 3 is in agreement with Machado et al. (2007) and with the empirical relation given for this estuary by Wolanski et al. (2006). The sediment transport sub-model is run until steady particle distribution and SR are obtained. The model is forced by residual currents provided by the 3D baroclinic model. Instantaneous tidal currents are used only to calculate the bed stress. This approach has been made since the model has been run tosimulate the dispersion of a passive tracer with and without tides (of course residual circulation is considered in both cases), and differences in the tracer concentration fields between both simulations could not be appreciated (the simulated time was just a few months). Suspended particle concentrations and sedimentation rates are stored in files that are read by the metal transport model. The vertical diffusion coefficient is taken equal to the vertical eddy viscosity, which is provided by the turbulence equation included in the 3D baroclinic model. The horizontal diffusion coefficient depends on the horizontal grid spacing. Following Dick and Schonfeld (1996): A=0:2055 × 10−3Δx1:15:ð11Þ The present grid resolution (2 min) gives A=2.0 m 2 /s. The model has been applied to simulate the dispersion of three metals: Zn, Cu and Ni. Values for parameters describing metal transport are summarized in Table 1 too. The correction factor ϕhas been fixed as 0.1 according to previous modelling works (Periáñez, 2000, 2003; Periáñez et al., 2005). An intermediate value of 0.5 was used for the sediment porosity. Indeed, it is 0.6 in the Alborán Sea model according to measurements (Periáñez, 2008). The sediment mixed depth is L=0.1m(Periáñez, 2000, 2003). The so-called equilibrium distribution coefficient, k d , describes the partition of a tracer between the dissolved and solid phases at equilibrium (Duursma and Carroll, 1996): kd=Cs Cd ð12Þ where C s and C d are tracer concentrations in the solid and dissolved phases respectively. This parameter is of course depending on the Table 2 Values of river discharge metal concentrations (nM). Metal Guadiana Guadalquivir Odiel–Tinto Cu 42 58 214 Ni 20 43 67 Zn 49 33 129 Table 1 Summary of model parameters. Parameter description Value Source Water kinematic viscosity ν=1.064×10 −6 m 2 /s Standard value Vertical diffusion coefficient K=variable Turbulence equation Horizontal diffusion A=2.0 m 2 /s Dick and Schonfeld (1996) Sediment mixing depth L=0.1 m Periáñez (2000, 2003) Particle density ρ=2600 kg/m 3 Standard value Particle radius R=4.0 µm Freitas and Abrantes (2002) Fraction of small particles in sediment f=variable Gonzalez et al. (2007) Correction factor ϕ=0.1 Periáñez (2000, 2003); Periáñez et al. (2005) Sediment porosity p=0.5 Periáñez (2008) (similar value) Erodability E=1.6×10 −3 kg/m 2 sTattersall et al. (2003) Critical erosion stress τ ce = 1.0 N/m 2 Tattersall et al. (2003) Critical deposition stress τ cd =0.06 N/m 2 Tattersall et al. (2003) Desorption kinetic coefficient k 2 =1.16×10 −5 s −1 Nyffeler et al. (1984) Zn distribution coefficient k d =7.0×10 4 IAEA (2004) Ni distribution coefficient k d =2.0×10 4 IAEA (2004) Cu distribution coefficient k d =2.6×10 4 Calibration In the cases of Eand the critical stresses for erosion and deposition, the selected values correspond to intermediate values within their range of variation commonly found in literature. A reference (Tattersall et al., 2003) is given as an example. Fig. 2. Map of the GoC showing all locations mentioned in the text. Table 3 Established, index obs,(NOAA, 1982) and computed, index comp, amplitudes (A, cm) and phases (g, deg) of tidal elevations at several locations indicated in Fig. 2. Station M 2 S 2 A obs g obs A comp g comp A obs g obs A comp g comp Faro 92 94 99 68 32 125 36 91 Chipiona 102 54 104 62 41 82 38 85 Rota 105 52 103 62 37 78 38 85 Cadiz 100 87 99 61 37 110 36 83 Ayamonte 100 59 101 65 32 88 36 89 Huelva 102 56 105 65 38 82 38 88 Casablanca 99 56 92 53 35 81 36 77 Rabat 88 59 98 57 35 83 36 78 tracer geochemical behaviour and on environmental conditions, thus it is a site-specific parameter. As described before (Nyffeler et al., 1984), k 2 is very similar even for elements with a rather different geochemical behaviour, being χ (see Eq. (2)) the essential parameter describing the tracer geochemical behaviour. Thus, the same value is given to k 2 for Ni, Zn and Cu, which is k 2 =1.16×10 −5 s −1 . This value has been successfully used in earlier simulations in the Odiel–Tinto estuary (Periáñez et al., 2005) and in the Strait of Gibraltar–Alborán Sea (Periáñez, 2008) for several elements (Ra, Cs, Pu). It was measured for Cs by Nyffeler et al. (1984). The exchange velocity for each metal is obtained from the equation relating this parameter with k 2 and its corresponding equilibrium distribution coefficient, k d (Periáñez, 2005a): kd=χ k2 3 ρR:ð13Þ The k d mean values recommended by the International Atomic Energy Agency (IAEA, 2004) for coastal waters have been used. This procedure has been widely used before (Periáñez 1999, 2000, 2003, 2008). Distribution coefficients are given in Table 1. Nevertheless, they may vary in more than one order of magnitude depending on environmental conditions. Consequently, the model sensitivity to this parameter has been studied and results are presented below (Section 3.3.2). The IAEA does not provide k d values for Cu, thus in the case of this metal it had to be obtained from a calibration procedure. The selected value, 2.6×10 4 is, however, of the same order of magnitude of the measured Cu k d (Wen et al., 1999) in USA coastal waters (1.3×10 4 ). Metal concentrations in the dissolved phase have been defined at the three main sources, Guadalquivir, Guadiana and Odiel–Tinto estuaries. Values, which are the same used by Beckers et al. (2007), are given in Table 2. The metal transport sub-model is started from background dissolved concentrations, which correspond to the open Atlantic Ocean metal concentrations reported in Elbaz-Poulichet et al. (2001). The corresponding background in the solid phase is obtained from the metal k d through Eq. (12). Equations are integrated until steady metal distributions are obtained in all phases. All the modelling work has been carried out in a sequential way. Thus, the hydrodynamic sub-model is developed first. Once that acceptable tides and residual currents are obtained, these are used to simulate sediment dynamics. Finally, metal transport is included when the sediment sub-model model has been validated. 3. Results and discussion A brief description of the hydrodynamic and sediment transport results is given. Next, model results for metal distributions are discussed. Fig. 3. Amplitude of the M 2 tide in m (left) and current amplitude in m/s (right). Fig. 4. Computed surface residual currents. Only one of each 4 vectors is drawn. 3.1. Hydrodynamics A map showing the different locations and sampling points mentioned in the paper is presented in Fig. 2. Computed tidal constants have been compared with established values for both M 2 and S 2 tides. Results are shown in Table 3, and it may be seen that, generally, there is a good agreement between both set of data. As an example, maps showing the M 2 tide amplitude and the amplitude of the tidal current are presented in Fig. 3. The amplitude of the tide is about 1 m over all the GoC, decreasing near the Strait of Gibraltar. Associated currents are weak, with amplitudes below 0.10 m/s over most of the GoC. Indeed at the RAP (Red de Aguas Profundas, Spanish Institute of Oceanography) buoy position (see Fig. 2), the M 2 barotropic tidal current is less than 0.03 m/s (García-Lafuente et al., 2006). The computed current at this position is 0.034 m/s. The current amplitude increases as approaching the Strait entrance, where currents about 0.8 m/s are produced (the color scale in Fig. 3 is limited to 0.5 m/s for more clarity). A similar behaviour is observed for the S 2 tide. The residual surface circulation in the northern GoC is characterized by a current directed to the SE (Criado-Aldeanueva et al., 2006; García-Lafuente and Ruiz, 2007) along the Spanish coast. This circulation is a rather constant pattern during most of the year. Part of the flow enters the Strait of Gibraltar and part is deflected to the south. The residual circulation computed with the baroclinic model at the sea surface is presented in Fig. 4. The current is effectively directed to the SE over the continental shelf and part of this flow enters the Strait of Gibraltar. Maximum currents are of the order of 0.3 m/s, in agreement with García-Lafuente et al. (2006). The anticyclonic eddy at the east of Faro (see Fig. 2) has been described by Machín et al. (2006). Also, the cyclonic eddy in front of the Strait of Gibraltar appears clearly in the models of Beckers et al. (2007) and Peliz et al. (2007). These last authors have attributed it to the strong convergence occurring in this area. Below the surface, the Mediterranean waters flow into the Atlantic and mainly direct to the NW (Serra et al., 2005; Ambar and Howe, 1979; Johnson et al., 2002). As an example, the computed circulation 590 m below the surface is presented in Fig. 5. Only the northern part of the GoC is shown to appreciate details more clearly. These currents are in agreement with the geostrophic velocities below 400 m, referenced to 300 m, provided by Criado-Aldeanueva et al. (2006) and with the calculations by Peliz et al. (2007) for summer. Water velocity is higher close to the Strait and then slows to about 0.1 m/s, in agreement with Ambar and Howe (1979). Temperature and salinity profiles in the water column calculated by the model have been compared with those obtained from observations carried out under the TOROS project (Elbaz-Poulichet et al., 2001) in summer 1997. Examples for 4 points shown in Fig. 2 may be seen in Fig. 6. There is an acceptable agreement between model results and experimentally obtained vertical profiles of Tand S. Finally, two vertical south–north sections of water salinity are presented in Fig. 7, at longitudes of −6.87° and −6.05°. In the first Fig. 5. Computed residual currents 590 m below the surface in the northern GoC. Only one of each 4 vectors is drawn. Fig. 6. Computed (lines) and measured (points) salinity and temperature profiles in points indicated in Fig. 2. case, the core of more dense Mediterranean water is clearly seen, with the salinity maximum below 600 m. This flow is aligned with the Spanish continental slope, in agreement with Jungclaus and Mellor (2000) and Peliz et al. (2007). Close to the Strait, the situation may be characterized as a two-layer exchange flow, with the interface tilted down southward (Jungclaus and Mellor, 2000). Although seasonally averaged values have been used for open boundary conditions and wind stress, the main features of water circulation in the GoC are, generally speaking, reproduced by the model for the considered season. 3.2. Suspended matter As expected, the river-discharged suspended matter is transported to the southeast along the Spanish coast, which is the residual current direction in this area of the GoC. As widely discussed (Gonzalez et al., 2007), there is a dominant eastward transport throughout the entire northern GoC. Some authors (Lobo et al., 2004) have postulated the existence of some westward transport of sediments released by the Guadalquivir River in the inner shelf, although little indication of it has been found in other works (Gonzalez et al., 2007). Maps of suspended matter concentrations at the surface and in the deepest water layer are presented in Fig. 8. Logarithmic scale is used to appreciate differences. The suspended particle plumes produced by the three rivers can be clearly seen in the maps. Indeed, concentrations of the order of 10 g/m 3 are obtained near the river mouths. Concentrations of the order of 10 −1 g/m 3 are obtained in part of the northern GoC and much smaller values are apparent to the south. Computed surface particle concentrations are in agreement, by order of magnitude, with measured surface concentrations. Indeed, concentrations of the order of 10 −1 g/m 3 have been measured at the surface in the northern GoC (Palanques et al., 1986–1987). Also, concentrations in the range 20–45 g/m 3 have been detected in the Guadiana River plume (Cravo et al., 2006) over a narrow area (b10 km from the coast). In the deepest water layer concentrations are slightly higher than at the surface over most of the northern GoC, in agreement with Palanques et al. (1986–1987). In some areas as the Strait of Gibraltar and the Spanish continental slope, concentrations are more significantly enhanced in the deep water. This is probably due to erosion produced by stronger tidal currents (in the area of the Strait) and also produced by the Mediterranean water current. A south–north section of suspended matter concentrations at the Atlantic entrance of the Strait of Gibraltar (longitude −6.10°) is shown in Fig. 9. The plume of particles released by the rivers may be seen in surface waters at the northern area of the section. Maximum concentrations, however, are obtained in Mediterranean waters, at depth about 200 m and along the Spanish continental slope. This indicates that erosion is occurring in this area, as commented above. Computed net sedimentation rates (not shown) are larger, as expected, in the vicinity of the river mouths, where they are of the order of 10 −3 g/cm 2 year. Much smaller values are obtained far from the river discharge influence. In particular, there is a region where no sedimentation occurs at the west of the Strait of Gibraltar. This is due to the strong currents in the Mediterranean outflow water (Fig. 5), which keep particles in suspension. Nevertheless, all these results have to be interpreted with care since several approximations have implicitly been made in the model: there are variations in water exchanges through the Strait of Gibraltar, as already discussed in previous modelling works (Periáñez, 2008). Of course, these variations will affect the SR. Although constant suspended matter concentrations have been defined in the three estuaries, there will also be seasonal variations depending on pluviometry, for instance. Only fine particles are considered in the model, and bed load transport is neglected. This approximation will lead to lower sedimentation rates. The same occurs due to the fact that organic particles are not included. Finally, atmospheric deposition events of particles coming from the Sahara Desert (Fabres et al., 2002) have not been considered in the model since they cannot be easily quantified (Periáñez, 2008). It must be taken into account that suspended particles will affect metal transport only close to the coast, given the extremely low suspended matter concentrations which are measured offshore (Palanques et al., 1986– 1987), and river supply is the main sediment source in this coastal area. 3.3. Heavy metal transport 3.3.1. Metal distributions As commented before, the model has been applied to simulate the transport of three metals in the GoC: Zn, Cu and Ni. Pathways of metals are essentially the same as for suspended particles. Thus, theyare transported to the east, along the Spanish coast, until the Strait of Gibraltar. The concentrations of several metals have been measured in the fine sediment fraction (b63 µm) along the Spanish coast from the Fig. 7. Vertical south–north salinity sections at longitudes −6.87° and −6.05° (left and right respectively). Guadiana to the Guadalquivir mouths (Sainz and Ruiz, 2006; Morillo et al., 2004) and also at some points closer to the Strait of Gibraltar (Riba et al., 2002). Samples were collected at an approximate distance of 500 m from the shoreline. A comparison of measured (points) and computed (lines) metal concentrations along the northern coast of the model domain can be seen in Fig. 10. In general, there is a good agreement between the measured and calculated concentrations, although in the case of Ni the model produces a concentration peak in the area of the Odiel–Tinto mouth that is not apparent in measurements. This discrepancy could be caused, at least partially, by a different nature of particles released from this river, which affects the partition on Ni (although does not in the cases of Zn and Cu). Also, it is possible that the source term is overestimated in the case of Ni for the Odiel–Tinto River. A clear conclusion cannot be obtained. Metal concentrations are very low westward from the Guadiana River. There is an increase in concentrations here since, as has already been mentioned,the three rivers drain the IberianPyrite Belt. Maximum concentrations exist in the mouth of the Odiel–Tinto rivers. Although river flows are much smaller than those of the Guadiana and Guadalquivir, Odiel–Tinto rivers are considerably more contaminated (Gonzalez et al., 2007) and, indeed, they have been recognized as the main source of metals along the coast (Morillo et al., 2004). There is not an appreciable metal enhancement in the area of the Guadalquivir mouth (except in the case of Ni where both measurements and calculations show such an enhancement). Metal concentrations at the east of the Odiel–Tinto rivers decrease much slower than at the west of the Guadiana River, which is due to the residual currents in the shelf, which flow towards the Strait of Gibraltar. Indeed, it has already been found (Elbaz-Poulichet et al., 2001) that coastal waters transport dissolved metals (which contaminate bed sediments as metals travel over them) from the Odiel–Tinto rivers to a Fig. 8. Computed surface suspended matter concentrations (g/m 3 ) in logarithmic scale at the surface (top) and bottom (down) of the water column. Fig. 9. Vertical south–north section of suspended matter concentrations (g/m 3 ) at the Atlantic entrance of the Strait (longitude −6.10°). Fig. 10. Measured (circles) and computed (lines) metal concentrations in sediments along the northern coast of the model domain. The three arrows indicate the location of the Guadiana, Odiel–Tinto and Guadalquivir mouths (from left to right).