Measuring and modelling temporal trends of 226Ra in waters of a Spanish estuary affected by the phosphate industry

Measuring and modelling temporal trends of 226 Ra in waters of a Spanish estuary affected by the phosphate industry R. Peria ´n˜ez * Dpto. Fı ´sica Aplicada 1. EUITA, Universidad de Sevilla, Ctra Utrera km 1, 41013-Sevilla, Spain Received in revised form 27 July 2004; accepted 16 August 2004 Abstract The presence and temporal evolution (1990–2001) of 226 Ra in a tidal estuary affected by the phosphate industry has been investigated. Water samples collected in the course of four separate sampling campaigns were analysed for 226 Ra content using a gas flow proportional counter following Ba coprecipitation. Two 226 Ra sources have been identified: direct discharges from the industrial complex and run-off from a phosphogypsum pile. Although activity levels are similar, or even higher, than those found in other environments affected by the phosphate industry, there has been a general decrease in contamination since direct discharges ceased in 1998 due to new regulations from the EU. However, sediments are now acting as a source of Ra to the water column due to redissolution processes. A numerical model of the estuary has been developed to describe quantitatively the experimental results. The model solves the hydrodynamics and the dispersion equation of 226 Ra including interactions with sediments. Model results are, in general, in good agreement with observations. 2004 Elsevier Ltd. All rights reserved. Keywords: Radium; Odiel-Tinto estuary; Phosphogypsum; Numerical modelling; Hydrodynamic; Sediment 0141-1136/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.marenvres.2004.08.003 * Tel.: +34 954486474; fax: +34 954486436. E-mail address: [email protected]. Marine Environmental Research 60 (2005) 35–49 www.elsevier.com/locate/marenvrev MARINE ENVIRONMENTAL RESEARCH 1. Introduction The Odiel and Tinto rivers, in the southwest of Spain, form a fully mixed tidal estuary which surrounds the town of Huelva (Fig. 1). Both rivers join at the Punta del Sebo. From this point, they flow together to the Atlantic Ocean. An industrial complex, containing a plant dedicated to the production of phosphoric acid and phosphate fertilizers, is located by the Odiel river. It is well known that the phosphate rock used as a raw material by this industry contains significant amounts of natural radionuclides, mostly U, Th and Ra. The industrial processing of the phosphate rock leads to a redistribution of radioactivity. For instance, during the wet process for phosphoric acid production, 86% of U and 70% of Th present in the rock are transferred to the phosphoric acid itself, while 80% of the Ra content follows the so-called phosphogypsum path (Guimond & Hardin, 1989). This is a form 10 20 30 40 50 60 10 20 30 40 50 60 70 80 01234 Canal del Burro km Tinto River Odiel River N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 fertilizer complex 16 P.Sebo HUELVA phosphogypsumpiles E D C B A Fig. 1. Map of the area of the estuary covered by the model. Numbered circles indicate the points where water samples were collected. Lettered triangles indicate the points where currents measurements were available. Units on the axes give the grid cell number. The sea is approximately 1 km to the south of point 1. 36 R. Peria ´n˜ez / Marine Environmental Research 60 (2005) 35–49 of impure calcium sulfate removed as a precipitate during the process. Phosphogypsum is usually disposed into piles in the open environment or discharged into rivers or estuaries, giving rise to a local radioactive impact. The raw material used in the plant is phosphate rock from Senegal, Morocco and Togo, which is treated by the wet process to obtain phosphoric acid. During 1990, for instance, 2 ·10 6 tons were processed and 3 ·10 6 tons of phosphogypsum were produced. These wastes were partially released directly into the Odiel river (20%), with the remaining 80% conducted with water through a pipeline to phosphogypsum piles located by the Tinto river (see Fig. 1), where such material is stored in the open air. The gypsum piles cover some 12 km 2 of the Tinto river margin. Since 1998, wastes are not released directly into the Odiel river due to new regulations from the EU, although phosphogypsum is still being disposed of in the piles by the Tinto river. A number of recent studies have investigated the radioactive impact of the phosphate industry in its close environment. Paridaens and Vanmarcke (2001) charted the 226 Ra contamination of the river Laak (Belgium) and some areas that are regularly flooded by the river. Haridasan, Paul, and Desai (2001) found that the surface runoff of water from a phosphogypsum pile placed by the Chitrapuzha river (India) appears to be the major pathway of transport of Ra to the river. McCartney, Davidson, Howe, and Keating (2000) studied the effect of a reduction of the discharges from a phosphate plant in the levels of U and Th isotopes along the Cumbrian coast (UK), while Poole, Allington, and Denoon (1995) conducted a similar study but for 226 Ra. The effect of the phosphate industry in the Odiel and Tinto estuary has also been investigated. Martı ´nez-Aguirre, Garcı ´a-Leo ´n, and Ivanovich (1994a) and Bolı ´var, Garcı ´a-Tenorio, and Vaca (2000) studied the presence of U and Th isotopes over different areas of the estuary. The presence of 226 Ra in waters of the estuary was investigated by Peria ´n˜ez and Garcı ´a-Leo ´n (1993), who determined 226 Ra activity concentrations in waters collected from the estuary in two different sampling campaigns in July 1990 and March 1991. The objective of the present work is to provide a summary of the work carried out at the University of Seville concerning the temporal evolution of the 226 Ra contamination of the estuary, presenting results for sampling campaigns carried out before and after the change in the waste policy of the phosphate plant in 1998. To this effect, two new sampling campaigns were conducted; in October 1999 and February 2001. Results obtained for the four campaigns are presented together. Although the results for the 1990 and 1991 samples have been published previously (Peria ´n˜ez & Garcı ´aLeo ´n, 1993), they will also be included here for ease of comparison. It has to be pointed out that sampling campaigns are non-ideal for a temporal trend study, since have been carried out in different months during the four years. Nevertheless, some interesting conclusions may be obtained from them, as will be shown. Interpretation of experimental results is not an easy task since a number of processes affect concentrations in the dissolved phase. Currents due to tidal oscillations produce advective transport of dissolved radionuclides along the estuary. Furthermore, turbulent diffusion is present and there are also transfers of radionuclides between the dissolved and solid phases (suspended matter particles and bottom R. Peria ´n˜ez / Marine Environmental Research 60 (2005) 35–49 37 sediments). The solid phase may act as a sink or as a source of radionuclides to the water column, depending on the concentration of radionuclides in each phase and on the rates governing the adsorption and release reactions. To account for these processes, a numerical model of the estuary, including all these factors, has been developed. The model solves the hydrodynamics of the estuary and the dispersion of 226 Ra, including the interactions between the dissolved and solid (bottom sediments) phases. The model application allows reproduction of the experimental results for the different campaigns and an estimate to be made of the input of 226 Ra from the phosphogypsum piles to the Tinto river. Details of the sampling and the experimental method for 226 Ra extraction and measurement are presented in the next section. This is followed by the presentation and discussion of results. Finally, the model is described and applied. 2. Sampling and analysis Sampling campaigns were carried out on 19th July 1990, 5th March 1991, 27th October 1999 and 14th February 2001, under low water conditions (1 h after low water in the case of 1991). Samples (see Fig. 1 for positions of the sampling points) were collected in 25 l plastic bottles. The water was filtered soon after sampling, using 0.45 lm pore size Nuclepore filters to remove suspended particles. Filtered samples were acidified to pH 2–3 with HNO 3 to avoid the growth of microorganisms and minimize water–wall interactions during storage. 226 Ra activities were determined from 0.5 l aliquots of the water samples. After neutralisation of the filtrate using NH 4 OH, 5 mg of BaCl 2 and 20 ml of 1 M H 2 SO 4 were added to the sample. Under these conditions Ba–Ra sulfate precipitates after some 20 min of continuous stirring and moderate heating. The precipitate, containing the radium, was collected by filtration through a 0.45 lm Millipore filter. Activity on the filter was measured using a low background gas flow proportional counter previously calibrated for total efficiency vs. precipitate mass thickness. These procedures have been widely validated and applied. Further details on the method can be found in Moro ´n, Garcı ´a-Tenorio, Garcı ´a-Montan˜o, Garcı ´a-Leo ´n, and Madurga (1986),Martı ´nez-Aguirre and Garcı ´a-Leo ´n (1994) and Peria ´n˜ez and Garcı ´a-Leo ´n (1993). 3. Experimental results 226 Ra concentrations for the four campaigns are presented in Table 1. The levels measured in the first three campaigns are broadly similar (or somewhat higher) to those found in other rivers affected by fertiliser industries. In the Boben River (Slovenia), for example, 226 Ra concentrations of 40 mBq/l have been reported in waters from the most affected area (Kobal, Brajnik, Kaluza, & Vengust, 1990), while values in the range 9–68 mBq/l have been reported for affected waters in the Schelde, Netherlands (Koster et al., 1991). More recently, 226 Ra concentrations of 30 mBq/l have 38 R. Peria ´n˜ez / Marine Environmental Research 60 (2005) 35–49 been measured downstream a phosphogypsum pile in the Chitrapuzha River, India (Haridasan et al., 2001). In the case of the 1990 samples, there is an intense peak of 670 mBq/l at point 8, located close to the pipeline through which direct discharges to the Odiel river were carried out. This shows a clear local radioactive impact in the estuary due to the operation of the fertilizer plant. For the remainder of the estuary, activity concentrations range from 64 to 106 mBq/l, also suggesting a clear contamination if they are compared with values reported above. For samples collected in 1991, the peak in the vicinity of the plant outlet is not apparent. However, activities ranging from 21 to 66 mBq/l are indicative of contamination still being present, although the input from the plant to the Odiel river seems less than during the 1990 sampling. It is interesting to note that there is an activity peak in sample 15, collected in the Tinto river close to the phosphogypsum pile. It seems that some run-off from the phosphogypsum pile to the river is occurring and, as a consequence, 226 Ra is being transported to the river, as has been found in other environments (Haridasan et al., 2001). This is related to the fact that the phosphogypsum storage area is crossed by a number of small streams of natural and artificial origin (Peria ´n˜ez & Garcı ´a-Leo ´n, 1993). Some of these may be possible remnants of the wet lands totally integrated in the estuarine system of the rivers only 40 years ago, before industrial activity started in the area (Borrego & Pendo ´n, 1988). Presumably, these streams are capable of transporting activity from the storage area to the main river. This effect was also probably taking place in 1990, although it may have been masked by the high contamination produced by the large discharges that were released into the Odiel. This point will be investigated later with the numerical model of the estuary. Table 1 226 Ra activities (mBq/l) in water samples of the Odiel and Tinto rivers Sample 1990 1991 1999 2001 1 3.6 ± 0.6 9.3 ± 1.2 2 19.1 ± 1.5 8.9 ± 1.4 3 16.4 ± 1.8 9.2 ± 1.2 4 71 ± 2 54.2 ± 1.8 19.6 ± 2.1 5 82 ± 1 46.6 ± 1.6 16.4 ± 1.6 12.5 ± 1.9 6 13.0 ± 0.6 9.4 ± 1.7 7 106 ± 3 46.5 ± 1.6 17.9 ± 2.2 8.9 ± 2.1 8 670 ± 13 53 ± 2 11.0 ± 0.6 9.4 ± 1.1 9 106 ± 2 25.3 ± 1.2 11.3 ± 1.0 6.6 ± 0.6 10 86 ± 2 21.3 ± 1.4 6.8 ± 0.6 7.2 ± 0.5 11 15.9 ± 1.7 9.7 ± 0.9 12 64 ± 2 33.6 ± 1.5 45 ± 5 10 ± 4 13 69.3 ± 1.9 36.0 ± 1.4 35.3 ± 1.8 8.4 ± 0.6 14 29 ± 4 9.7 ± 1.0 15 67 ± 1 65.6 ± 1.9 25.1 ± 2.3 8.6 ± 0.9 16 14.0 ± 0.9 8.6 ± 1.0 Uncertainties (1r) are due to counting statistics in sample measurements and in detector calibration. Empty spaces mean that such samples were not collected. R. Peria ´n˜ez / Marine Environmental Research 60 (2005) 35–49 39 Activity levels detected in the 1999 campaign, ranging from 3.6 to 45 mBq/l, are lower than those obtained in 1990 and 1991. This is likely to be related to the introduction of stricter waste policies in 1998, which resulted in the cessation of direct discharges into the Odiel River. As a consequence, the main source of 226 Ra to the estuary that year must have been run-off from the phosphogypsum pile. That this is the case is evidenced by the higher concentrations observed along the Tinto River in comparison to the Odiel River, and by the relative uniformity of measured concentrations in the latter. The pattern is very different from that observed in previous years, when direct discharges were taking place. Nevertheless, concentrations along the Odiel River are still higher than those found in non-perturbed rivers, which are typically in the range 0.09–3.4 mBq/l (Bhat & Krishnaswamy, 1969; Rona & Urry, 1952). The enhanced concentrations can be attributed to the redissolution of 226 Ra from a pool of underlying contaminated sediments. Indeed, there is now strong evidence to suggest that contaminated sediments can become a source of remobilised radionuclides when the external input is reduced and desorption reactions dominate over adsorption (Cook, MacKenzie, McDonald, & Jones, 1997; McCartney et al., 2000). This hypothesis will be studied in more detail with the help of the numerical model described below. Activity levels detected in samples collected in 2001 are rather uniform over all the estuary, ranging from 6.6 to 12.5 mBq/l. Unlike the 1999 data, there is no evidence of run-off of 226 Ra from the phosphogypsum pile to the Tinto River. This is not unexpected, given that earth dikes were built around the piles to prevent run-off in the period between the two sampling campaigns. It may, therefore, be assumed that the only important source of 226 Ra to the waters of the estuary at the time of this last campaign is remobilisation of this radionuclide from previously contaminated sediments. 4. Model description The system under study is divided into a number of grid cells or compartments. Two phases are present in each grid cell: dissolved and active bottom sediments (particles with a diameter <62.5 lm). The active sediments correspond to muddy sediments, following the Wentworth scale of sediment grain size (see for instance Pugh, 1987). Larger grain sizes are not considered since it has been shown that virtually all the radioactivity is associated with the muddy sediment (Aston, Assinder, & Kelly, 1985). Suspended matter particles have not been considered in the model, and thus deposition processes and erosion of the sediment have been neglected. This approximation is used since previous calculations have shown that the radionuclide adsorption capacity of suspended matter, given the typical suspended matter concentrations in the estuary (maximum concentrations of the order of 50 ppm) is very small compared with that of the sediment (Peria ´n˜ez, Abril, & Garcı ´a-Leo ´n, 1996a). Moreover, the erosion-deposition rates, obtained from a suspended matter model of the estuary (Peria ´n˜ez et al., 1996a), are also small (10 2 g/cm 2 /year). Thus, as an approximation, it has been considered that the most important phases controlling radionuclide 40 R. Peria ´n˜ez / Marine Environmental Research 60 (2005) 35–49 transport are the dissolved phase and the bottom sediment. This approximation seems realistic given the generally good agreement between model results and observations (see below). Adsorption and desorption reactions are described in terms of kinetic transfer coefficients. Thus, the adsorption process (transfer of radionuclides from water to the sediment) will be governed by a coefficient k 1 and the inverse process (desorption to the dissolved phase) by a coefficient k 2 . The adsorption process is a surface phenomenon that depends on the surface of particles per water volume unit into the grid cell. Following the notation of Peria ´n˜ez, Abril, and Garcı ´a-Leo ´n (1996b), the adsorption coefficient is written as: k1¼v1SE;ð1Þ where S E is the exchange surface and v 1 is a parameter with dimensions of velocity, denoted as the exchange velocity. As a first approach, assuming spherical sediment particles and a step function for the grain size distribution of particles, it can be shown (Peria ´n˜ez et al., 1996b) that SE¼3Lf /  rH ;ð2Þ where  ris the mean radius of sediment particles, His the total water depth, Lis the average mixing depth (the distance to which the dissolved phase penetrates the sediment), fgives the fraction of active sediment and /is a correction factor that takes into account that not all the exchange surface of the sediment particle is in contact with water since part of it can be hidden by other particles (thus, it is also related to sediment porosity). Consequently, it is implicitly assumed that radionuclide concentrations in sediment pore waters, considering a sediment layer of thickness Linside which the sediment is homogeneous, are equal to those in the water column. On the other hand, diffusion of radionuclides to deeper sediment layers has been neglected given the time scale of the simulations carried out. The transfer coefficient k 2 is considered constant. Of course, dissolved radionuclides will be transported along the estuary by advective and diffusive processes. Therefore, the hydrodynamic equations must be solved too. 4.1. Hydrodynamic equations and physical characteristics of the estuary The 2D shallow water hydrodynamic equations are (see for instance Pugh, 1987): oz otþo ox½ðDþzÞuþo oy½ðDþzÞv¼0;ð3Þ ou otþuou oxþvou oyþgoz oxXvþKuffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þv2 p Dþz¼0;ð4Þ ov otþuov oxþvov oyþgoz oyþXuþKvffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2þv2 p Dþz¼0;ð5Þ R. Peria ´n˜ez / Marine Environmental Research 60 (2005) 35–49 41 where uand vare the depth averaged water velocities along the xand yaxis, Dis the depth of water below the mean sea level, zis the displacement of the water surface above the mean sea level measured upwards, Xis the Coriolis parameter (X=2wsinb, where wis the earth rotational angular velocity and bis latitude), g is acceleration due to gravity and Kis the bed friction coefficient. The use of a 2D model is justified since the estuary is very shallow (maximum depth around 19 m) and well mixed vertically. Moreover, the river flows are very low and a fast dispersion of fresh water into a much larger volume of salt water occurs (Borrego & Pendo ´n, 1988). The Odiel River flows usually range from less than 1 m 3 /s during the summer to some 70 m 3 /s in November. In the case of the Tinto river the corresponding flows are even smaller, with some 3 m 3 /s in November and no flow during the summer months. This mixing between river and sea water takes place upstream of the studied area, and for this reason horizontal gradients of salinity are not considered here. Indeed, electrical conductivity along the estuary is constant and typical of sea water, with an average value of 37.1 ± 0.5 mS/cm (Bolı ´var et al., 2000). The solution of these equations provides the instantaneous values of the two components of the current and the water elevation over the model domain, information required to solve the advection–diffusion dispersion equation of dissolved radionuclides. 4.2. Radionuclide equations The equation that gives the temporal evolution of activity concentrations in the dissolved phase, C d (Bq/m 3 ), is: oðHCdÞ otþoðuHCdÞ oxþoðvHCdÞ oy¼o oxHKD oCd ox  þo oyHKD oCd oy  k1CdHþk2AsLqsf/;ð6Þ where k 1 is given by Eqs. (1) and (2), total depth is H=D+z,K D is the diffusion coefficient, A s (Bq/kg) is activity concentration in the active sediment and q s is the sediment bulk density expressed in kg/m 3 . The external source of radionuclides should be added to this equation where necessary. The equation for the temporal evolution of activity concentration in the active sediment fraction is: oAs ot¼k1 CdH Lqsfk2As/:ð7Þ 4.3. Computational scheme The hydrodynamic equations are solved using an explicit finite difference scheme. The grid cell size is Dx=Dy= 125 m and the time step is fixed as Dt= 6 s. The CFL criterion is satisfied with these selections. Water elevations are specified for each time step along the southern boundary from observations. A radiation condition is ap42 R. Peria ´n˜ez / Marine Environmental Research 60 (2005) 35–49 plied along the northern and eastern open boundaries. Along the coast, the current component that is normal to the boundary is set to zero. Water depths were introduced for each grid cell from bathymetric maps. Instead of solving the hydrodynamic equations simultaneously with the dispersion equations, the hydrodynamic model is calibrated and validated in advance to speed up simulations. Once the hydrodynamics have been validated, standard tidal analysis is used to determine the tidal constants (amplitude and phase) for each grid cell. These constants are evaluated for both components of the flow (uand v) and for the water elevation (z), and for all the tidal constituents included in the model. Once the tidal constants are known, computation of flow and water elevation just involves the calculation and addition of a few cosine terms since the constants are stored in files that are read by the dispersion model. The net residual flow over the estuary must also be calculated by the hydrodynamic model and added to the instantaneous flow obtained from the tidal constants, since a residual transport cannot be generated with the pure harmonic currents that are given by the tidal analysis. The Monotonic Second Order Upstream (MSOU) explicit scheme is used to solve the advective transport in the dispersion equation of dissolved radionuclides. A second order accuracy scheme has also been used to solve the diffusion terms. It is considered that there is no flux of radionuclides through land boundaries. Along open boundaries, the boundary condition described in Peria ´n˜ez (1998) is applied. 5. Model results Only the two main tidal constituents, M 2 and S 2 , have been included. As will be shown below, this is enough to have a realistic representation of the dispersion patterns of 226 Ra. The calibration of the hydrodynamic model consisted of selecting the optimum value for the bed friction coefficient K. After some model runs, it was selected as K= 0.040 for the Tinto river and K= 0.005 for the rest of the estuary. With these selections, a reasonable agreement between observed and computed currents has been achieved. A comparison between observed and computed magnitude and direction of the maximum currents for several locations in the estuary (see Fig. 1)ispresented in Table 2 for a situation of medium tides (coefficient 74.4). Once that water circulation is reproduced by the hydrodynamic model, results are treated with standard tidal analysis to calculate the tidal constants to be used by the dispersion model to compute current and water elevation at any position and instant of time, as discussed previously. A value that depends on a horizontal length scale (set as the grid spacing) is chosen for the diffusion coefficient. Indeed, following Breton and Salomon (1995) such coefficient was taken as 0.61 m 2 /s. Since the dispersion model is not restricted by the CFL stability criterion, the time step in the dispersion model has been increased to 30 s. However, stability conditions imposed by the dispersion equation are satisfied with this value. R. Peria ´n˜ez / Marine Environmental Research 60 (2005) 35–49 43