Three-dimensional modelling of the tidal dispersion of non-conservative radionuclides in the marine environment. Application to 239,240 Pu dispersion in the eastern Irish Sea

Ž. Journal of Marine Systems 22 1999 37–51 Three-dimensional modelling of the tidal dispersion of non-conservative radionuclides in the marine environment. Application to 239,240Pu dispersion in the eastern Irish Sea R. Perianez ) ´˜ Dpto. Fısica Aplicada, E.U. Ingenierıa Tecnica Agrıcola, UniÕersidad de SeÕilla, Ctra. Utrera km 1, 41013-SeÕilla, Spain ´´´´ Received 15 July 1998; accepted 27 November 1998 Abstract A three-dimensional model to simulate the tidal dispersion of non-conservative radionuclides in a marine environment has been developed. The model solves the three-dimensional hydrodynamic equations using normalized s coordinates in the vertical and a flow dependent formulation for the eddy viscosity. Simultaneously, the suspended matter equation is also solved in s coordinates, which includes advective–diffusive transport of particles plus settling. Erosion and deposition are incorporated into the sea bed boundary condition of the equation. It is considered that radionuclides can be present in three phases: solution, suspended matter particles and bottom sediments. Three equations have been formulated, whose solutions give the time evolution of specific activity into each phase. The exchanges of radionuclides between the liquid and solid Ž. phase suspended matter and sediments have been described in terms of kinetic transfer coefficients, so that the model can be applied under non-equilibrium conditions. The model has been applied to simulate the dispersion of 239,240 Pu in the eastern Irish Sea, which is released from Sellafield nuclear fuel reprocessing plant. Computed and measured suspended matter concentrations have been compared for several points in the sea to test the suspended matter sub-model. Also, computed and observed 239,240 Pu distributions in water, suspended matter and sediments have been compared for several years. kdistribution coefficients have been computed as well. Model results are, in general, in good agreement with d observations. q1999 Elsevier Science B.V. All rights reserved. Keywords: s -coordinates; tides; sediments; suspended matter; transfer coefficients 1. Introduction Over the past years, there has been an increasing interest in developing mathematical models to simulate the dispersion of radionuclides in aquatic envi- )Tel.: q34-95423-3669; Fax: q34-95423-2644; E-mail: [email protected] ronments since these models can be used as predictive instruments that can be applied in the assessment of radioactive contamination following an accidental release of radionuclides. On the other hand, Ž useful oceanographic information for instance, . flushing and transit times can be obtained by means Ž of the application of such models Prandle, 1984; . Salomon et al., 1995 . 0924-7963r99r$ - see front matter q1999 Elsevier Science B.V. All rights reserved. Ž. PII: S0924-7963 99 00029-9 () R. PerianezrJournal of Marine Systems 22 1999 37–51 ´˜ 38 Radionuclides are said to behave conservatively in water if they remain in solution, that is, no fraction is removed from the water column due to geochemical processes. In this case, radionuclide transport is governed by an advection diffusion equation. When radionuclides have a certain affinity to be Ž fixed to the solid phase suspended matter particles . and bottom sediments , they are said to be non-conservative. Thus, adsorption–desorption reactions must be included in models developed to simulate the dispersion of non-conservative radionuclides. On the other hand, deposition of suspended particles and erosion of the sediment will produce an exchange of radionuclides between suspended matter and the sediment. Thus, the suspended matter dynamic must be solved also. A detailed discussion of these processes Ž. can be found in Perianez et al. 1996a . ´˜ The first models developed to simulate the dispersion of non-conservative radionuclides in the sea Ž were long-term dispersion models Gurbutt et al., . 1987; Abril and Garcıa-Leon, 1993a in which the ´´ time steps used allowed to describe the transfers of radionuclides between the liquid and solid phases in terms of equilibrium distribution coefficients, k. d Ž. The same approach was used by Ng et al. 1996 to simulate the dispersion of Cd and Zn in an English Ž. estuary. Perianez et al. 1996a developed a mathe- ´˜ matical model to simulate the dispersion of non-conservative radionuclides in tidal waters based upon kinetic transfer coefficients, thus the model could be used under non-equilibrium conditions, and applied this model to study the dispersion of 226Ra,238U and 232 Ž Th in a Spanish estuary Perianez and Martınez- ´˜ ´ . Aguirre, 1997a; Perianez et al., 1996b . ´˜ Ž All these models are two-dimensional integrated .Ž. in the vertical . Prandle et al. 1993 have pointed out that vertical structure can be of primary concern, Ž. even in shallow water depth ;50 m , when typical simulated times are in the order of one month or less if tidal mixing keeps the vertical diffusion coefficient y32y1Ž smaller than 10 m s this value can be consid- . ered representative of a strong tidal action . However, few three-dimensional dispersion models have Ž. been developed. Nies et al. 1997 and Baxter et al. Ž. 1998 studied the three-dimensional transport of radionuclides in the Arctic Ocean, but the authors could not compare the model results with observaŽ. tions. Perianez 1998a developed a three-dimen- ´˜ sional model of the Irish Sea to study the three-dimensional dispersion of dissolved 137Cs, released from Sellafield nuclear fuel reprocessing plant. Some interesting qualitative results were obtained from this work: for instance, it was found an alternating, with the same frequency as tidal oscillations, of wellmixed and vertically structured states. This model was improved using normalized s -coordinates in the vertical and a flow dependent eddy viscosity Ž. 137 Perianez, in press . Observed and computed Cs ´˜ distributions in the eastern Irish Sea were compared for a number of years. Model results were, in general, in good agreement with observations. These models can be used for conservative radionuclides. The objective of this paper is to present a threedimensional model for non-conservative radionuclides dispersion in tidal waters. The model solves the three-dimensional hydrodynamic equations in s -coordinates and with a flow dependent eddy viscosity. Simultaneously, it solves the suspended matter dynamic equation, which is also written in s -coordinates and includes the deposition and erosion terms in the sea bed boundary condition. Finally, the equations for radionuclides are solved. The model includes three phases: radionuclides in solution, in suspended matter and in the active fraction of the Ž. sediment as denoted by Benes et al., 1994 , that is, Ž in muds particles with a diameter -62.5 mm that can be eroded and incorporated into the water col- . umn as suspended matter . Thus, three equations more are solved. These equations include advective Ž and diffusive transport in the case of radionuclides . in solution and in suspended matter plus the transfer of radionuclides between phases, which is formulated using kinetic transfer coefficients. The model has been applied to simulate the dispersion of 239,240Pu in the eastern Irish Sea, which is released from Sellafield nuclear fuel reprocessing plant. Observed and computed plutonium distributions in water, suspended matter and sediments have been compared. Also, kvalues over the sea and kprofiles dd in the water column have been obtained with the model. The model equations are presented in the Section 2, where the methods used to solve them are also briefly described. After, the values for all the parameters involved in the model are given and finally model results are presented and discussed. () R. PerianezrJournal of Marine Systems 22 1999 37–51 ´˜ 39 2. The model 2.1. Hydrodynamics The three-dimensional hydrodynamic equations are written using normalized s coordinates in the vertical. This way a constant number of grid boxes is used in the vertical at each horizontal grid point and vertical resolution is not lost in the shallower areas. Ž The transformation to s coordinates is see Davies, . 1985a, for instance : zq z s s1 Ž. hq z Ž. where his the undisturbed mean depth of water, z is sea surface displacement from the mean level due to tidal oscillations and zcoordinate is measured from the mean sea level to the bottom. Thus, the hydrodynamic equations are transformed from the interval — z FzFhinto the constant interval 0F s F1. The form of the equations, for a homogeneous sea, in normalized coordinates has been given previŽ. ously Davies, 1985a; Davies and Stephens, 1983 and will not be repeated here. The boundary conditions applied at the sea bottom and sea surface can Ž. also be seen in Perianez in press . Wind effects are ´˜ not considered and a linear law for bottom friction has been used. Although a quadratic law for bottom friction is now more extended than a linear formulation, the linear law has been adopted since a faster convergence of the equations is obtained. Also, there are no appreciable differences in dispersion patterns Ž. for a dissolved radionuclide obtained when the hydrodynamic model is calibrated using a linear and Ž a quadratic law for bottom friction Perianez, in ´˜ . press . A flow-dependent eddy viscosity has been used in the model. This formulation has been used previously and has given good results for tidal flow Ž studies Davies and Lawrence, 1994; Jones and . Davies, 1996; Davies et al., 1997a : 22 ( NsCuqÕh2 Ž. N where Cs0.0025 is a dimensionless experimenN tally measured coefficient and uand Õare depth mean currents along the xand yaxis, respectively. Eddy viscosity has been taken constant in the vertical. This approximation has also been used by Davies Ž. Ž . and Lawrence 1994 and Davies et al. 1997b . 2.2. Suspended matter dynamics As usual in suspended matter studies, we will consider that only particles with a diameter -62.5 mm can remain in the water column as suspended matter, since larger particles will sink rapidly to the bottom and, as a consequence, their horizontal moveŽ ment is negligible Belderson, 1964; Gurbutt et al., 1987; Abril and Garcıa-Leon, 1994; Clarke, 1995; ´´ .Ž. Perianez et al., 1996c . Indeed, Eisma 1981 has ´˜ pointed out that for all practical purposes muds can be regarded as synonymous of suspended matter. Although bed load transport of coarse particles may exist, it is not included in the model since, as will be shown below, specific activity in large sediment particles can be neglected when compared to specific activity in the active sediment. The three-dimensional suspended matter equation has been obtained by transforming the equation proŽ. posed by Nicholson and O’Connor 1986 into s coordinates. Such equation includes advective–diffuŽ. sive transport plus vertical fall settling of suspended particles: E m E m E m E m quqÕqw) E t E x E y Es EE m EE m sKqK hh ž/ ž/ E x E x E y E y 1 EE m qKv 2 ž/ EsEs hq z Ž. 1 E wm Ž. s y3 Ž. hq zEs where mis the suspended matter concentration in Ž. ppm parts per million , w)is the vertical water velocity along the s axis, Kand Kare the hv horizontal and vertical diffusion coefficients respectively and wis the settling velocity. s The vertical diffusion coefficient can be written as Ž a function of the eddy viscosity Launder and Spald- . ing, 1972 : Ks e N4 Ž. v () R. PerianezrJournal of Marine Systems 22 1999 37–51 ´˜ 40 where the non-dimensional number ´ ranges from 0.1 to 0.5. It is assumed that there is no flux of particles through the sea surface. Deposition and erosion of the sediment are incorporated into the sea bed Ž. boundary condition Nicholson and O’Connor, 1986 : 1 E mM KsEfq ywbmb Ž. Ž. vs ž/ hq zEs s s1 =q 1y5 Ž. ž/ qcd The first term of the right member of the equation Ž represents erosion: Eis the erosion constant Nichol- . son and O’Connor, 1986 that also acts as a scaling factor, Mis some power of the velocity typically in Ž. the range 2–5 Prandle, 1997 and fgives the fraction of active sediments. It is considered that erosion occurs only when the near-bed current mag22 ( nitude, qsubqÕbis larger than a critical Ž. Ž. erosion velocity, q. Otherwise this term is set to ce zero. The second term represents deposition of partiŽ. cles on the sea bed: wbis the near-bed settling s velocity and qis a critical deposition velocity. cd Thus, deposition occurs only if q-q; otherwise cd this term is set to zero. Ž. Note that bmeans that the corresponding magnitude is evaluated in the water layer that is in contact with the sediment. The settling velocity of suspended particles, at Ž. low concentration smaller than 1000 ppm , increases as the suspended matter concentration increases. A standard formula used to represent this Ž process is Pejrup, 1988; Mehta, 1989; Nicholson . and O’Connor, 1986; Clarke, 1995 : wsam a26 Ž. s1 It has been taken as1.7=10y6and as1.6 if 12 wis measured in m sy1and min ppm. Since s suspended matter concentrations in the Irish Sea are of the order of 1 ppm, settling velocities of the order of 10y6ms y1are obtained, which is the order of Ž magnitude that can be found in literature Prandle et . al., 1993; Prandle, 1997 for the settling velocity of mud particles. Although the accepted interval for a2 Ž. is 0.5 – 1.2 Pejrup, 1988 , a value of 1.6 was Ž. obtained by Mehta 1989 and was also used in Ž. previous modelling work Perianez et al., 1996c . It ´˜ can be found in literature that aranges from ;10y3 1 Ž. y9Ž. Clarke, 1995 to ;10 Mehta, 1989 . A value of 0.18 m sy1has been used for the Ž critical deposition velocity Eisma, 1993; Perianez et ´˜ .y1 al., 1996c . A typical value of 0.21 m s was taken Ž for the critical erosion velocity Heathershaw, 1981; . Pugh, 1987 . A source term must be included to the suspended matter equation in those grid cells located along the coastline. This term represents the input of particles Ž from runoff of continental waters Abril and Garcıa- ´ . Leon, 1994 and can be written as: ´ a I7 Ž. i,j where a is a constant to be fixed that represents the mass of suspended load per continental water volume unit and Iis the total volume of continental i,j water that enters compartment i,j.Iwill be differi,j ent from zero when the compartment is situated along the coastline. The solution of the suspended matter equation provides the suspended matter concentration at each position in the model domain and at each instant of time, but also gives information about the sedimentation processes that take place in the area under study. Thus, the sedimentation rate is obtained as the net balance between the deposition and erosion terms. Ž. It has been shown Davies and Lawrence, 1995 that wave-current interaction must be considered when studying sediment transport, and this is particuŽ larly important during major wind events when . wave heights are large and swell is significant . In this situation suspended matter concentration will Ž increase and the near bed current magnitude influ- . enced by wave-current interaction will be the primary mechanism moving the suspended sediment. However, at this stage of model development, wave-current interaction has not been considered. Nevertheless, the model reproduces historically deduced net sedimentation rates, and realistic plutonium concentrations in water, suspended matter and sediments are given by the model, as will be shown. It is our objective to study, in the next future, the influence of wind in the dispersion patterns of ra- () R. PerianezrJournal of Marine Systems 22 1999 37–51 ´˜ 41 dionuclides. Then wave-current interaction will have to be considered. 2.3. Radionuclides dispersion As was commented in Section 1, the model considers that radionuclides can be present in three phases: solution, suspended matter and active bottom Ž. sediments particles with a diameter -62.5 mm.In some preliminary calculations, the larger grain size Ž. sediment particles with diameter )62.5 mm was Ž. also considered, as in Perianez et al. 1996a . How- ´˜ ever, it was found that specific activity in the large grain size fraction of the sediment could be neglected when compared to specific activity in the active sediment. Thus, only the active sediment was included in the final model. This is not in agreement with the calculations of Abril and Garcıa-Leon ´´ Ž. 1993b , that found that specific activities in both sediment fractions are comparable. Nevertheless, our model confirms previous experimental work: Aston Ž. et al. 1985 found that virtually all the activity is associated with the active sediment, and McKay and Ž. Walker 1990 found that Pu specific activity in the active sediment is of the order of 20 times greater than that of the sand. We will consider that the transfer of radionuclides Ž from water to the solid phase suspended matter and . bottom sediments is governed by a coefficient k1 and the inverse process by a coefficient k, which 2 are denoted kinetic transfer coefficients. It is known that actinide adsorption tends to be a surface pheŽ. nomenon Ramsay and Raw, 1987 and will depend on the surface of particles per water volume unit into the grid cell. This quantity has been denoted as the Ž exchange surface Perianez and Martınez-Aguirre, ´˜ ´ .Ž 1997a; Perianez et al., 1996a . Thus Perianez et al., ´˜ ´˜ . 1996a : ks x SqS8 Ž. Ž. 11ms where Sand Sare the exchange surfaces for ms suspended matter and bottom sediments, respectively, and x is a parameter with the dimensions of 1 a velocity. It is denoted as the exchange velocity Ž. Perianez et al., 1996a . As a first approach, assum- ´˜ ing spherical particles and a step function for the grain size distribution of particles, it can be easily Ž. obtained Perianez et al., 1996a that: ´˜ 3m Ss9 Ž. m r R 3Lf f Ss10 Ž. sR Ds H where r is the suspended particles density, Ris the mean radius of suspended matter and active sediment particles, Hshq z is total water depth, Lis the Ž average mixing depth the distance to which the . dissolved phase penetrates the sediment , D s is spacing in the vertical direction and f is a correction factor that takes into account that not all the mass of the sediment is in contact with water. Then the kinetic transfer coefficient kis written as: 1 3m3Lf f ks x q x skqk11 Ž. 11 1 1112 r RR Ds H The transfer coefficient kis considered constant. 2 This description of the transfer of radionuclides between the dissolved and solid phases has been used Ž successfully in previous modelling studies Perianez, ´˜ 1998b; Perianez and Martınez-Aguirre, 1997a,b; Pe- ´˜ ´ rianez et al., 1996b; Martınez-Aguirre and Perianez, ´˜ ´ ´˜ . 1998 . Thus, the equation which gives the time evolution of radionuclide concentration in the dissolved phase Žy3. C, in Bq m in each grid cell is: d E Cdsadvqdif ykq p kC Ž.Ž . 3D 11 12 d E t AL r f f sm 3 qkmCq p =10 y l C 2s d. ž/ Ds H 12 Ž. The time derivative is interpreted as in a forward Ž. finite differences scheme. advqdif represents 3D the three-dimensional advective plus diffusive transport of dissolved radionuclides from one grid cell to another. p s0 unless we are solving the equation for a grid cell that is in contact with the sediment; in this case p s1 to allow the transfer of radionuclides between water and the bottom sediment. Cis the s concentration of radionuclides in the suspended matŽy1. ter particles Bq g , Ais the concentration of sŽy1. radionuclides in the active sediment Bq g and the sediment bulk density, r , is expressed in kg m () R. PerianezrJournal of Marine Systems 22 1999 37–51 ´˜ 42 my3. Finally, l is the radioactive decay constant. The external source of radionuclides should be added to this equation if it exists. The equation for the time evolution of specific Žy1. activity in suspended matter Bq g is: E CC sd sadvqdif ysettqk Ž. 3D 11 E tm) mm ykC q p resydep y l C Ž. 2s s m)m) 13 Ž. Ž. where sett represents settling see Eq. 3 and the ) means that the corresponding suspended matter concentration must be evaluated in the new time step. Žy1y1. The deposition term Bq g s is written as: wbCbmb q Ž. Ž. Ž. ss deps1y.14 Ž. ž/ H Ds m)bq Ž. cd Žy1y1. Similarly, the resuspension term Bq g s is written as: fEqMAs ress.15 Ž. Ds Hm)b Ž. Details about how the expressions for these terms Ž can be deduced in the case of a two-dimensional .Ž. model can be seen in Perianez et al. 1996a . ´˜ The equation for the time evolution of specific Ž activity in the active fraction of the sediment, ABq s y1. g is: E ACbH Ds Ž. sd y3 sk=10 ykA f 12 2 s E tL r f m qdepyres y l A.16 Ž. Ž. s The deposition and erosion terms are now written as follows: wbmbCb q Ž. Ž. Ž. ss y3 deps1y=10 17 Ž. ž/ L r fq mcd EqMAsy3 ress=10 . 18 Ž. L r m Again, details about how these terms are obtained, in the case of a two-dimensional model, can be seen in Ž. Perianez et al. 1996a . ´˜ 2.4. Numerical solution The hydrodynamic equations are solved using a explicit finite difference scheme, although the Saul’ev Ž. 1957 implicit method is applied to the vertical Ž. diffusion term to retain stability Davies, 1985b . Some boundary conditions must also be given. At land boundaries, the normal component of the current and the flux of suspended particles and radionuclides are zero. Along open boundaries, amplitude and phase of the surface elevation at the Mfre2 quency are specified from observations and the surface current component that is normal to the boundŽ ary is obtained from a radiation condition Kowalick . and Murty, 1993; Glorioso and Davies, 1995 . Of course, better results would be obtained if other tidal Ž. components at least Sand Nare included, since 22 they induce currents that may affect the advection of radionuclides. However, we are mainly interested in tuning and testing the three-dimensional formulation of the dispersion of non conservative radionuclides: as will be shown below, realistic results are obtained although only the Mtide is used. This is also the 2 reason why surface wind stress is not considered in the model. On the other hand, inclusion of the S2 and N tides would influence bed stress, and thus, 2 sediment resuspension. However, as will be shown, slight variations in computed suspended matter concentrations are obtained when Eand Mare changed. This suggests that suspended matter concentrations in the eastern Irish Sea are not dominated by erosion events. Thus, it is not expected that the increase in bed stress produced by these other components produce a significant change in erosion and in suspended matter concentrations. The open boundary condition described in Ž. Perianez 1998a, in press is applied to the dispersion ´˜ equation of suspended matter and radionuclides: Cs b C19 Ž. iiy1 where Cis the concentration in the open boundary i and Crepresents the concentration just inside the iy1 computational domain. The non-dimensional number b is obtained from a calibration exercise. The advection and horizontal diffusion terms are treated using explicit second order accuracy finite () R. PerianezrJournal of Marine Systems 22 1999 37–51 ´˜ 43 Ž. difference schemes Kowalick and Murty, 1993 and vertical diffusion is again solved using the Saul’ev method. AFORTRAN code was developed to solve the equations involved in the model. The code was implemented on a HP SPP 2000 X-Class computer. It has 24 PA-8000 processors at 180 MHz rate. Memory is 3 GB. The equations are solved, for each time step, in the following sequence. Ž. 1 Hydrodynamic equations are solved to obtain the water elevation and components of the water velocity for each grid cell. Eddy viscosity is also computed. Ž. 2 Suspended matter equation is solved to obtain the suspended matter concentration into each grid cell. Deposition and erosion are also evaluated. Ž. 3 Radionuclide dispersion: Equations for dissolved radionuclides, radionuclides in suspended matter and radionuclides in the active sediments are solved in this sequence. Ž. 4 Sources: External sources of suspended matter, dissolved radionuclides and radionuclides in particulate form are introduced into the grid cells in which these sources exist. One month of simulation takes about 2.5 h of CPU time. Model results will be compared with observations. 3. Application of the model The model has been applied to study the dispersion of 239,240 Pu in the eastern Irish Sea. These radionuclides are discharged from a nuclear fuel reprocessing plant at Sellafield. The model has a horizontal resolution DxsDys5000 m. Ten layers are used in the vertical direction, thus, D s s0.1. Time step is fixed as Dts60 s. Stability conditions Ž. Kowalick and Murty, 1993 are satisfied with this selection. The computational domain is presented in Fig. 1, where the location of Sellafield plant is also shown. Water depths are introduced from bathymetric maps and range from 55 m in the west of the computational domain to a shallower area around the Ž. British coast. The fraction of active sediments muds Ž. has been obtained from Pentreath 1985 and Abril Ž. and Garcıa-Leon 1994 . ´´ Fig. 1. Map of the computational domain. The star is Sellafield nuclear fuel reprocessing plant and circles are the points where Ž. suspended matter concentrations were measured see Section 4 . Ž. Each unit in the xand yaxis is 5000 m grid element number . Surface elevations and phases along the open boundary have been taken from the observations of Ž. Howarth 1990 . The mean value 0.3 has been taken for ´ in Eq. 4. The horizontal diffusion coefficient has been fixed as Ks500 m2sy1since results in good agreement h with observations are obtained with this value. InŽ. deed, Bowden 1950 suggested that for the Irish Sea 500-K-900 m2sy1and, on the other hand, this h value has been successfully used in previous modŽ. elling work Perianez, in press . ´˜ After some model runs, a value of 2.7=10y4 y11 2 5 2 gm s was selected for the erosion constant Ž and it was taken Ms3.5. On the other hand a Eq. .6y3 7 was fixed as 25=10 kg km . The input of continental waters along the coast was obtained from Ž. MAFF 1987 , and is also presented in Abril and Ž. Garcıa-Leon 1994 . Results in agreement with ob- ´´ servations are obtained with these selections. The mixing depth Lwas taken a 0.1 m, following Ž previous modelling works Abril and Garcıa-Leon, ´´ 1993b; Perianez and Martınez-Aguirre, 1997b; ´˜ ´ . Perianez et al., 1996b . The density of suspended ´˜ matter particles is taken as r s2600 kg my3, which () R. PerianezrJournal of Marine Systems 22 1999 37–51 ´˜ 44 is the established value for soil particles density Ž. Baver et al., 1972 . The mean radius of suspended matter and active sediment particles has been selected as Rs15 mm, since only particles with a diameter -62.5 mm can be present in the water Ž column as suspended matter Gurbutt et al., 1987; .y3 Perianez et al., 1996b . A mean value of 950 kg m ´˜ Ž was selected for the sediment bulk density Kershaw . and Young, 1988; Abril and Garcıa-Leon, 1994 . ´´ After a calibration exercise, the geometry correction factor was taken as f s0.1 and a value of 0.9 was Ž selected for b in equation 19 Perianez, 1998a, in ´˜ . press . The source of virtually all plutonium to the Irish Sea is the discharge from Sellafield through a pipeline which extends 2.5 km beyond high water. The annual discharges are presented in Fig. 2, compiled Ž. from McKay and Pattenden 1993 and Hunt et al. Ž. 1997 . It has also been reported that about 99% of the discharged Pu is associated with particulate matŽ. ter Pentreath, 1985 . Finally, to simulate the dispersion of 239,240Pu the values of the exchange velocity, x , and of the 1 kinetic transfer coefficient, k, must be known for 2 this radionuclide. These parameters can be obtained from the measurements carried out by Nyffeler et al. Ž. 1984 and the equilibrium distribution coefficient, k, following the method presented in detail in dŽ. Perianez and Martınez-Aguirre 1997a and Perianez ´˜ ´ ´˜ Ž. 1998b . However, the behaviour of plutonium in aquatic systems is of considerable complexity due to the fact that, in solution, it can exist in different Ž. oxidation states simultaneously. Thus, Pu III and 239,240 Žy1. Fig. 2. Annual Pu discharges Bq year to the sea from Sellafield. Ž. Ž. Pu IV predominate as the reduced and Pu V and Ž. Pu VI as the oxidized forms. The reduced Pu is highly particle-reactive and has been shown to possess a distribution coefficient which is approximately two orders of magnitude higher than that of the more Ž soluble oxidized Pu McKay and Pattenden, 1993; . Mitchell et al., 1995 . Hence, the values observed in field measurements represent the properties of the mixture of oxidation states that is present in the Ž. particular sample. Pentreath et al. 1986 have found that oxidation of Pu occurs quickly following release from Sellafield to sea water, and Mitchell et al. Ž. 1995 report that the oxidation state of Pu is the same all over the Irish Sea, result that has also been Ž. confirmed by Boust et al. 1996 . Thus, since there are no differences in the speciation of Pu over the sea, a constant value for the equilibrium distribution coefficient can be used. We have taken the recommended value for Pu in coastal waters, 1.0=105L y1Ž. kg IAEA, 1985 , since this value is in agreement Ž with measurements in the Irish Sea McKay and Walker, 1990; Garland et al., 1990; Mitchell et al., . 1995 . As commented above, from the value of the Pu distribution coefficient and the experiments of Ž. Ž Nyffeler et al. 1984 , it can be obtained Perianez, ´˜ . 1998b that: x s1.51=10y5ms y1 1 ks1.16=10y5sy1 2 This values have been successfully used in previŽ. ous modelling work Perianez, 1998b . ´˜ 4. Results and discussion The hydrodynamic part of the model was tested by comparing observed and computed values of tidal elevation amplitudes and phases, current profiles at several locations and semi-major axis magnitude and orientation of the Mtidal current ellipse at several 2 depths and locations. The difference between computed and observed tidal amplitudes is always - Ž 10%, only in one point is 17.4% comparisons have . been made for 12 points . In the case of tidal phases, absolute values of the differences between observed and computed phases range from 08to 278, being the mean value 10.58. The difference between computed () R. PerianezrJournal of Marine Systems 22 1999 37–51 ´˜ 45 and observed semi-major axis magnitude of the tidal current ellipse ranges from 1.5% to 45%. Only in Ž three points errors are above 25% 12 points used for . comparison . The absolute value of the difference between observed and computed orientation of the axis ranges from 0.18to 30.58, being the mean value Ž. 7.68. Details can be seen in Perianez in press . In ´˜ that work, the advection–diffusion dispersion equation was calibrated by studying the dispersion of dissolved 137Cs. The calibration consists of selecting the optimum boundary conditions and values of the diffusion coefficients, which were presented in Section 3. Ž. Kershaw and Young 1988 measured the suspended matter concentration at several locations in the eastern Irish Sea and for both surface and sea bottom samples. Sampling points are shown in Fig. 1. Computed suspended matter profiles, along the west-to-east direction, following the points of KerŽ. shaw and Young 1988 for surface and bottom waters are presented in Fig. 3a,b, respectively, toFig. 3. Measured and computed suspended matter concentration Ž. Ž . ppm profiles following the points of Kershaw and Young 1988 , Ž. Ž. shown in Fig. 1, in surface a and bottom b waters. Ž. Fig. 4. Computed suspended matter ppm distribution in surface Ž waters. Each unit in the xand yaxis is 5000 m grid element . number . gether with the measured concentrations. As can be seen, the model gives a good estimation of the observed suspended matter concentrations for both set of samples. As expected, maximum concentrations are obtained close to the coast, and these concentrations decrease when entering into the open sea. A map showing the computed surface suspended matter concentration is presented in Fig. 4, which is Ž in agreement with previous modelling work Abril . and Garcıa-Leon, 1994 . Averaged sedimentation ´´ Žy2 rates are small range from y0.02 to 0.05 g cm y1. year , thus the model does not violate historically Ž. deduced net sedimentation rates: Kirby et al. 1983 suggest that no major erosion or accretion of sediments is taking place in the Irish Sea at the present time, so that a state of equilibrium exists. Kershaw et Ž. al. 1988 have also found very low sedimentation Žy2y2y1. rates of the order of 10 g cm year in the eastern Irish Sea. Observed and computed distributions of 239,240Pu in water, suspended matter and bottom sediments Ž have been compared for several years 1974, 1977, .Ž. 1979 and 1989 . The input from Sellafield Fig. 2 was introduced in the model for each year. However, as can be seen in Fig. 2, the input has been taking place since the fifties. Thus, instead of starting the model from zero concentrations, we have assumed a uniform background in water, suspended matter and sediments that represents the effect of previous dis-