A three dimensional σ-coordinate model to simulate the dispersion of radionuclides in the marine environment: application to the Irish Sea

A three dimensional s-coordinate model to simulate the dispersion of radionuclides in the marine environment: application to the Irish Sea R. Peria´n˜ez * Departamento Fı´sica Aplicada, E.U. Ingenierı´a Tecnica Agrı´cola, Uni6ersidad de Se6illa, Ctra. utrera Km 1 , 41013 Se6illa, Spain Abstract A three dimensional model to simulate the tide induced dispersion of radionuclides in the sea has been developed. The model uses normalized s-coordinates in the vertical so that resolution is not reduced in the shallower regions. The hydrodynamic equations are solved and, simultaneously, the three dimensional advection diffusion dispersion equation (also written in normalized coordinates) is solved too. An instantaneous flow and depth dependent eddy viscosity has been used. The model has been applied to study the dispersion of l37 Cs in the eastern Irish Sea, where a nuclear fuel reprocessing plant releases radionuclides. The hydrodynamic part of the model has been tested by comparing observed and computed values of tidal elevations, phases and currents. The model gives, in general, a good representation of the water circulation in the sea. Also, it gives results in agreement with observations when measured and computed levels of l37 Cs are compared. In order to show how the model can be used to obtain water quality parameters of interest, it has been applied to obtain the turn over time of a region of the sea. It has also been applied to simulate the dispersion of a pollutant after an hypothetical accidental discharge. Keywords : Tides; Currents; Advection; Diffusion; 137 Cs 1. Introduction Radionuclides are being released to the sea from nuclear fuel reprocessing plants (McKay and Baxter, 1985; Guegueniat et al., 1994; Herrmann et al., 1995; Cook et al., 1997). There has been an increasing interest in developing and improving models to simulate the dispersion of these radionuclides in the marine environment (Prandle, 1984; Breton and Salomon, 1995) due to the fact that useful oceanographic information (like flushing and transit times) can be obtained by means * Tel.: +34 95 4233669; fax: +34 95 4232644. of the application of the models (Prandle, 1984; Salomon et al., 1995). On the other hand, these models can be used as predictive instruments which can be applied in the assessment of radioactive contamination following an accidental release of radionuclides. Although these models are two dimensional, 2D (integrated in the vertical direction), there has been an increasing use of three dimensional models (3D) to solve oceanographic problems (Davies and Lawrence, 1994; Proctor and James, 1996; Davies et al., 1997) due to the limitations of 2D models, which for instance, can not give information about the vertical profile of currents. A 3D model can also be applied to simulate the dispersion of pollutants in the sea. Indeed, Prandle et al. (1993) have shown that vertical structure can be of primary concern when typical simulated times are in the order of one month or less. Moreover, the vertical variability is important for such typical times, even in shallow waters (depth :50 m), if tidal mixing keeps the vertical diffusion coefficient smaller than 10 −3 m 2 /s (Prandle et al., 1993) (a value in the order of 10 −3 m 2 /s can be considered representative of a strong tidal action). Thus, a 3D approach should be used to build models that can be used as predictive tools in the assessment of the near field radioactive contamination of tidal waters, since the time scale involved is of the order of weeks and vertical structure can be important, as said above. Only in the case of a very strong tidal action a 2D approach can give enough accurate results. However, few 3D dispersion models have been developed: the model of Nies et al. (1997), for instance, studies the 3D transport of radionuclides in the Arctic Ocean, but the authors could not compare the model predictions with observations. In Peria´n˜ez (1998) some preliminary results on a 3D dispersion model for radionuclides were presented. However, this model worked with layers (in the vertical) of the same thickness, which implies a loss of resolution in the vertical direction in the shallower regions. On the other hand, the model used a very simple parametrization of the eddy viscosity coefficient (it was considered constant in time and with the same value over the whole computational domain). This work focused on the study of some qualitative aspects of three dimensional tidal dispersion. Now the model has been improved to overcome the above mentioned problems. Thus, normalized scoordinates are used in the vertical direction, in such a way that resolution is not lost in the shallower regions. On the other hand, eddy viscosity has been formulated in such a way that it depends on the instantaneous current and depth. Thus, a more realistic model has been developed. The model solves the 3D hydrodynamic equations and, simultaneously, the advection diffusion dispersion equation, which has also been written in scoordinates. This way, the tide induced dispersion of radionuclides is obtained. An implicit numerical scheme has been adopted to solve the hydrodynamic and dispersion equations in the vertical direction in order to retain stability. This is the first time, to the author’s knowledge, that scoordinates are used to solve the dispersion equation. The model has been applied to the Irish Sea since there is enough oceanographic information to test the model results (tide amplitudes, currents and current profiles), there is a well known source of radioactivity (Sellafield nuclear fuel reprocessing plant) and there are measurements of radionuclide concentrations over the sea. The hydrodynamic part of the model has been tested by comparing observed and computed tidal amplitudes, tidal phases and magnitudes and directions of tidal currents. Observed and computed 137 Cs concentrations have also been compared. The turn over time of an estuary of the Irish Sea has been calculated to show how the model can be applied to obtain water quality parameters of interest. Finally, a numerical experiment has been carried out to show the potential predictive power of the model in the assessment of contamination following a hypothetical accident. The model equations are presented in the next section. Next the numerical methods used to solve them are described briefly and finally model results are presented and discussed. 2. Model equations If a fixed finite difference grid in the vertical is used, the number of vertical grid boxes decreases in the shallow regions. This has the effect of reducing the vertical resolution in the shallow water areas. This problem can be solved transforming the 3D hydrodynamic equations into depth following scoordinates. This way a constant number of grid boxes is used in the vertical at each horizontal grid point. The transformation to scoordinates is (see Davies, 1985a, for instance): s=z+z h+z(1) where his the undisturbed (mean) depth of water, zis sea surface displacement from the mean level due to tidal oscillations and zcoordinate is measured from the mean sea level to the sea bottom. Thus, the hydrodynamic equations are transformed from the interval −z5z5hinto the constant interval 05s51. The transformed equations for an incompressible flow and for a homogeneous sea are (Davies, 1985a): (z (t+( (x(h+z)& 0 1 udsn +( (y(h+z)& 0 1 ndsn=0 (2) (u (t+g(z (x−V6=1 (h+z) 2 ( (sN(u (s(3) (6 (t+g(z (y+Vu=1 (h+z) 2 ( (sN(6 (s(4) w*=1 h+z(z (t(1−s)n +1 h+z ( (x(h+z)& s 1 udsn +1 h+z ( (y(h+z)& s 1 6dsn(5) where u,6and w* are the components of the water velocity along the directions of x,yand s axis, respectively. gis gravity, Vis the Coriolis parameter, V=2vsin f;(vbeing the earth rotational angular velocity and fthe latitude), and N is the coefficient of eddy viscosity. The non linear advective terms have been removed from these equations since dimensional analysis has shown (Charnock and Crease, 1957) that they are important only when zis comparable with the mean depth, and this is not the case. These non linear terms generate a residual current that may affect the radionuclide transport when studying longterm dispersion (time scale of the order of years). However, this weak current can be neglected if time scale of interest is in the order of weeks, as is the case. Moreover, Prandle (1984) excludes the advective terms in his long-term dispersion model for the European shelf seas since their effect is to add additional structure to residual distributions, but this structure is often exaggerated due to poor topographic resolution (Prandle, 1984). To solve the equations, surface and sea bed boundary conditions must be specified. The surface boundary conditions are: rN(u (s s=0 =−(h+z)F s (6) rN(6 (s s=0 =−(h+z)G s (7) where F s and G s denote the components of wind stress acting on the water surface along the xand ydirections, which can be written as in Pugh (1987) and Peria´n˜ez et al. (1994), and ris the water density. Similarly, at the sea bed: rN(u (s s=1 =−(h+z)F b (8) rN(6 (s s=1 =−(h+z)G b (9) where F b and G b are the components of bottom stress. Assuming a linear law of bottom friction: F b =rku b (10) G b =rk6 b (11) where kis a friction coefficient and u b and 6 b are the components of water velocity at a given height above the bottom, which usually is 1 m (Davies and Stephens, 1983). 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. Indeed, a linear law is more appropriate in linear models (Davies, 1985b). The 3D advection-diffusion dispersion equation for dissolved radionuclides, which has been written in scoordinates, is: (C (t+u(C (x+6(C (y+w* (C (s =( (xK x (C (x+( (yK y (C (y +1 (h+z) 2 ( (sK 6 (C (s−lC(12) where Cis the radionuclide concentration, lthe radioactive decay constant and K x ,K y and K 6 are the diffusion coefficients along the x,yand vertical directions, respectively. The external source of radionuclides, where it exists, should be included in this equation. 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., 1997) N=C N u¯ 2 +6¯ 2 h(13) where C N =0.0025 is a dimensionless experimentally measured coefficient and u¯ and 6¯are depth mean currents: u¯ =& 0 1 uds(14) 6¯=& 0 1 6ds(15) Eddy viscosity decreases in the region close to the sea bed (Davies et al., 1997). However, it has been taken constant in the vertical since no attempt has been made to solve the high shear region close to the sea bed. This approximation has also been used by Davies and Lawrence (1994). The vertical diffusion coefficient can be written as a function of the eddy viscosity (Kowalick and Murty, 1993): K 6 =oN(16) where the non-dimensional number oranges from 0.1 to 0.5. 3. Numerical solution All the equations are solved using finite differences. A staggered grid is used in the horizontal, with uniform spacing Dxand Dy, and a grid in s coordinates, with spacing Ds, is employed in the vertical direction. Time step is fixed as Dt. The model is started from rest. If the hydrodynamic equations are integrated using an explicit method besides the CFL criterion, there is another stability condition imposed by the vertical diffusion term, which is related to the magnitude of eddy viscosity and the water depth: DtB(hDs) 2 2N(17) Thus, in shallow water, this condition can lead to a time step which is smaller than that required by the CFL criterion. This problem can be avoided using an implicit method to treat the eddy viscosity term. In this work, the method developed by Saul’ev (1957) has been used. In this method an alternating direction sweep is employed at alternate time steps. Details can be seen, for instance, in Davies (1985b). Some boundary conditions are also required. For closed borders, a no flux condition is imposed: q=0 (18) where qis the current component which is normal to the boundary. Along open boundaries, water elevations are specified from observations and the normal component of the surface water current, q, is obtained from a radiation condition (Kowalick and Murty, 1993; Glorioso and Davies, 1995): q=c hz(19) where c=gh. Some boundary conditions are also required by the dispersion equation. There is no flux of radionuclides through a closed boundary, thus: (C (x i =0 (20) where x i is the normal direction to the boundary. Along open boundaries, the condition described in Peria´n˜ez et al. (1994) was applied: C i =aC i−1 (21) where C i is the concentration in the open boundary and C i−1 represents the concentration just inside the computational domain. The non dimensional number ais obtained from a calibration exercise. To solve the dispersion equation, a centred scheme is used for the horizontal diffusion terms, upwind differences are used for the advective terms and the Saul’ev method (Saul’ev, 1957) is again employed to solve the vertical diffusion term. It is well known that upwind differences introduce numerical diffusion. It has been shown (Prandle, 1984) that the magnitude of numerical diffusion is equivalent to increasing the diffusion coefficient K i by K% i , where: K% i =1 2(u i Dx i −u i 2 Dt) (22) where the subindex irepresents the three directions in space. Numerical diffusion in the vertical direction can be neglected due to small values of the vertical velocity (maximum value is of the order of 10 −4 m/s). In the horizontal directions, numerical diffusion has been reduced by subtracting the instantaneous value of K% h to K h , where the index h represents xor y(French, 1988). 4. Application of the model As said above, the model has been applied to study the dispersion of radionuclides in the Irish Sea. These radionuclides are discharged from a nuclear fuel reprocessing plant at Sellafield. The model has a horizontal resolution Dx=Dy=5000 m. Ten layers are used in the vertical, thus, Ds= 0.1. Time step is fixed as Dt=60 s. Stability conditions are satisfied with this selection. The computational domain is presented in Fig. 1, where the location of the nuclear fuel plant is also shown. Water depths have been introduced from bathymetric maps and range from 55 m in the west of the computational domain to a shallower area around the British coast. Water elevations are specified along the open boundary from observations (Howarth, 1990). Only the main tidal component, M 2 , has been considered. The friction coefficient has been taken as k=0.0112, a similar value to that used by Jones and Davies (1996). Wind effects have not been considered: F s =G s =0. The mean value 0.3 has been taken for ein Eq. (16) and good results are obtained selecting a=0.9 in Eq. (21). The horizontal diffusion coefficients have been fixed as K x =K y =500 m 2 /s since results in agreement with observations are obtained with these values. Indeed, Bowden (1950) suggested that for the Irish Sea 500BK x ,K y B900 m 2 /s. Fig. 1. Map of the computational domain. Letters indicate the points where tidal amplitudes and phases have been measured, numbers indicate the points where tidal currents have been measured and circles denote the points where 137 Cs concentrations have been obtained. The star is Sellafield nuclear fuel reprocessing plant. Each unit in the xand yaxis is 5000 m (grid element number). Table 1 Observed and computed tide amplitudes and phases for the points shown in Fig. 1 Point Amplitude (m) Phase (degrees) Computed ObservedComputed Observed Difference Difference 285 300 −5.0a −3.91.98 2.06 311 317b 2.95 2.92 1.0 −1.9 312 323c 2.99 3.11 −3.8 −3.4 326328 0.6d 5.93.23 3.05 6.3 329 325 1.2e 3.35 3.15 331 331f 3.10 3.08 0.6 0.0 332357 7.5g −0.042.75 2.73 339 8.0h 2.76 2.75 0.04 366 303 326I 1.90 2.30 −17.4 −7.0 324336 3.7j −4.62.51 2.63 −5.3 319 318k 0.32.48 2.62 7.8332m 2.50 2.55 −2.0 358 Differences between computed and observed values are given in % relative to the observed value. 5. Results and discussion 5.1. Water circulation Tidal amplitudes and phases, and magnitude and direction of tidal currents calculated with the model have been compared with the measured values for a number of points inside the computational domain. These points are shown in Fig. 1. Computed and observed tidal amplitudes and phases are presented in Table 1. The difference between computed and observed values in % relative to the observed value is also given. It can be seen that the difference in amplitudes is B7% for all the points and only in the case of point ‘i’a difference \10% is obtained. In the case of tidal phases, errors are B8% for all points. Observed and computed semi-major axis magnitude and orientation of the M 2 tidal current ellipse have also been compared for a number of locations and depths (Table 2). It can be seen that the model gives, in general, a good representation of current magnitude and direction in the sea. Computed and observed current profiles at points 8 and 9 (see Fig. 1) are presented in Fig. 2. The shape of the profile is reproduced by the model at both points and for both the uand 6 components of the water velocity. Thus, it seems that, in general, the model gives a good representation of the water circulation in the studied area since good agreement between observed and computed tidal amplitudes, phases and currents has been obtained. Moreover, our main objective is to study the dispersion of radionuclides and the results of the hydrodynamic part of the model seem good enough to allow an adequate description of the dispersion processes. 5.2. Radionuclide dispersion The dispersion of 137 Cs released from the nuclear fuel reprocessing plant at Sellafield has been simulated. As a first approach, l37 Cs was considered to be perfectly conservative, i.e. no fraction is removed from the water column due to biological or geochemical processes. This is usual in some models (Prandle, 1984; Abril and Garcı´a-Leo´n, 1992). Indeed, the mean value of the Cs distribution coefficient, k d , in coastal waters is 3×10 3 l/kg (IAEA, 1985). Since typical suspended matter concentrations in the eastern Irish Sea are of the order of 1 ppm (Kershaw and Young, 1988), it can be calculated that only 0.3% of the total Cs content in a given water volume is fixed to solid particles. The major source of 137 Cs to the Irish Sea has been the discharges from Sellafield. Other sources, such as nuclear weapon test fallout contributed Table 2 Observed and computed semi-major axis and orientation of the M 2 tidal current ellipse at several depths and locations shown in Fig. 1 Computed valuesObserved valuessh (m)Point Direction (degrees) Axis (m/s)Axis (m/s) Direction (degrees) 9 1.05 0.11 1.1050 0.30 0.91 10 0.932 1.550 0.44 14 0.412 52 0.92 0.62 6.6 0.52−32 −1.53 0.4920 0.60 0.40 6 0.58 8.74 45 0.98 −5 0.995 45 0.50 0.79 −7.1 0.67−11 −12.36 0.5830 0.50 −14 0.587 20 0.76 −17.90.48 −9 0.898 55 0.42 0.86 2.5 0.52−6 6.78 0.6655 0.88 0.72 3 0.919 3.145 0.62 4.60.679 45 0.80 0.66 6 Orientation is given in degrees measured anticlockwise from east (thus range from −90 to 90°). B1% of the total input (Jefferies and Steele, 1989). Observed and computed 137 Cs distributions have been compared for a number of years. The real input from Sellafield (Jefferies and Steele, 1989) was introduced in the model for each year. This input was 2970 TBq/year for 1980, which is equivalent to 5.6×10 9 Bq per time step. However, the input has been taking place since the 1960s. Thus, instead of starting the model from zero concentrations, we have assumed an uniform background of 1900 Bq/m 3 . This background represents the effect of previous discharges. In Peria´n˜ez et al. (1994) it was shown that model results do not depend upon the way the background is created. Thus, the same results would be obtained if a large discharge is performed and some time is allowed to elapse so that the discharge is distributed over the sea. To save CPU time, the uniform background option was chosen. Thus, discharges from Sellafield are carried out over this uniform background and results are obtained after a simulation period of 25 days. These results are compared with observations. Observed and computed 137 Cs concentrations in surface waters along the British coastline (Fig. 1), north and south from Sellafield, can be seen in Fig. 3. It can be seen that the model gives the general distribution pattern of 137 Cs in the sea. An intense peak is obtained at Sellafield (point 0 in the xaxis) and concentrations decrease as we move south or north from Sellafield. The activity levels measured in the sea have been reproduced by the model. A distribution map over the sea is presented as an example, in Fig. 4. This map is not significatively different from that obtained from observations (Jefferies and Steele, 1989). The shape of the 1500 Bq/m 3 isoline suggests that there is an input of non-con- taminated water between Anglesey and the Isle of Man. This is consistent with the fact that strong currents (of the order of 1.5 m/s) are observed in this area (Howarth, 1990) and these currents produce a residual flow which enters the eastern Irish Sea between Anglesey and the Isle of Man. The input from Sellafield for 1982 was 2000 TBq/year, which is equivalent to 3.8×10 9 Bq per time step. The uniform background was now selected as 1300 Bq/m 3 . As can be seen in Fig. 5, the general distribution of 137 Cs along the coast is again reproduced by the model. Fig. 2. Current profiles at points 8 (a) and 9 (b) of Fig. 1. The input for 1984 was 434 TBq/year (8.2×10 8 Bq per time step) and the background was taken as 1200 Bq/m 3 . In the case of year 1985, the input was 325 TBq/year (6.2×10 8 Bq per time step) and the background was selected as 900 Bq/m 3 . Observed and computed distributions of 137 Cs along the coast for years 1984 and 1985 can be seen, respectively, in Figs. 6 and 7. The general behaviour of 137 Cs is again reproduced by the model. Thus, it seems that the model gives a realistic representation of the dispersion processes of radionuclides in the Irish Sea. 5.3. Predicti6e studies Once the model has been tested, it can be used as a predictive tool that can be applied, for instance, in the assessment of contamination follow- Fig. 4. Distribution map of 137 Cs concentrations (Bq/m 3 )in surface water for year 1980. Each unit in the xand yaxis is 5000 m (grid element number). Fig. 3. Observed and computed 137 Cs concentrations (Bq/1) at several locations along the British coast (see Fig. 1) for year 1980 north (positive distances) and south (negative distances) from Sellafield. Fig. 5. Same as Fig. 3 but for year 1982. Fig. 7. Same as Fig. 3 but for year 1985. ing an accidental release of radionuclides in any site of the sea. It must be noted that the model has been developed for radionuclides, but can be used for other dissolved conservative pollutants by setting l=0. Also, the model can be applied to obtain some water quality parameters. As an application example, the turn-over-time of the Solway estuary (Fig. 1) has been calculated. This parameter is defined as the time in which concentration inside a bounded region decreases by a factor e −1 (Prandle, 1994). It is computed assuming an arbitrary concentration in the region of interest and obtaining the time evolution of this concentration. In our calculation, an initial concentration of 10 units/m 3 was considered in the Solway. The time evolution of this concentration is shown in Fig. 8. From this time evolution, the turn-over-time can be estimated as 2.6 days. It is interesting to note that some oscillations appear in concentrations. They are due to tidal oscillations and were already observed in a similar experiment carried out with a two dimensional model (Peria´n˜ez et al., 1996a). If an accidental discharge of a dissolved conservative pollutant occurs in the Solway, the model has estimated that concentration should decrease by a factor 0.37 in a time of :2.6 days. Two Fig. 8. Time evolution of pollutant concentration (arbitrary units) in the Solway estuary.Fig. 6. Same as Fig. 3 but for year 1984.