A particle-tracking method for simulating the dispersion of non-conservative radionuclides in coastal waters

a A particle-trackingmethod for simulating the dispersion of non-conservative radionuclides in coastal waters R. Peri!aa*nnez a, *, A.J. Elliott b Departamento Fı´sica Aplicada I, University of Sevilla, EU Ingenierı´aTe ´cnica Agrı´cola, Ctra. Utrera km 1, 41013 Sevilla, Spain b Centre for Applied Oceanography, University of Wales, Bangor, Marine Science Laboratories, Menai Bridge, Anglesey LL59 5EY, UK Abstract A particle-tracking method has been used to simulate the dispersion of non-conservative radionuclides in the sea. Three dimensional turbulent diffusion and the interactions between water, suspended matter and bottom sediments are simulated using a stochastic method. Kinetic transfer coefficients, as in finite difference models, are used to describe the transfers between the liquid and solid phases. Deposition of suspended matter and erosion of sediment are also included in the model. The method has been applied to simulate the dispersion of 137 Cs and 239,240 Pu in the English Channel and the results have been compared with those of a finite difference model. The results from both techniques are, in general, in good agreement. Keywords: Model; Dispersion; Particle-tracking; Finite difference; Radionuclides; English Channel 1. Introduction The state-of-the-art models used to simulate the dispersion of radioactivity in the sea consist of finite elements or, mainly, finite difference models that solve the hydrodynamic equations together with the advection–diffusion dispersion equation, in either a two dimensional or a three dimensional form (Harms, 1997; Peri! aa* nnez & *Corresponding author. Tel.: +34-95-4233669; fax: +34-95-4232644. E-mail addresses: [email protected] (R. Peri! aa* nnez), [email protected] (A.J. Elliott) Reguera, 1999). If the model is applied to simulate the dispersion of nonconservative radionuclides, then the interactions with the solid phases (suspended matter and bottom sediments) must also be considered. This implies that the suspended matter equation, including the settling of particles, deposition and erosion of the sediment, must be solved too (Aldridge, 1998; Margvelashvily, Maderich, & Zheleznyak, 1997; Piasecki, 1998; Peri! aa* nnez, 1999, 2000). Absorption/desorption reactions are usually described using kinetic transfer coefficients instead of the less appropriate equilibrium distribution coefficients, kd. These models are computationally expensive due to the fact that a large number of equations must be solved. Also, since finite difference models work with radionuclide concentrations, the whole computational grid must be swept each time step. Finally, the CFL condition (Kowalick & Murty, 1993) limits the size of the time step that can be used if an explicit scheme is used to solve the hydrodynamic equations. In consequence, large CPU times may be required even using recent supercomputers (Peri! aa* nnez, 1999). An alternative approach is to make use of particle-tracking techniques. In this kind of model, the impact of the CFL criterion can be reduced by making the hydrodynamic computations off-line. Also, advection can be simulated to a high degree of accuracy since numerical dispersion (which appears when the advection equation is solved with finite differences) is not introduced. Particle-tracking models have already been used to simulate the dispersion of conservative tracers and oils spills (Hunter, 1987; Elliott, Dale, & Proctor, 1992; Proctor, Elliott, & Flather, 1994a, b). In such applications, a release of radioactivity to the sea is modelled as a number of discrete particles, each particle being equivalent to a number of units (Bq, moles, atoms, etc.). Then the path followed by each individual particle is computed, turbulent diffusion being modelled as a three dimensional random walk (Monte Carlo) process. The density of particles is calculated to obtain the radioactivity concentrations at the end of the simulation. The main difficulty that appears in the simulation of the dispersion of non-conservative radionuclides is the treatment of absorption and desorption: how to decide if each particle is fixed to suspended matter or bottom sediments (if initially dissolved) or if it is redissolved (if initially present in the suspended matter or the bottom sediment). The main contribution of this paper consists of a new method developed to solve this problem: a formulation (suitable for a particle-tracking model) to describe the transfers of radionuclides between water, suspended matter and bottom sediments, based upon kinetic transfer coefficients and a stochastic method, is presented. It is important to point out that exactly the same physical parameters as in the equivalent finite difference models are used. The particle-tracking modelling technique is well suited to problems in which high contaminant gradients are involved, since numerical diffusion is not introduced. This, together with the fact that it can give very fast answers, even in a PC, allows the technique to be considered as a very useful predictive tool in the assessment of contamination following accidental or deliberate releases of radionuclides. The particle-tracking model is presented in the next section, then an application to the English Channel is shown. The model has been used to simulate the dispersion of an instantaneous hypothetical release of radionuclides from La Hague nuclear fuel reprocessing plant, located on the French shore of the English Channel. Simulations have been carried out for two radionuclides with very different geochemical behaviours: the relatively conservative 137 Cs and the high reactive 239,240 Pu. A finite difference model of the Channel has been previously developed and validated through the study of the dispersion of these radionuclides, comparing observed and computed concentrations in water, suspended matter and bottom sediments (Peri! aa* nnez and Reguera 1999; Peri! aa* nnez, 2000). Thus results of the particle-tracking model will be compared with the output of the finite difference model in equivalent simulations and the relative advantages of each technique will be assessed. 2. The particle-tracking model 2.1. Advective transport The position vector of a given particle, rtþDtðÞ, at time tþDtis computed from rtþDtðÞrtðÞ Dt¼qtðÞ;ð1Þ where Dtis the time step used in the model and qthe current vector of components u and valong the xand yaxes, respectively. Currents are obtained by running a hydrodynamic model in advance. Standard tidal analysis is used to determine the tidal constants (tide amplitude and phase) for each grid cell of the hydrodynamic model. These constants are evaluated for both components of the flow and can be derived for as many tidal constituents as desired. In this work, only the two main semidiurnal tides, M2and S2, will be considered. Once the tidal constants are known, computation of the flow vector, q, just involves the calculation and addition of a few cosine terms. As a consequence, the evaluation of the tidal advective transport of particles is very fast and is not limited by the CFL criterion. For real applications, first order accuracy in the particle-tracking scheme is adequate. In simulations of the movement of drogues in an estuarine environment, Elliott and Clarke (1998) found no improvement in the results when a second order accuracy scheme was used to simulate the movement of surface drifters by the particle-tracking technique. Moreover, in ocean dispersion problems, the effects of turbulence will mask any small errors in the advection scheme. The net residual current in the modelled area must be added to qsince a residual transport cannot be generated with the pure harmonic tidal currents that are used in particle-tracking calculations. The residual flow vectors can also be obtained from the (previously run) hydrodynamic model. Wind-induced transport can be included in the model by assuming that the surface wind-induced current is a percentage of the wind speed, generally 2–3% (Proctor et al., 1994b). This current decreases logarithmically below a depth z1(the thickness of the wind-driven surface layer) to zero at a depth z2. Eckman theory (Pugh, 1987) predicts that the surface windinduced current due to a steady wind blowing over deep water should be deflected to the right of the wind direction (in the northern hemisphere). However, observational evidence suggests that this deflection angle can be neglected (Proctor et al., 1994b) in shallow coastal waters. 2.2. Turbulent diffusion Three dimensional diffusion is simulated using a random walk method. It has been shown (Proctor et al., 1994b; Hunter, 1987) that it is a simulator of Fickian diffusion provided that the maximum size of the horizontal step given by the particle, Dh,is Dh¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 12KhDt pð2Þ in the direction 2pRAN, where RAN is a random number between 0 and 1. This equation gives the maximum size of the step. In practice, it is multiplied by RAN to obtain the real size at a given time and for a given particle. Similarly, the size of the vertical step is Dv¼ffiffiffiffiffiffiffiffiffiffiffiffiffi 2KvDt pð3Þ given either towards the sea surface or the sea bottom. Khand Kvare the horizontal and vertical diffusion coefficients, respectively. 2.3. Radioactive decay Consider the radioactive decay equation: qC qt¼lC;ð4Þ where lis the radioactive decay constant. This equation can be treated using a stochastic method if it is assumed that the probability pof removal of a particle at each time step is (Hunter, 1987; Proctor et al., 1994b) p¼1elDt:ð5Þ In practice, a random number is generated for each particle on each time step. If RAN4pthen the particle is removed from the computation. 2.4. Transfers between water, suspended matter and bottom sediments Consider a two phase system. If the transfers of radionuclides between the two phases are described through the kinetic transfer coefficients k1and k2, the equations that give the time evolution of activity in the two phases are qA1 qt¼k1A1þk2A2;ð6Þ qA2 qt¼k1A1k2A2: These equations are easily solved using finite differences. In particle tracking, a label is given to each particle to differentiate if it is in phase l or 2. If the particle is in phase 1, the probability p1that the particle goes to phase 2 in each time step is p1¼1ek1Dt:ð7Þ Similarly, if the particle is in phase 2, the probability p2that it goes to phase 1 each time step is p2¼1ek2Dt:ð8Þ Thus, in particle tracking, the exchanges between two phases can be modelled as two decay processes with probabilities p1and p2. These processes are treated as the radioactive decay process described above. If a given particle goes from one phase to the other, its label is changed and the new corresponding decay process is considered at the next time step. This method has been compared with the finite difference solution of the system of Eq. (6). It has been considered that all radionuclides are, at t¼0, in phase 1. Thus, the solution given by each method refers to the percentages of radionuclides that are in phase 1 at each following time step. Results obtained by both methods are presented in Fig. 1, using 200 and 10,000 particles in the stochastic simulation. High fluctuations occur with 200 particles, but the finite difference solution is well modelled if 10,000 particles are used in the particle-tracking calculation. The mean value and standard deviation of the difference between finite difference and stochastic solutions are presented in Table 1 for different numbers of particles in Fig. 1. Time evolution of the percentage of radionuclides in phase 1 given by the solution of the differential equations using finite differences and the particle-tracking method, with 200 and 10,000 particles. the stochastic method. It can be seen that both the mean value and standard deviation of the difference decrease as the number of particles considered in the stochastic simulation increases. Acceptable results are obtained for numbers of particles of the order of 10,000. The stochastic method can be extended to the case in which there are three different phases: water, suspended matter and active bottom sediments (particles with a diameter 562.4 mm). It is considered that the transfer of radionuclides from the solid phases to water is governed by a kinetic transfer coefficient k2, the transfer from water to suspended matter by k1m and the transfer from water to the sediment by k1s. The absorption of radionuclides depends on the surface of particles per water volume unit. Thus the exchange surface and exchange velocity concepts have been used (Peri! aa* nnez, Abril, & Garc! ııa-Le! oon, 1996; Peri! aa* nnez & Mart! ıınez-Aguirre, 1997; Peri! aa* nnez, 1999, 2000). Following these papers: k1m ¼w1 3m rR;ð9Þ k1s ¼w1 3Lf f RH ;ð10Þ where w1is the exchange velocity, mis the suspended matter concentration, rand R are the density and mean radius of suspended matter particles, Lis the average mixing depth (the distance to which the dissolved phase penetrates the sediment), f is a correction factor that takes into account that not all the mass of the sediment is in contact with water and Hgives the thickness of the water layer above the sea bottom that interacts with the sediment. In a two dimensional depth-averaged model, His equal to the water depth. Since particle tracking is three dimensional, H is left as a free parameter to be calibrated. The decay equations that are equivalent to the differential equations that describe transfers between the three phases, presented for instance in Peri! aa* nnez (2000), are qCd qt¼k1mCdk1sCd;ð11Þ Table 1 Mean value and standard deviation of the difference between the percentage of radionuclides in phase 1 given by finite differences and the stochastic method a NP hDis 100 0.654 2.861 200 0.637 2.220 500 0.296 1.557 1000 0.335 1.021 10,000 1.56 10 2 0.309 100,000 6.17 10 2 0.114 a NP is the number of particles used in the simulation. qCs qt¼k2Cs;ð12Þ qAs qt¼k2fAs:ð13Þ where Cd,Csand Asare radionuclide concentrations in water, suspended matter and active bottom sediments, respectively. A label is given to each particle to classify in which phase it is present. Depending on the label of the particle, the corresponding decay equation is treated. If the particle is in suspended matter, the probability that it goes to the dissolved phase, in each time step, is p¼1ek2Dt:ð14Þ Similarly, if the particle is in the bottom sediment, the probability that it is redissolved is p¼1ek2fDt:ð15Þ If the particle is initially dissolved and its distance to the sea bottom is smaller than H, it can go to any of the two solid phases with a probability p¼1ek1mþk1s ðÞDt:ð16Þ A random number is generated to decide if the particle is effectively removed from solution. If it is, the normalized probability that the particle goes to the sediment is calculated as p¼ps pmþps ;ð17Þ where pm¼1ek1mDt;ð18Þ ps¼1ek1sDt:ð19Þ A second random number is then generated. If RAN5p, then the particle goes to the sediment. If RAN >p, then it goes to the suspended matter. Of course, if the distance of the particle to the sea bottom is larger than H, only the decay to suspended matter is considered since such particles cannot interact with the sediment. A numerical experiment has been carried out to test the method in which a volume of water with a given suspended matter concentration and active sediment on the bottom is considered. A dissolved radioactive tracer is added and the equations that give the time evolution of activities in the three phases are solved using finite differences and the stochastic method. The following realistic parameters have been used: w1¼2:1108m/s, k2¼1:2105s 1 ,m¼10 ppm, R¼15 mm, r¼2600 kg/m 3 ,L¼0:01 m, f¼1, f¼0:1, H¼0:2 m with 10,000 particles being used in the particle-tracking simulation. The comparison between both methods is presented in Fig. 2, where the time evolution of the fraction of tracer that is dissolved, in suspended matter and in the sediment is presented. The simulation shows that the stochastic method solution is in very good agreement with the finite difference solution for the three phases. Indeed, solutions corresponding to both methods cannot be distinguished in the case of water and sediment. 2.5. Suspended matter deposition and sediment erosion The suspended matter concentrations over the model domain can be obtained by running in advance a finite difference suspended matter model. Results are then analyzed in a similar way to currents, so that the suspended matter concentration at each point and for any time can be obtained as the simple calculation of cosine functions. Suspended matter falls to the sea bottom with a settling velocity ws. If a particle (in the particle-tracking sense, not a suspended matter particle) is fixed to suspended matter, its position above the bottom, z, at time tþDtis obtained from ztþDtðÞztðÞ Dt¼ws:ð20Þ If ztþDtðÞ40, then the particle is considered to fix to the sediment and its label is appropriately changed. A standard formula for flocculation has been used to represent the increase in the settling velocity as the suspended matter concentration increases (Clarke & Elliott, 1998; Eisma, 1993; Peri! aa* nnez, 2000) ws¼a1ma2;ð21Þ where a1and a2are obtained from measurements or from model calibration. The particle-tracking model is three dimensional. If the finite difference suspended matter model is depth-averaged, its output is the depth-averaged suspended matter concentration. A Rouse profile is then used to resolve the vertical structure of suspended matter. This allows the calculation of the suspended matter concentration at height zabove the bottom, mz, from the depth-averaged suspended matter concentration m(Clarke & Elliott, 1998): mz¼m1ws=bku*  h z  ws=bku *;ð22Þ where his water depth, kis the von Karman constant (0.4), bis an arbitrary constant usually taken as 1 (Eisma, 1993; Clarke & Elliott, 1998) and u*is the scalar friction velocity u*¼kq jj ln h=z0  1;ð23Þ where z0is the bottom roughness. The corresponding value of mzis used to calculate k1m at the position of each particle from Eq. (9). Erosion of the sediment has been described in terms of the erosion constant concept (Nicholson & O’Connor, 1986). Thus, the probability that a particle is removed from the sediment and incorporated to the water column as suspended matter is p¼1eEf jqjMDt;ð24Þ where Eis the erosion constant and Mis some power of the water velocity typically in the range 2–5 (Prandle, 1997). Thus, the erosion description used in previous finite difference models (Peri! aa* nnez, 1999, 2000) has been converted into a stochastic form. It is considered that erosion can only take place if the water velocity is larger than a Fig. 2. Time evolution of the fraction of radionuclides in water, suspended matter and bottom sediments given by finite differences and by the particle-tracking method with 10,000 particles. Solid lines correspond to the stochastic solution and dashed lines to the finite difference solution. 3.2. Discussion Particle tracking is a powerful tool that can be applied in the assessment of radioactive contamination following an accidental release in aquatic environments in general. Also, the method can be applied to both conservative and non-conservative radionuclides, using the same conceptual approach for the interactions between Fig. 6 (continued). Fig. 7. 239,240 Pu activity concentrations in water (mBq/l), suspended matter (Bq/g) and active sediment (Bq/g) given by the finite difference (a) and the particle-tracking (b) models. liquid and solid phases as used in finite difference models. The particle-tracking technique, combined with tidal analysis and mean flow databases for the hydrodynamics, presents two clear advantages over finite differences: speed of Fig. 7 (continued). computation and the fact that it does not introduce numerical diffusion. The number of particles used in simulations depends on the speed of computation and/or the accuracy required. For example, Proctor et al. (1994b) used 4000 particles to simulate an oil spill during real-time forecasting of an incident in the Arabian Gulf. However, it is possible that, due to the geochemical behaviour of certain radionuclides, the activity concentrations in some of the phases may not be resolved. This is the case for 137 Cs in suspended matter. The number of particles used in the simulation should be of the order of 10 6 to be able to calculate, with reasonable accuracy, activity concentrations in suspended matter but then the computing time would be similar to that of the finite difference model. However, if we are interested in the assessment of contamination following an accident, it is probably enough to calculate activities in water and bottom sediments, which are the phases where almost all the radioactivity is present. If the most important aspect is the speed of computation, the number of particles can be reduced. For instance, if 137 Cs dispersion is simulated using 20,000 particles, good results are still obtained in water and bottom sediments and the computation time is reduced to 2.6% of that required by the finite difference model. In the case of 239,240 Pu, the number of particles could be reduced more and activities in the three phases may still be calculated since Pu is more widely distributed between the three phases than Cs. The particle-tracking method is fully 3-D and takes account of horizontal advection plus turbulent diffusion in the x–y–zdirections. In the present application, the tidal and residual flows have been computed by a depth-averaged 2-D hydrodynamic model. A 2-D hydrodynamical model is adequate due to the very strong tidal currents and the vertically well-mixed character of the coastal waters. However, the particle-tracking method is suitable for use with a 3-D hydrodynamic model (e.g. Harms, Karcher, & Dethleff, 2000) and the manner in which vertical stability can inhibit vertical mixing can be parameterised in such applications by making the vertical step-size a function of the stability of the water column. 4. Conclusions The purpose of this paper is to demonstrate the viability of new particle-tracking algorithms for the simulation of processes that involve chemical speciation. The particle-tracking method is commonly used for oil and chemical spill applications where there is a need for rapid response simulations. In such applications, a significant increase in computational speed is achieved by computing the hydrodynamics ‘off-line’ and then using tidal prediction and mean flow databases to reconstruct the water movement. This avoids the CFL criterion during the rapid response application, although the particle-tracking time-step must be kept reasonably short in order to maintain the accuracy of the first order advection scheme that is used to compute the particle movement. While it would also be feasible to solve finite difference speciation equations using a pre-computed flow field, the proposed technique has the advantage that the computational time will be very short if a relatively small number of particles is used in the simulation. Moreover, the method will be more efficient if the contaminant patch covers only a small fraction of the grid domain (Hunter, 1987) and the accuracy of the results can be improved by releasing larger numbers of particles}although at the expense of an increased simulation time. One of the main advantages of the particle-tracking method is the ad-hoc manner in which particle characteristics can be defined. For example, non-Fickian dispersion can be readily simulated by attaching an age to each particle (equal to the time since the release of the particle into the computation) and then making the diffusive step of the particle a function of its diffusion time. In the new model, the interactions between the liquid and solid (suspended matter and bottom sediment) phases have been described in terms of kinetic transfer coefficients, as in actual finite difference models. A stochastic method has been developed to simulate such interactions. Deposition of particles and erosion of the sediment are also included in the model. A Rouse profile is also used to estimate the vertical structure of suspended matter concentrations. As a demonstration of the method, the particle-tracking model has been used to simulate the dispersion of radionuclides in the English Channel. A finite difference model was previously developed and validated for the Channel. Thus the output from both models has been compared. Two radionuclides with different geochemical behaviours, 137 Cs and 239,240 Pu, have been used for the comparisons. The agreement between both models is, in general, rather good. However, it is possible that, due to the particular behaviours of certain radionuclides, activity concentrations in one of the phases cannot be calculated. Acknowledgements R Peri! aa* nnez was supported by ENRESA and the EU under contract FIGECT2000-00085. A J Elliott was funded by the Chief Scientist’s Group of MAFF under contract AE1021 and by the Procurement Executive of the MoD under contract NPC1B/2012. References Aldridge, J. N. (1998). CSERAM: A model for prediction of marine radionuclide transport in both particulate and dissolved phases. Radiation Protection Dosimetry,75, 99–103. Clarke, S., & Elliott, A. J. (1998). Modelling suspended sediment concentrations in the Firth of Forth. Estuarine, Coastal and Shelf Science,47, 235–250. Eisma, D. (1993). Suspended matter in the aquatic environment. Berlin: Springer. Elliott, A. J., & Clarke, S. (1998). Shallow water tides in the Firth of Forth. Hydrographic Journal,87, 19– 24. Elliott, A. J., Dale, A. C., & Proctor, R. (1992). Modelling the movement of pollutants in the UK shelf seas. Marine Pollution Bulletin,24, 614–619. Guegueniat, P., Herrmann, J., Kershaw, P. J., Bailly du Bois, P., & Baron, Y. (1996). Artificial radioactivity in the English Channel and the North Sea. In P. Guegueniat, P. Germain, & H. Metivier (Eds.), Radionuclides in the oceans (pp. 121–154). Les Ulis: Les Editions de Physique. Harms, I. H. (1997). Modelling the dispersion of 137 Cs and 239 Pu released from dumped waste in the Kara Sea. Journal of Marine Systems,13, 1–19. Harms, I. H., Karcher, M. J., & Dethleff, M. (2000). Modelling Siberian river runoff}implications for contaminant transport in the Arctic Ocean. Journal of Marine Systems,27, 95–115. Hunter, J. R. (1987). The application of Lagrangian particle tracking techniques to modelling of dispersion in the sea. In J. Noye (Ed.), Numerical modelling. Applications to marine systems (pp. 257– 269). North-Holland: Elsevier. Kowalick, Z., & Murty, T. S. (1993). Numerical modellingof ocean dynamics. Singapore: World Scientific. Margvelashvily, N., Maderich, V., & Zheleznyak, M. (1997). Threetox: A computer code to simulate three dimensional dispersion of radionuclides in stratified water bodies. Radiation Protection Dosimetry,73, 177–180. Nicholson, J., & O’Connor, B. A. (1986). Cohesive sediment transport model. Journal of Hydraulic Engineering,112, 621–640. Peri! aa* nnez, R. (1999). 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, 37–51. Peri! aa* nnez, R. (2000). Modelling the tidal dispersion of 137 Cs and 239,240 Pu in the English Channel. Journal of Environmental Radioactivity,49, 259–277. Peri! aa* nnez, R., Abril, J. M., & Garc! ııa-Le! oon, M. (1996). Modelling the dispersion of non conservative radionuclides in tidal waters. Part 1: Conceptual and mathematical model. Journal of Environmental Radioactivity,31, 127–141. Peri! aa* nnez, R., & Mart! ıınez-Aguirre, A. (1997). U and Th concentrations in an estuary affected by phosphate fertilizer processing: Experimental results and a modelling study. Journal of Environmental Radioactivity,35, 281–304. Peri! aa* nnez, R., & Reguera, J. (1999). A numerical model to simulate the tidal dispersion of radionuclides in the English Channel. Journal of Environmental Radioactivity,43, 51–64. Piasecki, M. (1998). Transport of radionuclides incorporating cohesive/non cohesive sediments. Journal of Marine Environmental Engineering,4, 331–365. Prandle, D. (1997). Tidal characteristics of suspended sediment concentrations. Journal of Hydraulic Engineering,123, 341–350. Proctor, R., Flather, R. A., & Elliott, A. J. (1994a). Forecast and hindcast simulations of the Braer oil spill. Marine Pollution Bulletin,28, 219–229. Proctor, R., Flather, R. A., & Elliott, A. J. (1994b). Modelling tides and surface drift in the Arabian Gulf: application to the Gulf oil spill. Continental Shelf Research,14, 531–545. Pugh, D. T. (1987). Tides, surges and mean sea level. Wiley: Chichester (p. 15).