Using oceanography to control and forecast nuclear accidents and other passive particles problems

Using Oceanography To Control And Forecast Nuclear Accidents And Other Passive Particles Problems M. Toscano-Jimenez, J.M Abril & R.Garcia-Tenorio Engineering School. Applied Physics Department, University of Sevilla, Spain. e-mail. mtoscano d),esi. us.es Fax. +34 95 4 48 60 03 Abstract.- This article is the last improved version of a previously published model 1181 for the transport of the nuclear contamination and other passive particles in the ocean. Two interesting advances have been developed during last two years for my PhD thesis to be finished in the next months: (a) A Suspended Particulate Matter (S.P.M.) submodel, including erosion, transport and sedimentation. (b) A new advection-diffusion approach with numerical and computational improvements: Finite Elements (FE), Finite Differences (FD) and Monte Carlo (MC) methods have been compared and calibrated. These studies will be submitted to two different scientific journals to become an added guarantee of my PhD work The Baltic Sea has been elected as the validation scenario of the model and the radioisotope Cs137 is the radiotracer to be analysed. This scenario was the most contaminated ecosystem out of the Soviet Union due to the Chernobyl accident occurred at the end of April 1986, and the elected radiotracer Cs137 was the main long-lived radioisotope emitted to the environment. However, an important aim of this model is its potential usefulness in other oceanic scenarios affected by a nuclear disaster in the future. It could be an interesting tool to predict and minimize the ecological and economical impacts offuture accidents. This model can also be extended easily to non-nuclear contamination problems such as: oil accidents, nutrients dynamics and other biologicalproblems. I. INTRODUCTION It must be noted that the model is three-dimensional and its horizontal resolution is 10 km, while for the vertical resolution a total of six layers are considered. A computation time of approximately 12 hours was necessary to simulated 1 meteorological year, using a Matlab code in a personal computer (AMD- 1.4 GHz). Some approaches adopted allow to save a lot of computational time [17] and minimize numerical errors. The Circulation sub-model is based on some studies developed in the S.M.HI. (Swedish Meteorological and Hydrological Institute). The influence of the winds, tides, and inertial forces is obviously very important in the numerical simulations. Eddy-like motions with variable intensities and scales have been modelled: small (kilometres), medium (decades of km), and large (hundreds of km) Experimental information on current spectra has been analysed to be included in the circulation, diffusion and erosion-sedimentation sub-models. The model as a whole has been validated by comparing the evolution of the Cs'3 concentrations with experimental data taken from the literature. The Cs'3 model predictions are in an acceptable agreement with the experimental Cs'3 data, either in the water column or in the sediments, thereby guaranteeing the validity of the model. II. A SUSPENDED PARTICULATE MATTER (S.P.M.) SUB-MODEL An important fraction of passive particles can be fixed to the suspended matter phase present in the water column, and consequently can be deposited on the seabed. The logarithmic grain size has been used to classify the sediments: d ( = -10g2- d = diameter do = 1mm do (1) The grain size distribution (figure ?) has been modelled with a five class scheme (Bobertz et. al, 2004), and a modified Fermi function: |F(X) = IF + exp -A asktanh(x xd) F(x)~ ~~~~ Lk tX(Xd-Xne F(x): Relative weight of grains below a given size x; A=1.7 (2) { Xmd, a, sk}: {mean,variance,assimetry} parameters FEFV FUNCMN for five sediment classes in the Batic Sea .an-sze in loganthmicscae Fig. 1. Grain-size distribution (Femi Function) for five sediment classes (coarse-medium sand, fine sand, very fine sand, coarse-medium silt, medium-fine silt) in the system. A good correlation between the depth and class has been found in the Baltic Sea, and the next map (figure 2) has been calculated for this work. 1-4244-0635-8/07/$20.00 )2007 IEEE Authorized licensed use limited to: Universidad de Sevilla. Downloaded on October 21,2022 at 11:14:55 UTC from IEEE Xplore. Restrictions apply. The velocity v* is defined in relation to the bottom stress , X = p (v7)2,- and the fine particles (silt and clay) have modelled with three diameters 5={4, 7,10}. The S.PM. concentration SPM(g/m3), together with the sedimentation rate SR(glm2ly) have been modelled and validated for this article(figures 3 and 4). Fig. 2. Logarithmic fine fraction calculated for the Baltic Sea bottom. On the one hand, the resuspension events take place [11] if the velocity vf>t1,2,31cmls for the diameters O={4,7,10}. On the other hand, a mean sinking velocity Wsink=4.10 4cmls was chosen, with the next information about the erosion rate ER ([7],[14]): ER = E.f -1); E: erosion constant E=0.04gm2 Is; f: fraction of smallparticles 3 |Ce critical erosion stress : l E [0.1 ,1.5]N /m2 ; 0.4AN /m2 A published work [9] about seasonal spectra of v* have been used to calculate next tables (I and II). TABLE I: EVOLUTION OF V max(CmIs) DURING A METEOROLOGICAL YEAR. Dec- Feb - Apr - Jun - Aug - Oct - Jan Mar May Jul Sep Nov 6 4 3 2 3 4 5 3 2 1.5 3 4 2.5 1.75 1.25 1.75 2.5 3 1.5 1 1.5 2 2 1.5 1 0.5 1 1.5 TABLE II: RESUSPENSION FREQUENCY IN THE SYSTEM V*max(CMIs) 2 2.5 3 4 5 6 Events/year | 1 3 7 31 44 69 An interesting map of v*max(x,y) in the Baltic bottom [9] together with last two tables has completed the resupension dynamics of this work. The initial conditions were SPM(x,y,t=O)= 4g/m 3. The shore line, the rivers (mainly Neva, Vistula and Oder), and the organic primary production have played an important role as S.P.M. sources in accordance with the scientific literature. Fig. 3. Modelled SPM(g/m3) in the Baltic Sea for a typical January. EAt. Sea Sedimented mte(/)1 y- and 6 months oftransport. | NEVP _20 * 100 1COO 11CO Fig. 3. Modelled transport. Sedimentation Rrate (g/m2) after 1.5 years of III. NUMERICAL METHODS The Monte Carlo (MC) method has been calibrated (4) and compared with the Finite Differences (FD) method. r2-D Calibration functions: l A(x,y,t)= A Ferfr x+ h2 _ ( x h2 4 2a, (t ) 2 x (t) y + 5 ( ( ,/i: Le rf e rf Y) wih a t = [ ] (t) 2K t ; erf (z) 2 "dx -t C f 0 (4) 12M 102: 12DDC 100C 80C 600 1 2CO 1 OOCO _ FINLAND fl 300 BAL-nc REPUBLICS 0 100 200 300 400 5DO 600 700 800 900 Authorized licensed use limited to: Universidad de Sevilla. Downloaded on October 21,2022 at 11:14:55 UTC from IEEE Xplore. Restrictions apply. The MC method has been found to be better than FD method for dispersion problems with high gradients in concentrations functions, such as point-source accidents, spillages and biological blooms. It must be underlined that the MC method has been the best one to solve the numerical dispersion problem. The 2D diffusion (Kdijf25Um2 Is ) of a square (h =20km) spot is being simulated and calibrated (figures 5,6 and 7). controlled by the Prandle formula ([1],[13]) K num = (Ax - vAt)v / 2 , for non-magical conditions. Another original question of this work is the development of MC Eulerian method based on a "quasiparticles" and splines scheme (figure 8). Similar properties as the FE method have been found. A u Xi-i Xi Xi+i h 0 Fig. 8.- Scheme of MC Eulerian method, similar to FE method. IV. CURRENTS MODELLING Fig. 5: Calibration of a 2D diffusion problem using 10000 particles, Kd1f 250m2 /s, 3 months of transport. A. Circulation Sub-Model The velocity (i v) of a general point of our system can be described in a straightforward way as the sum of its annual mean velocity (7m ) plus its associated fluctuation (v ). v = v + v (5) Fig. 5: Calibration of the FD method, 1 month of diffusion. CPU- time=3.1 s. Fig. 9. Some circuits related to the annual mean currents in the Baltic Sea. Upper currents take place in the depths [0, 1 O]m, while down currents in [10, 20]m. Fig. 5: Calibration of the MC method, 1 month of diffusion. CPU- time=0.36s. The calibration of the advection has been successfully solved; the numerical dispersion has been Special attention needs to be paid to [4], which modelled the average wind-driven currents in the Baltic Sea with horizontal resolution of 10 km. This model is based on a typical meteorological year with the most probable events and their results have been validated through comparison with some experiments [5]. V PARTICLES.OJibrEd on of agrangian4knte-arloM Profi e)c-Q 3DdEVs of diffusion. Kef--250m 2/s. CPJ--Q36s rFINTE-DIFFEFENES2D Profilex=O.3Ddaysofdiff sio.Kef=25a0m2/s. CPL) 2 a E I. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on October 21,2022 at 11:14:55 UTC from IEEE Xplore. Restrictions apply. The thickness of the six layers considered are the following: [0,5]m, [5,10]m, [10,20]m, [20,40]m, [40,60]m and [60m,bottom]. The figure 1 shows the annual mean velocities for the layer [0, 5]m: some circuits have been designed in this publication, keeping in mind the vertical currents influenced by up-welling and down-welling processes studied in [12]. Measurements in the Baltic Sea ( [3] and [10] ) show that current spectra have peaks for periods similar to those of wind spectra, i.e. in the order of days. In Figure 2, and as an example, a representative current spectrum together with the different mean flows used in this work are shown. In this way, the fluctuations with regard to the annual mean velocity, v, can be expressed as the sum of two terms: V = U +W (6) where the field vi7 represents the small-scale fluctuations shorter than two days, and where wT" represents the large-scale fluctuations affecting the system longer than two days. Velocity (cm/s) East st A Annual man Daily mean flow I flow Ir + i Hourl mean flow Cit'W+ a Time (days) Fig. 10. A representative current spectrum (based on [3]) with the different mean flows proposed in this manuscript. where K, and Ky denote the horizontal diffusion coefficients, and T, and Ty are the Eulerian integral timescales in the horizontal directions. Based on equations (7) we have used a set of average experimental values of K&, Ky, Tx, and Ty compiled in [5]. The values of the average horizontaldiffusion coefficients which are compiled in Table I. TABLE III HORIZONTAL DIFFUSION COEFFICIENTS, {K ,Ky (l2 s), USED IN THE SMALL-SCALE DIFFUSION {K ,KY}m2 s 30 10 8 6 4 depth(m) [0,10] [10,20] [20,40] [40,60] [60,bott.] B. Small-scale Vertical diffusion The modelling of the vertical component of the small-scale fluctuation velocity, u,, is based on the same theoretical fundamentals as those used for the horizontal simulation. The average vertical-diffusion (table IV)coefficients are applicable to every location of the system analysed. TABLE IV VERTICAL DIFFUSION COEFFICIENTS, K, (cm2ls), USED IN THE DIFFUSION SUB-MODEL. Kz (cm2 Is) 1.00 0.50 0.20 0.20 0.05 depth(m) 5 10 20 40 60 As regards to the set of experimental values of T, and Ty, we can indicate that they were in the interval [1.5,2.5] hours, with a clear maximum at 2.0 hours. In fact, through equations (7) we can determine ux and u' at any point of the system, by the assignment of values from the Monte Carlo method, according to Gaussian distributions of probability. C. Large-scale diffusion. - A. Small-scale horizontal diffusion. - The first objective therefore is to simulate the value of the small-scale velocity fluctuation terms U' and u' x y for all the points of the ecosystem analysed. For the performance of these calculations, the variances of the two terms can be expressed in the following way: o(u 'y) u~~~ Observations and analyses carried out at several stations distributed over the Baltic Sea indicate that it can be assumed (8) for the simulation of large-scale fluctuations: w w' = WI = 0 ; U(w ) \w = Z- x ~ ~ ~ ~ W Tvx{x(8) tS(wi -F Therefore, based on (8), the calculations of the velocities {wx, w } for each point of the system were also based on the Monte Carlo method in the same way Authorized licensed use limited to: Universidad de Sevilla. Downloaded on October 21,2022 at 11:14:55 UTC from IEEE Xplore. Restrictions apply. as {vI, v } . The timescales and the variance of the velocities because of winds were in the intervals: {a(w Ix y) E- [ 3, 1 0] cm / s ;ETwx- [ 1, 31 days } (9) V. OTHER VALIDATIONS AND RESULTS The model as a whole permits the simulation of the evolution of the 137Cs concentrations in the Baltic Sea for the time interval June 1986- June 1987. 1200 1 000 800 = FINLAND SWED EN 600 H 400C BALTIC _PUBLICE no= 400 600 800 I 000 km Fig. 11.- Experimental 137Cs specific activity distribution (Bq/m3) at the surface of the Baltic Sea in summer 1986 (initial conditions), some weeks after the accident of CHERNOBYL. The experimental spatial distribution maps of the 137Cs concentrations on the initial and the final simulation dates, (see Figures 11 and 12) have been interpolated from experimental values found in several radiological journals ([6]and [15]). All the simulation was carried out using the code Matlab 7 on a personal computer AMD-1.4 GHz. The computation time needed for the simulation of the complete year (summer'86-summer'87) was approximately 12 hours. In the performed simulation, two different time steps have been chosen for the execution of the model: At, =1 days and , At, = 6 hours were adopted in agreement with the decay times of the large- and smallscale velocity fluctuations, respectively, observed in the spectra. The sedimentation sub-model was run by using Ats = 6 hours as time step. As an example, only the modelled 137Cs distribution map for the superficial layer, in June 87, is going to be shown, which needs to be compared with the experimental distribution map shown in Figure 12. Figure 13 corresponds to K 150y= 5Om2 Is (Equations 8). The satisfactory agreement is evident. Fig. 12. Experimental 137Cs specific activity distribution (Bq/m3) at the surface of the Baltic Sea in summer 1987 (final conditions). 70 14 SD 20 Fig. 13. Modelled" Cs. specific activity distributin (q In) at te200 supericia lae of th Bati Se in sume 197 wit an an 200 -b __lll ~~~~~~~~~~100 aveag vau 150' (mil winds). 100 200 300 400 500 600 700 800 900 1000)0 km Fig. 13. Modelled 137CS specific activity distribution (Bq/m3) at the superficial layer of the Baltic Sea in summer 1987, with and an average value Kw,Y=Sm /s (mild winds). VI. SOME CONCLUSIONS The Monte Carlo method has been found to be better than Finite Differences method for dispersion problems with high gradients in the concentrations functions, such as point-source accidents (oil or nuclear accidents), spillages and biological blooms. It must be underlined that the MC method has been the best one to solve the numerical dispersion problem. In the present work, the first S.P.M. model of the Baltic Sea as a whole has been presented, and also the first 3D radio-ecological model [18] of this important ecosystem. The Cs'37 model predictions are in an acceptable agreement with the experimental CS'37 data, either in the water column or in the sediments, thereby guaranteeing the validity of the model. REFERENCES [1] J.M. Abril and M.G. Leon. 1992. "A Marine Dispersion Model for Rasdionuclides and its Calibration from Non-radiological Information" J. of Environmental Radioactivity, vol.2, pp. 127-46 Fl NLAND E SWEDEN BALTIC REP UB LICS Fl NLAND SWEDEN BALTIC REP UB LICS Authorized licensed use limited to: Universidad de Sevilla. Downloaded on October 21,2022 at 11:14:55 UTC from IEEE Xplore. Restrictions apply. [2] B.Bobertz and J.Harff. 2004. "Sediment facies and hydrodynamica setting: a study of the Baltic Sea"Ocean Dynamics, vol.54, pp.39-48 [3] E. Franke. "A contribution to the investigations of the current conditions in the surface layer in the area of Darss Sill". [Proceedings 12th Baltic Oceanographic Conference, Leningrad, 1980] [4] L. Funkquist and L. Gidhagen. "A model for pollution studies in the Baltic Sea". SMHI-Report RHO-39, 1984 [5] L. Gidhagen, L. Funkquist and R. Murthy, "Calculations of horizontal exchange coefficients series current meter data from the Baltic Sea". SMHI-Report RO- 1, 1986. [6] Z.G. Gritchenko et al., "Radiation situation in the Baltic Sea in 1986 in sea water and sediments". Three years observations of the levels of some radionuclides in the Baltic Sea after the Chernobyl accident. Baltic Sea environment proceedings, vol.31 pp.10-30, Helsinki, 1989(a). [7] J.T. Holt and I.D.James. 1999. Continental Shelf Research, vol. 19, pp. 1617-1642. [8] D.J. Jansen et al. 2003. "Deposition of organic matter and paticulate nitrogen and phosphorus at the North Sea -Baltic Sea transition- a GIS study", Oceanologia, vol.45, No.1, pp. 283-303. [9] A. Jonsson, A. Danielsson and L. Rahm. 2005. "Bottom type distribution based on wave friction velocity in the Baltic Sea". Continental Shelf Research, vol. 25, pp. 419-435 [10] ,J. Kielmann, W.Krauss and K.H. Kennecke, "Currents and stratification in the Belt Sea and the southern Arkona Basin during 1962-68". Kieler Meeresforsh., vol. 29, 1973. [11] C. Kuhrts, W. Fennel & T. Seifert. 2004. "Model studies of transport of sedimentary material in the western Baltic". Journal of Marine Systems, vol. 52, pp.167-190 [12] A.Lehmann and H.H. Hinrichsen. 2000. "On the wind driven and thermohaline circulation of the Baltic Sea", Physics and Chemistry of the Earth(B), vol.25, No.2, pp.183-189. [13] D.Prandle. 1984. Phil. Trans. R. Soc. Lond A., vol. 310, pp. 407-36 [14] W.Puls and J.Sundmann. 1990 "Simulation of suspended sediment dispersion in the North Sea".In: Residual currents and long term transport. pp.356- 372 [15] A. Stigebrandt, A. 1983. "A model for the exchange of water and salt between the Baltic Skagerrak", Journal of Physical Oceanography, vol.13,pp.411-27, 1983. [16] A. Svansson, "The water exchange of the Baltic", Ambio Special Report, vol.1, pp.15-19, 1972. [17] M. Toscano-Jimenez and R. Garcia-Tenorio, "Modelling The Dispersion Of 137Cs In Marine Ecosystems With Monte Carlo Methods", Nuclear Instruments And Methods B, vol.213, pp.789-793, Application To The Baltic Sea After The Chernobyl Disaster", Ocean Engineering, vol.31, pp.999-1018, 2004. [19] A. Viopio, "The Baltic Sea". Elsevier, Amsterdam, pp.162-167, 1981. [20] D. Weiss, "The distribution of radionuclides in bottom sediments of the Open Baltic Sea and Greifswald Bodden", Three years observations of levels of some radionuclides in the Baltic Sea after the Chernobyl accident. Baltic Sea environment proceedings, vol.31 pp.94-122, Helsinki, 1989. 2003. [18] M. Toscano-Jimenez and R. Garcia-Tenorio, "A Three-Dimensional Model For The Dispersion Of Radioactive Substances In Marine Ecosystems. Authorized licensed use limited to: Universidad de Sevilla. Downloaded on October 21,2022 at 11:14:55 UTC from IEEE Xplore. Restrictions apply.