Tidal circulation, sediment and pollutant transport in Cádiz Bay (SW Spain): A modelling study

Tidal circulation, sediment and pollutant transport in Cádiz Bay (SW Spain): A modelling study R. Periáñez a, n , M. Casas-Ruíz b , J.P. Bolívar c a Dpt Física Aplicada I, ETSIA, Universidad de Sevilla, Ctra Utrera km 1, 41013-Sevilla, Spain b Dpt Física Aplicada, Universidad de Cádiz, 11510-Puerto Real, Cádiz, Spain c Dpt Física Aplicada, EPS La Rábida, Universidad de Huelva, Ctra Palos s/n, 21819-Huelva, Spain article info Article history: Received 5 November 2012 Accepted 12 May 2013 Available online 10 June 2013 Keywords: Cádiz Bay Numerical modelling Hydrodynamics Sediment Metals Radionuclides abstract A numerical model to simulate the dispersion of particle-reactive tracers in Cádiz Bay (SW Spain) has been developed. It includes a hydrodynamic submodel to provide water currents, a sediment transport submodel, which provides suspended matter concentrations and sedimentation rates and the pollutant dispersion model. Pollutant exchanges between the liquid and solid phases are described in a dynamic way, using kinetic transfer coefficients. Results of the hydrodynamic and sediment transport models have been compared with observations. The contamination of sediments of the bay by fallout 137 Cs has been simulated. The existence of an accumulation area in the inner bay, found in field measurements, has been reproduced by the model. Dispersion experiments for Zn and 226 Ra have also been carried out. Flushingtimes of the inner bay have been finally determined through numerical experiments. &2013 Elsevier Ltd. All rights reserved. 1. Introduction Cádiz Bay is located in the southwest of Spain. It is a semienclosed water body, landlocked except along its northwestern limit where it faces the Atlantic Ocean. It may be divided (Fig. 1) into two regions: the shallower inner bay and the deeper outer bay, open to the Atlantic Ocean, which are connected through a narrow channel. Anyway the whole bay is shallow, with maximum depth around 20 m in its connection to the open sea. The inner bay is characterized by the presence of extensive intertidal areas, which are covered by water during part of each tidal cycle. Tides are predominantly semidiurnal, with a M 2 amplitude of approximately 1 m. Water circulation inside the bay is mainly controlled by tidal dynamics (Araujo et al., 2009). Several cities surround the bay: Rota, El Puerto de Santa María, Puerto Real, San Fernando and Cádiz; thus some 400 000 people live in the area. Main industrial activities are related to ship, offshore, car and aerospacing manufacturing (Carrasco et al., 2003). As a result of this population and industrial activity contamination in water and seabed sediments of the bay has been reported (Araujo et al., 2009; Carrasco et al., 2003;Gómez-Parra et al., 1984;Ligero et al., 2002). Some interesting modelling works concerning tidal dynamics in the bay have already been carried out (Alvarez et al., 1999,2003; Kagan et al., 2001;Alvarez, 1999). However, modelling studies on the dynamics of pollutants in Cádiz Bay have not been published before. These studies are useful, through the knowledge of water and sediment dynamics, to predict the effects of planned authorized releases of pollutants, as well as to assess the impact of potential accidental spills in the system. Studies of similar nature have been recently reported for Algeciras Bay (south Spain), also affected by industrial and urban waste releases (Periáñez, 2012), and for southern Spain coastal waters (Periáñez, 2011). Thus, the objective of this paper consists of presenting a model for the Cádiz Bay suitable for the description of the transport of non-conservative pollutants, which are those that do not remain dissolved in the water column, but are fixed to suspended particles and to seabed sediments. Actually, many pollutants as heavy metals and some radioisotopes are predominantly attached to these solid phases. A hydrodynamic model, which provides the currents that induce advective transport and turbulent mixing in the water column, as well as a sediment transport model is required for describing the transport processes of these particle-reactive pollutants. The sediment model will calculate suspended matter concentrations and sedimentation rates, which are relevant for describing adsorption of pollutants by suspended matter and pollutant exchanges between suspended particles in the water column and seabed sediments. The model is described in the next section. Later, results are presented and discussed. 2. Model description 2.1. Tidal simulations As mentioned above, water circulation in the bay is essentially controlled by tides, thus our focus has been simulating tidal dynamics. Contents lists available at SciVerse ScienceDirect journal homepage: www.elsevier.com/locate/oceaneng Ocean Engineering 0029-8018/$ - see front matter &2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.oceaneng.2013.05.016 n Corresponding author. Tel.: +34 954486474. E-mail address: [email protected] (R. Periáñez). Ocean Engineering 69 (2013) 60–69 Some authors have stated that using a 2D model for simulating tides is areasonableapproach(Dyke, 2001;Yanagi, 1999).Also,itmustbe considered that the bay is rather shallow. For instance, mean depth in the inner bay is only about 3 m. Also, flows of channels and rivers discharging in the bay are extremely low and have been neglected in previous modelling studies of tidal dynamics (Alvarez, 1999). Significant horizontal or vertical density gradients are not expected. Consequently, depth-averaged models have already been successfully used to study the hydrodynamics of Cádiz Bay (Alvarez, 1999;Alvarez et al., 2003;Kagan et al., 2001). The influence of sediment load on tidal dynamics has also been studied with a depth-averaged model (Alvarez et al., 1999). Moreover, the theoretical analysis in Prandle et al. (1993) indicates that vertical homogeneity of substance concentrations is maintained for water depths lower than 50 m if tidal action ensures that the vertical eddy viscosity is at least of the order of 10 −3 m 2 /s, a value which may be considered representative of a relatively weak tide. Consequently, the use of a 2D depth averaged barotropic model is justified as a first approach. It is based on the following equations: ∂ζ ∂tþ∂ ∂xðHuÞþ ∂ ∂yðHvÞ¼0;ð1Þ ∂u ∂tþu∂u ∂xþv∂u ∂yþg∂ζ ∂x−fv þτ u ρ w H¼A∂ 2 u ∂x 2 þ∂ 2 u ∂y 2  ;ð2Þ ∂v ∂tþu∂v ∂xþv∂v ∂yþg∂ζ ∂yþfu þτ v ρ w H¼A∂ 2 v ∂x 2 þ∂ 2 v ∂y 2  ;ð3Þ where uand vare the depth averaged water velocities along the x-andy-axes, his the depth of water below the mean sea level, ζis the displacement of the water surface above the mean sea level (positive upwards), H¼hþζis the total water depth, fis the Coriolis parameter (f¼2Ωsin λ,whereΩis the Earth rotational angular velocity and λis latitude), gis acceleration due to gravity, ρ w is water density and Ais the horizontal eddy viscosity. Friction stresses τ u and τ v are written in terms of a quadratic law: τ u ¼kρ w uffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 þv 2 p τ v ¼kρ w vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 þv 2 p;ð4Þ 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 used by the dispersion model to calculate the advective transport. The model includes the following tidal constituents, M 2 ,S 2 ,N 2 ,M 1 ,O 1 and M 4 . 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. Tidal residuals for each constituent are calculated from the equation: q r !¼〈Hq t !〉 〈H〉;ð5Þ which corresponds to the Eulerian residual transport velocity (Delhez, 1996). In this equation 〈〉 is the time averaging operator, q r !is the tidal residual and q t !is the instantaneous tidal current. Tidal constants and residual are stored in files which are later read by the sediment and pollutant transport models for a extremely fast calculation of current in each grid cell in the domain. Due to the large intertidal areas existing in the inner bay, a flood-dry algorithm is required. The one reported in Kampf (2009) has been applied in this study. Wet grid cells are defined as those with a total water depth Hlarger than a threshold value H min typically set as a few centimeters. Dry cells are defined as cells where H≤H min . Flooding and drying is implemented in the code via the calculation of the water velocity normal to the interface between wet and dry cells. The calculation is performed when the pressure gradient force is directed towards the dry cell. Otherwise velocity is set to zero at this point. In the case of a non-zero velocity, water level in the dry cell will increase and the cell turns into a wet one once the water depth is larger than H min . 2.2. Sediment transport Sediment dynamics is described in a very simple way, but still retaining its main aspects. This has proven to be enough to satisfactorily simulate reactive pollutant dispersion. 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. The formulation of these processes is based upon standard formulae. Thus, the erodability constant is used for theerosionterm.Particledeposition is described using the settling velocity, which is obtained from Stokes'law. Critical erosion and deposition stresses are applied as usually (Periáñez, 2005a,2005b; Liu et al., 2002;Lumborg and Windelin, 2003). The equation for suspended sediment transport is ∂ðHmÞ ∂tþ∂ðuHmÞ ∂xþ∂ðvHmÞ ∂y¼∂ ∂xHK h ∂m ∂x  þ∂ ∂yHK h ∂m ∂y  þðER−DEPÞþS;ð6Þ where mis the suspended matter concentration, K h is an effective horizontal diffusion coefficient, Sis the external particle source and ER and DEP are the erosion and deposition terms, respectively. 10 20 30 40 10 20 30 40 50 60 70 80 90 100 110 X Y • Carraca • Pto Real • Pte Carranza Pto Cadiz • • Pto Sherry • Rota Outer BayInner Bay ←Guadalete River ←San Pedro River ← Santi Petri Channel •S1 •S2 •S3 •BI4 •BI1 •BI6 •BE4 −18 −16 −14 −12 −10 −8 −6 −4 −2 Fig. 1. Model domain. The color scale gives water depths in m. Locations mentioned in the text are displayed. Arrows indicate main rivers. Each unit in the axis gives the grid cell number. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.) R. Periáñez et al. / Ocean Engineering 69 (2013) 60–69 61 The deposition term is written in the following form: DEP ¼w s m1−jτ !j τ cd  ;ð7Þ where τ !is the bottom friction stress (whose components are given by Eq. (4)) and τ 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' law as mentioned above: w s ¼ρ−ρ w ρ w gD 2 18ν;ð8Þ where ρand Dare suspended particle density and diameter, respectively, and νis the kinematic viscosity of water. The erosion rate is written in terms of the erodability constant E: ER ¼Ef p jτ !j τ ce −1  ;ð9Þ where f p gives the fraction of fine particles in the bed sediment and τ 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. Three phases are then considered: dissolved, suspended matter particles and active bed sediments. These correspond to the fine sediments (particles with diameter o62:5μm) since metals and radionuclides are predominantly fixed to them (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 described by a single reversible reaction, a coefficient k 1 characterizes the transfer from the liquid to the solid phases and a coefficient k 2 characterizes the inverse process. Dimensions of these coefficients are [T] −1 . The adsorption process is a surface phenomenon that depends on the surface of particles per water volume unit. This quantity has been denoted as the exchange surface (Periáñez, 2008,2009). Thus in general k 1 ¼χðS m þS s Þ¼k spm 1 þk sed 1 ;ð10Þ where S m and S s are the exchange surfaces for suspended matter and bottom sediments, respectively (dimensions [L] −1 ). χis a parameter with the dimensions of a velocity. It is denoted as the exchange velocity (Periáñez, 2008,2009). Assuming spherical particles, the exchange surfaces are written as (see references cited above) S m ¼3m ρRð11Þ and S s ¼3Lf p ð1−pÞϕ RH ;ð12Þ where Ris particle radius, pis sediment porosity, Lis the sediment mixing depth (the distance to which the dissolved phase penetrates the sediment) 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,butwiththisapproachitispossibletoobtainananalytical 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þ∂ðvHC d Þ ∂y ¼∂ ∂xHK h ∂C d ∂x  þ∂ ∂yHK h ∂C d ∂y  −Hðk spm 1 þk sed 1 ÞC d þk 2 HmC s þk 2 Lρ s ϕf p A s ;ð13Þ where A s and C s are concentrations in the active fraction of bed sediments and suspended matter, respectively. ρs is the sediment bulk density. The temporal evolution of pollutant concentration in suspended particles is given by ∂ðHC s mÞ ∂tþ∂ðuHC s mÞ ∂xþ∂ðvHC s mÞ ∂y¼∂ ∂xHK h ∂ðC s mÞ ∂x  þ∂ ∂yHK h ∂ðC s mÞ ∂y  þk spm 1 C d H−k 2 mC s HþSED;ð14Þ where SED expresses the pollutant exchange between suspended particles and the bed sediment resulting from erosion/deposition: SED ¼ −SRC s ;SR40 −SRA s ;SRo0: (ð15Þ The equation for the temporal evolution of concentration in the bed sediment is ∂A s ∂t¼k sed 1 C d H Lρ s f p −k 2 A s ϕþSED;ð16Þ where now the exchange due to erosion/deposition of suspended particles is written as SED ¼ SRC s Lρ s f p ;SR40 SRA s Lρ s f p ;SRo0: 8 > > > < > > > : ð17Þ The total concentration of pollutants in the sediment, A tot ,is computed from A tot ¼f p A s :ð18Þ 2.4. Numerical solution All equations are solved using explicit finite difference schemes (Kowalik and Murty, 1993) on a grid with resolution of approximately 200 m. Second order accuracy schemes are used for advective and diffusiveterms.Timesteptosolvethehydrodynamicequations, limited by the CFL condition (Kowalik and Murty, 1993), is Δt¼5s. Since the sediment transport and pollutant dispersion models are run off-line, time step in these calculations was increased to 30 s (stability conditions for hydrodynamic calculations are more restrictive than for the transport ones). The model domain may be seen in Fig. 1, where locations of interest are also shown. Open boundary conditions consist of, in hydrodynamic calculations, specifying sea surface elevations compiled from observations (Alvarez, 1999). In transport calculations, zero gradient is assumed across the open boundary. As commented before, the hydrodynamic model is run off-line. Once it has been calibrated and validated through comparisons of calculated tidal constants with observations, these stored tidal constants are applied to simulate sediment processes. Once realistic R. Periáñez et al. / Ocean Engineering 69 (2013) 60–6962 results on sediment transport are obtained, the pollutant transport model may be simulated. The pollutant transport model is run coupled with the sediment transport one, since instantaneous values of suspended matter concentration and sedimentation rates are required in each point of the domain by the pollutant model. It is worth commenting that the simple addition of tidal and residual currents for different constituents implies that non-linear interactions between flow components are neglected. However, this is a common practice in pollutant transport modelling (see for instance Proctor et al., 1994). 3. Results Results from the hydrodynamic, sediment transport and pollutant dispersion models are presented separately in the following subsections. 3.1. Tide simulations The model has been run for the M 2 ,S 2 ,N 2 ,M 1 ,O 1 and M 4 constituents separately. A calibration procedure has been carried out with the M 2 constituent to select the optimum values of the bed friction coefficient and of horizontal eddy viscosity. Once selected, the same values are used for the other constituents. They are k¼0.007 and A¼10 m 2 /s. As usual, the hydrodynamic model is started from rest and run until stable oscillations are obtained. Then tidal analysis is carried out. As an example, maps showing the tide amplitude and current amplitude for the main constituent, M 2 , are presented in Fig. 2. Residual currents for this constituent are shown in Fig. 3. These maps are in agreement with those presented in Alvarez (1999). It may be seen that the tide amplitude increases by a few cm f rom the bay connection with the Atlantic towards the inner Bay. Maximum currents, as should be expected, are found in the narrow channel connecting the outer and inner bays. Tide residual circulation presents a number of gyres, in agreement with those previously found (Alvarez, 1999). A quantitative comparison between observed and calculated tidal constants for the considered semidiurnal constituents is presented in Table 1, and for the diurnal ones in Table 2. Finally, results for the M 4 constituent are given in Table 3. Generally speaking, there is a good agreement between the observed and calculated tidal constants for all considered constituents. 3.2. Sediment dynamics The sediment transport model is started from a uniform suspended matter concentration over all the domain and, as in the case of tides, run until stable oscillations in suspended matter are obtained. One of the difficulties of modelling sediment transport is defining parameters appearing in the model. In the present application, the deposition and erosion threshold stresses are fixed, respectively, as 0.2 and 0.35 N/m 2 , and the erodability constant as 1.1 10 −4 kg/m 2 s. They are within typical ranges found in literature (Tattersall et al., 2003). The characteristic diameter of suspended particles has been defined as 4 μm. This value has been selected after model calibration, and is of the same order as the value used in a similar modelling study for the Alborán Sea (Periáñez, 2008), which was 8 μm. A standard value of 2600 kg/m 3 has been used for particle density ρ. The effective diffusion coefficient has been set as 15 m 2 /s. With this set of parameters, results in agreement with available data are obtained. Three suspended sediment sources are included: San Pedro River, Guadalete River and Santi Petri Channel (see Fig. 1). As stream flows are extremely low, input of particles is produced only during the ebb tide. It is known that sediments in the inner bay are predominantly muds, since this area is protected from winds and waves (Ligero et al., 2005). On the other hand, sediments are predominantly sandy in areas with strong tidal currents, as the entrance of the bay and the channel connecting the inner and outer bays (Araujo et al., 2009;Carrasco et al., 2003). Thus, in a schematic way, parameter f p (which gives the fraction of muddy 10 20 30 40 0 10 20 30 40 50 60 70 80 90 100 110 X Y 1 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 10 20 30 40 0 10 20 30 40 50 60 70 80 90 100 110 X Y 0 0.1 0.2 0.3 0.4 0.5 Fig. 2. Calculated M 2 amplitude (m, left) and current amplitude (m/s, right). Blues in the inner bay (left panel) correspond to intertidal areas. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this article.) R. Periáñez et al. / Ocean Engineering 69 (2013) 60–69 63 sediments) has been fixed as 0.2 in the channel, 0.45 in the outer bay and 1.0 in the inner bay. As will be shown, this schematic representation of the seabed composition is enough to obtain realistic results. As an example, suspended particle concentrations at four instants during a tidal cycle (high water, ebb, low water and flood) are presented in Fig. 4. Minimum concentrations are observed during high water, since there are not riverine inputs and already present particles are diluted in a larger water volume. During the ebb and at low water, inputs from river are apparent. Finally, during the flood period the outer bay is relatively clear due to the input of water from the sea, while higher suspended matter concentrations are observed in the inner bay. Indeed, the flood tide will erode sediments of the channel connecting the inner and outer bays, and these resuspended sediments will be transported towards the inner bay, where they are concentrated. Suspended matter concentrations have been compared with measurements in 4 points of the domain (see Fig. 1). Sampling started 1.5 h before high water and ended 1.5 h after high water, but the exact sampling time is not known for each sample. Thus, we have compared the range of suspended matter concentrations calculated in such time frame with measurements. Results are presented in Table 4. It may be seen that computed values are reasonable, similar to measured concentrations. Sedimentation rates have also been compared with those obtained from measurements in sediment cores (Ligero et al., 2002). Location of sampling points is indicated in the map in Fig. 1 and a comparison between measurements and model results is presented in Table 5. Calculated results in Table 5 have been obtained averaging the instantaneous sedimentation rates provided by the model (each 30 s) over a period of 60 days. The model produces sedimentation rates of the order of 0.20 cm/y, similar to measured values. Table 1 Observed, subscript obs, and computed, subscript comp, amplitudes (A, m) and phases (g, deg) of semidiurnal tidal elevations at several locations indicated in Fig. 1. Station M 2 S 2 N 2 A obs g obs A comp g comp A obs g obs A comp g comp A obs g obs A comp g comp Carraca 1.08 60 1.08 60 0.35 86 0.38 85 0.27 41 0.28 35 Pto Real 1.07 58 1.08 60 0.39 82 0.37 85 0.29 43 0.28 35 P. Carranza 1.07 58 1.07 58 0.36 87 0.37 84 0.27 38 0.27 35 Pto Cádiz 1.03 55 1.04 54 0.36 84 0.36 82 0.26 37 0.26 33 Pto Sherry 1.03 53 1.04 54 0.35 81 0.36 82 0.26 33 0.26 33 Rota 1.01 56 1.04 54 0.33 81 0.35 82 0.25 48 0.26 34 Table 2 Observed, subscript obs, and computed, subscript comp, amplitudes (A, m) and phases (g, deg) of diurnal tidal elevations at several locations indicated in Fig. 1. Station K 1 O 1 A obs g obs A comp g comp A obs g obs A comp g comp Carraca ––0.062 43 0.055 302 0.051 301 Pto Real ––0.062 43 0.055 302 0.051 301 P. Carranza 0.063 46 0.062 43 0.053 303 0.051 301 Pto Cádiz 0.067 41 0.061 42 0.052 302 0.050 300 Pto Sherry 0.064 42 0.061 42 0.054 300 0.050 300 Rota 0.064 38 0.062 42 0.058 302 0.051 301 Table 3 Observed, subscript obs, and computed, subscript comp, amplitudes (A, m) and phases (g, deg) of M 4 tidal elevations at several locations indicated in Fig. 1. Station A obs g obs A comp g comp Carraca 0.028 302 0.024 324 Pto Real 0.026 309 0.024 324 P. Carranza 0.028 309 0.023 323 Pto Cádiz 0.021 315 0.020 319 Pto Sherry 0.020 308 0.020 320 Rota 0.018 312 0.019 320 10 20 30 40 10 20 30 40 50 60 70 80 90 100 110 5 cm/s X Y Fig. 3. Calculated M 2 residual currents. R. Periáñez et al. / Ocean Engineering 69 (2013) 60–6964 Once it is known that the model produces realistic suspended matter concentrations and sedimentation rates, it may be applied to carry out particle reactive pollutant dispersion experiments. 3.3. Pollutant dispersion 3.3.1. Dispersion experiments Some additional parameters have to be defined for the pollutant transport model. Following previous work (Periáñez, 2000, 020 40 0 20 40 60 80 100 X Y High 0 10 20 30 40 020 40 0 20 40 60 80 100 X Y Ebb 0 10 20 30 40 020 40 0 20 40 60 80 100 X Y Low 0 10 20 30 40 020 40 0 20 40 60 80 100 X Y Flood 0 10 20 30 40 Fig. 4. Suspended matter concentration (mg/L) maps at 4 tidal states. Table 4 Measured spm concentrations 1.5 h around high water and calculated range of values in the same time frame. Station Measured spm (mg/L) Calculated range (mg/L) BI4 7.4 1.0–5.6 BI1 8.0 1.7–12.3 BI6 12.2 6.0–11.8 BE4 16.0 9.0–13.4 Table 5 Measured and computed sedimentation rates in three points of the domain. Station Measured (cm/y) Calculated (cm/y) S1 0.16 0.12 S2 0.22 0.20 S3 0.27 0.24 R. Periáñez et al. / Ocean Engineering 69 (2013) 60–69 65 2003,2009), L¼0.1 m and ϕ¼0:1. Also, a constant porosity has been assumed for sediments, fixed as 0.6 from Ligero et al. (2005). Kinetic rates must be defined for each simulated pollutant. As widely described before (Periáñez, 2005a,2008,2009), the same value may be used for k 2 for different elements since the main parameter controlling pollutant geochemical behavior is χ. Following the procedure described in such references, the exchange velocity may be deduced from k 2 and the pollutant equilibrium distribution coefficient k d since the following relation holds (see the same references): k d ¼χ k 2 3 ρR:ð19Þ The value k 2 ¼1.16 10 −5 s −1 has been fixed. It has been successfully used in earlier simulations in the Odiel-Tinto estuary (Periáñez et al., 2005) for Ra, in the Strait of Gibraltar-Alborán Sea (Periáñez, 2008) for Cs and Pu, and in the Gulf of Cádiz (Periáñez, 2009)forheavy metals. It was measured for Cs by Nyffeler et al. (1984). Afirst numerical experiment has been carried out to simulate the adsorption of fallout 137 Cs (no other anthropogenic sources exist) on seabed sediments. A k d of 2.0 m 3 /kg has been used since it agrees with the established value by IAEA (2004) for coastal waters and also was successfully applied in Periáñez (2008). Since fallout has been occurring since several decades and with varying intensity, we have just made a qualitative comparison of the calculated 137 Cs inventories (total amounts) in the sediment mixed layer with the measured inventories in Ligero et al. (2005). Thus, we have made a 3 month simulation with a constant and arbitrary atmospheric flux of 137 Cs. The calculated inventories in each grid cell have been normalized in such a way that the maximum calculated inventory is the same as the maximum measured one. Then, both spatial distributions are compared. Results are presented in Fig. 5. Although measurements are available only in the inner bay, measured and calculated distributions are similar, indicating a significant accumulation of caesium in the zone of Carraca station (see Fig. 1). Model results are normalized, thus the significant aspect of this test is that the model predicts the accumulation of a pollutant introduced via atmospheric deposition in a given area of the bay, and this is in agreement with observations. Nevertheless, a second, but weaker, accumulation area in the western inner bay is not reproduced by the model. Calculated levels of 137 Cs in the outer bay are lower, due to the higher water exchange with the open sea and subsequent sediment cleaning. As commented before, a number of pollutant sources exist in Cádiz Bay (industrial and urban wastewater discharges), whose magnitudes are not known. Thus, it makes no sense trying to compare calculated metal distributions with measurements. This would imply setting the magnitude of each source, by trial and error, until results in agreement with observations are obtained. Consequently, we have just carried out some simple numerical experiments to study the pollutant distributions resulting from some sources. An experiment has been carried out with Zn, which is a much more reactive element than Cs. Its distribution coefficient is (IAEA, 2004)70m 3 /kg, which has been used in previous simulations for the Gulf of Cádiz (Periáñez, 2009). Carrasco et al. (2003) have found high Zn contents in waters of both the inner and outer bays, without significant differences between tidal conditions. These Zn concentrations were particularly high in the channel connecting the inner and outer bays, in the area of Guadalete River (see Fig. 6 for its location) and in Pto. Cádiz site (Fig. 1). A numerical experiment has been carried out introducing a constant Zn release, of arbitrarily defined magnitude, from the Guadalete River. Calculated Zn concentrations in water, suspended matter and bed sediments after 60 days of release are shown in Fig. 6. It seems clear that contaminants are pumped into the inner bay, where they remain trapped, confirming the existence of an accumulation area in the zone of Carraca station, as already found with the Cs experiment. This accumulation area must be related to residual circulation in the inner bay, where several gyres (see Fig. 3) retain substances making them to recirculate. Thus any pollutant release, even if it occurs in the outer bay, may have a significant environmental impact on the inner bay and has to be planned with extreme care. In this sense, it is relevant to investigate flushing times of pollutants for the inner bay, which are analyzed in the following section. Recently, an extensive sampling campaign has been carried out in Cádiz Bay (Casas-Ruíz et al., 2012), and 226 Ra concentrations have been measured in sediments. A map of the measured concentrations is presented in Fig. 7. A simulation has been carried out introducing a source of arbitrary magnitude from Guadalete River, where the maximum concentrations have been measured. The recommended k d value for Ra by IAEA (2004) is 2 m 3 /kg, 10 20 30 40 10 20 30 40 50 60 70 80 90 100 110 X Y −1 −0.5 0 0.5 1 1.5 2 2.5 3 Fig. 5. Calculated (top, normalized logarithms) and measured (bottom, Ligero et al., 2005) 137 Cs inventories in bed sediments (Bq/m 2 ). Note that the model grid is rotated 30.81clockwise. The black rectangle on the top panel indicates the region where measurements indicated in the lower panel were carried out. R. Periáñez et al. / Ocean Engineering 69 (2013) 60–6966 which has been used in our simulation. The computed 226 Ra distribution in sediments is also presented in Fig. 7. Measured and computed distributions are somewhat similar, with lower concentrations in the channel connecting the inner and outer bays and higher values again in the inner bay. As in the case of the 137 Cs experiment, the calculated concentrations have been normalized. Even so, calculated concentrations cannot be compared with the measured ones (note for instance that calculated concentrations in the inner bay are lower than measured ones). This is due to the fact that 226 Ra is a natural radionuclide, also having a geological (thus diffuse) origin (Casas-Ruíz et al., 2012). We are not including this diffuse source, but only a point source in Guadalete River. It is not surprising that calculated concentrations are lower than measured ones. Nevertheless, the overall calculated 226 Ra distribution seems to be reasonable. 3.3.2. Flushing time of the inner bay Flushing time is defined as the time in which the tracer inventory in the water column (total amount of tracer in the considered region) decreases in a factor e(Prandle, 1984). Flushing time is relevant for the water quality of a system, and it is important to know the time scale for a pollutant discharged into a water body to be transported out of the system (Shen and Haas, 2004). Thus, flushing times have been recently determined for a number of water bodies using numerical experiments (Shen and Haas, 2004;Choi and Lee, 2004;Patgaonkar et al., 2012;Periáñez, 2012). Flushing time has been determined considering a perfectly conservative (dissolved) contaminant. Model initial conditions consist of setting a 100 units/m 3 concentration at every grid cell inside the inner bay and zero concentration elsewhere. The model is integrated and the time evolution of the pollutant inventory inside the bay written to an output file. Fitting to exponential decay curves provides the flushing time. Thus, the system-wide (as defined in Choi and Lee, 2004)flushing time has actually been determined. These authors found that this parameter is useful for determining the long-term water quality of a system. Choi and Lee (2004) have found that the system-wide flushing time is better obtained from a double exponential decay curve rather than a single exponential decay. The temporal evolution of the mass within the system, M(t), is thus given by the following equation: MðtÞ M 0 ¼ð1þβÞe −α 1 t −βe −α 2 t ;ð20Þ where M 0 is the initial mass within the system. If the three parameters β,α 1 and α 2 are determined from numerical fitting, the system-wide flushing time is given by (Choi and Lee, 2004) T f ¼1þβ α 1 −β α 2 :ð21Þ As an example, the temporal evolution of the normalized tracer inventory, M/M 0 in the inner bay may be seen in Fig. 8. Numerical fitting to the double exponential provides a flushing time of 63.1 days, with a determination coefficient for the fitting r 2 ¼0.9966. This time is considerably longer than in Algeciras Bay, for instance (Periáñez, 2012). This is not surprising since, although tides are stronger in Cádiz Bay than in Algeciras Bay, the inner Cádiz Bay is a much closed water body. Also, pollutants are trapped by residual circulation as commented above. Indeed, an additional numerical experiment has been carried out in which advective terms in the dispersion equation have been removed, i.e., like there is not any water circulation at all and the only transport process is diffusion. The computed flushing time is of the same order as that obtained with all tidal constituents, confirming the extremely low capacity of tides to flush pollutants off the inner bay. 4. Conclusions A numerical model to simulate the dispersion of particlereactive contaminants in Cádiz Bay (SW Spain) has been developed. 020 40 0 10 20 30 40 50 60 70 80 90 100 110 X Y Dissolved −4 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 020 40 0 10 20 30 40 50 60 70 80 90 100 110 X SPM −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 020 40 0 10 20 30 40 50 60 70 80 90 100 110 X Sediment −2.5 −2 −1.5 −1 −0.5 0 0.5 1 Fig. 6. Calculated Zn concentrations in water (units/m 3 ), suspended matter and bed sediments (both in units/kg). The black dot indicates the Gualdalete River mouth, where a constant release of arbitrary magnitude occurs. Logarithms of concentrations are drawn. R. Periáñez et al. / Ocean Engineering 69 (2013) 60–69 67 It consists of a hydrodynamic model, a sediment transport model and a pollutant transport model. Hydrodynamics is calculated in advance, and the sediment and pollutant models are run off-line. Calculated water circulation and suspended sediment dynamics is generally reproduced by the model. The contamination of seabed sediments with fallout 137 Cs has been studied through a numerical experiment. It has been found that there is a accumulation area in the inner bay. This was already found in field measurements. Thus, the model is correctly describing transport processes inside the bay. If a continuous source of pollutants is present, even if it is in the outer bay, there is a pumping of the contaminant to the inner bay, where they are accumulated in the sediments of the area already found with the 137 Cs experiment. The existence of such accumulation area must be attributed to water circulation in the inner bay. Finally, flushing time of the inner bay has been determined. A value of the order of 60 days has been obtained. This time is significantly longer than in other environments with a weaker tidal forcing (for instance Algeciras Bay), since this is a much closed environment. Also, residual circulation produces an accumulation area of pollutants in the inner bay, where they remain trapped. This increases the flushing time. Acknowledgment Work supported by the project “Determination of scavenging rates and sedimentation velocities using reactive-particle radionuclides in coastal waters; application to pollutant dispersion modeling”, Refs. CTM2009-14321-C02-01 and CTM2009-14321-C02-02, funded by the Spanish Ministerio de Ciencia e Innovación. References Alvarez, O., 1999. Simulación numérica de la dinámica de marea en la bahía de Cádiz: análisis de las constituyentes principales, interacción marea-brisa e influencia del sedimento en suspensión. Ph.D. Thesis, University of Cádiz (in Spanish). Alvarez, O., Izquierdo, A., Tejedor, B., Mañanes, R., Tejedor, L., Kagan, B.A., 1999. The influence of sediment load on tidal dynamics, a case study: Cádiz Bay. Estuarine Coastal Shelf Sci. 48, 439–450. Alvarez, O., Tejedor, B., Tejedor, L., Kagan, B.A., 2003. A note on sea-breeze induced seasonal variability in the K 1 tidal constants in Cádiz Bay, Spain. Estuarine Coastal Shelf Sci. 58, 805–812. Araujo, C.V.M., Diz, F.R., Laiz, I., Lubián, L.M., Blasco, J., Moreno-Garrido, I., 2009. Sediment integrative assessment of the Bay of Cádiz (Spain): an ecotoxicological and chemical approach. Environ. Int. 35, 831–841. Carrasco, M., López-Ramírez, J.A., Benavente, J., López-Aguayo, F., Sales, D., 2003. Assessment of urban and industrial contamination levels in the Bay of Cádiz, SW Spain. Mar. Pollut. Bull. 46, 335–345. Casas-Ruíz, M., Ligero, R.A., Barbero, L., 2012. Estimation of annual effective dose due to natural and man-made radionuclides in the metropolitan area of the Bay of Cádiz (SW of Spain). Radiat. Prot. Dosimetry 150, 60–70. Choi, K.W., Lee, J.H.W., 2004. Numerical determination of flushing times for stratified water bodies. J. Mar. Syst. 50, 263–281. Delhez, E.J.M., 1996. On the residual advection of passive constituents. J. Mar. Syst. 8, 147–169. Duursma, E.K., Carroll, J., 1996. Environmental Compartments. Springer, Berlin. Dyke, P.P.G., 2001. Coastal and Shelf Sea Modelling. Kluwer, The Netherlands. Eisma, D., 1993. Suspended Matter in the Aquatic Environment. Springer-Verlag, Berlin. Gómez-Parra, A., Establier, R., Escolar, D., 1984. Heavy metals in recent sediments from the Bay of Cádiz, Spain. Mar. Pollut. Bull. 15, 307–310. IAEA, 2004. Sediment Distribution Coefficients and Concentration Factors for Biota in the Marine Environment. Technical Reports Series 422, Vienna. Kagan, B.A., Tejedor, L., Alvarez, O., Izquierdo, A., Mañanes, R., 2001. Weak wavetide interaction formulation and its application to Cádiz Bay. Cont. Shelf Res. 21, 697–725. Kampf, J., 2009. Ocean Modelling for Beginners. Springer-Verlag, Heidelberg. Kowalik, Z., Murty, T.S., 1993. Numerical Modelling of Ocean Dynamics. World Scientific, Singapore. Ligero, R.A., Barrera, M., Casas-Ruíz, M., Sales, D., López-Aguayo, F., 2002. Dating of marine sediments and time evolution of heavy metal concentrations in the Bay of Cádiz, Spain. Environ. Pollut. 118, 97–108. Ligero, R.A., Barrera, M., Casas-Ruíz, M., 2005. Levels of 137 Cs in muddy sediments of the seabed of the Bay of Cádiz, Spain. Part I. Vertical and spatial distribution of activities. J. Environ. Radioact. 80, 75–86. Liu, W.C., Hsu, M.H., Kuo, A.Y., 2002. Modelling of hydrodynamics and cohesive sediment transport in Tanshui River estuarine system, Taiwan. Mar. Pollut. Bull. 44, 1076–1088. Lumborg, U., Windelin, A., 2003. Hydrography and cohesive sediment modelling: application to the Romo Dyb tidal area. J. Mar. Syst. 38, 287–303. 010 20 30 40 10 20 30 40 50 60 70 80 90 100 110 X Y −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 -6.3 -6.25 -6.2 -6.15 36.45 36.5 36.55 36.6 2 4 6 8 10 12 14 16 18 20 22 24 Ra 226 Bq/kg Longitude ( W) Latitude ( N) Fig. 7. Calculated (top) 226 Ra concentrations in sediments (Bq/kg), and measured (bottom) concentrations (Bq/kg). Logarithms of calculated concentrations are drawn. 010 20 30 40 50 60 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 Time (days) M/M0 r2=0.9966 Fig. 8. Time evolution of the normalized inventory of a dissolved contaminant in the inner bay. Fitting to the double exponential in Eq. (20) is also shown (solid line). R. Periáñez et al. / Ocean Engineering 69 (2013) 60–6968