CC BY
49
Probabilistic Dispersion Models In Large Water Areas. Statistics And Stochastic Computational Algorythms
O. Sorokovikova, D. Dzama, D. Asfandiyarov, D. Blagodatskikh
The Nuclear Safety Institute (IBRAE) olga_sorokov@mail .ru,_dzama. dv@mail .ru, dasfandiyarov@ibrae.ac.ru, blagodat@ibrae.ac.ru
Abstract
A Monte Carlo method for calculation of dispersion in large water areas with complex coastlines is presented. Having utilized a multi-year database of sea currents and of distributions of the mixed layer depth, a model of statistical ranking of water areas based on contamination levels is proposed. Techniques for parallelization of the Monte-Carlo method and statistical postprocessing of results have been realized. The method and the techniques have been validated against releases into a particular large water area.
Keywords: Monte Carlo method, ocean dispersion model, pseudorandom number generator, parallel calculations
I. Introduction
Now a significant number of objects related to the Cold War nuclear legacy are present in Russian coastal water areas [1-2]. Such objects are concentrated in the Arctic and Far-East regions. The inventory of these objects comprises coastal radioactive waste repositories (Andreeva Bay, Gremikha Bay, Sayda Bay in the Arctic region and Razboynik Bay in the Far-East region) as well as about 18 000 radiation hazardous objects resting on the seafloor of seas of the Arctic ocean including containers with solid radioactive waste and sunken nuclear submarines.
Therefore, the Arctic and the Far-East regions are considered to pose potential hazards. This is mainly due to nuclear submarines which are kept afloat near coastal storage facilities. Though, the probability of occurrence of pertinent emergency situations is very low a radiation safety analysis of possible consequences of emergency situations for the environment was conducted. It was demonstrated that under unfavorable conditions high concentrations of contamination can be maintained in plume [3].
According to the International Atomic Energy Agency (IAEA) recommendations state-of-the-art software packages are to be used in optimization of monitoring strategies near radiation hazardous objects through discerning dominant contamination pathways in the environment [4].
In the Nuclear Safety Institute, a software package allowing calculation of probabilistic characteristics of possible levels of contamination and tracing of key water areas for establishing monitoring sites has been developing. The exceedance of concentration in those water areas may be considered as a trigger for extending of the routine monitoring strategy.
The necessity of a probabilistic approach to modelling of consequences of radioactive releases into coastal waters is justified by the seasonal variability of ocean currents as well as uncertainty in the source parameters.
Ensemble approaches based on multy-year meteorological data have been widely used for modelling of dispersion of contaminants in the atmosphere and allow conservative estimates to be obtained for various localizations of the source. Finnish software package SILAM [5] can be considered as an example of such an approach.
As regards ocean currents, similar approaches are still at early stage of development due to the fact that for large water areas pertinent databases have been compiled relatively recently. Four-dimensional data in the ocean can only be obtained through reanalysis of the coupled ocean -atmosphere circulation based on data assimilation. Furthermore, information on ocean currents is sparse even in the vicinity of coastlines. Yet, in recent years the situation has changed for good due to ever increasing amount of corresponding data. Sophisticated ocean models allowing free access to reanalysis data have been developing by scientific groups in various countries.
Therefore, in order to tackle such problems, a real-time mathematical model and input data for a long period of observations with a high space- and time- resolution are needed. Furthermore, a probabilistic analysis is a computationally demanding task that cannot be performed without introducing an efficient parallelization method. Thus, an appropriate database of ocean currents for a particular water area, efficient usage of computing resources and tools for postprocessing are at the core of practical implementation of any probabilistic approach for estimation of contamination in large water areas.
The probabilistic model presented in this article is intended for performing analysis of ocean currents and for determining various contamination pathways in large water areas in the Arctic and Far-East regions of Russia. As a test case a hypothetical release into Razboynik Bay in the Far-East region has been simulated. The choice of that particular water area is not arbitrary due to the presence of one of the largest coastal radioactive waste repositories.
I. General Assumptions
In this article a quasi three-dimensional ocean dispersion model for calculation of surface contamination is presented. It is assumed that in the upper layer of the ocean contaminants drift with the current. Contamination is regarded either as a solution or fine-disperse aerosol affected by gravitational deposition.
The main assumption of the model is that a homogeneous distribution in the mixed layer (above the seasonal thermocline) is present. Thus, it is supposed that turbulent mixing is negligible below the seasonal thermocline whereas above it intense vertical turbulent mixing occurs due to wind forcing. Vertical currents into deeper layers of the ocean are attributed to large-scale horizontal divergence of surface currents.
The mathematical model employs the transport equation for Lagrangian particles and one which determines the increase in particles dimension. In addition the model utilizes a parametrization of the turbulent diffusion coefficient, takes into account large-scale downwelling and interaction of particles with a coastline.
II. Governing equations
The governing equation of the model is the convection-diffusion transport equation in the spherical coordinates:
II. Mathematical formulation of the model
dHC dHCU dHCV 5
-+-+-=-
dt Rcosddfi Rdd R2 cos
+
where the following notation is used: C - concentration (Bq/l), $ and d - longitude and latitude, R - Earth radius (m), H - the mixed layer depth (m), U and V - the zonal and meridional components of wind velocity (m/s), / - the turbulent diffusion coefficint (m2/s), <?(t) - the time -dependent source rate (Bq-m2/s). H, C, U, V , / depend on $, d and t. The source can also occupy a finite volume.
As a computational algorithm of the proposed model, a Monte Carlo method utilizing Lagrangian particles is used. At any moment contamination is regarded as a set of Lagrangian particles. Each particle is characterized by the longitude and the latitude of the center of the particle, the horizontal dimension and the activity of the particle. The distribution of concentration within each particle is described by the two-dimensional normal distribution with the dispersion equal to the horizontal dimension of the particle. The vertical dimension of the particle depends on the mixed layer depth. The distribution of concentration in the Eulerian coordinates is calculated as the total contribution from the particles.
The equation (1) is solved applying the Lagrangian coordinates. Hence, the probability density function of Lagrangian particles is a solution for the equation if the motion of the particles is governed by the Ito stochastic differential equation:
dX = Udt + B-dW (2)
where dX is the change in the position a particle over the time interval dt, U - the velocity of the particle, B - diagonal matrix with non-zero elements determining the intensity of turbulence, W -the Wiener process.
The advection and diffusion members of the equation (2) are numerically discretized as follows:
* Ut
1
n+—
i 2 =$n +-
1
n+-2
n-
2
Rcosd"
Vt
+ —
R (3)
Rcosd"
1
R
where t - the time integration step; ¡3 is a variable model parameter ranging from 0 to 1, a - the horizontal dimension of the particle. When ¡ = 1 a classical Lagrangian stochastic model is realized, i.e. turbulent diffusion affects only the displacement of particles, whereas when ¡ = 0 turbulent diffusion has no influence on the positions of the centers of particles. In this model the coefficient 3 is equal to 0,9.
The increase of the horizontal dimension of particles is given by the following equation:
(a2 p =(a2 )"+2(l -¡)jUt (4)
As parametrization schemes of the turbulent diffusion coefficient ^ two models are employed - those of Ozmidov [6] and Smagorinsky [7]. In the Ozmidov model the turbulent diffusion coefficient is ascribed to the whole plume taking into account its dimensions. In the Smagorinsky model this coefficient is determined through calculation of the local value of the shear stress tensor thereby being different for various computational cells. At any moment the maximum of the two magnitudes is ascribed to every particle.
The distribution of concentration in the Eulerian coordinates is integrated as contribution from the clouds of particles. The horizontal dimension of the cloud of a particle is assumed to be finite and is equal to three radii of the particle ("the three-sigma rule").
The application of Gauss particles with finite dimensions, as opposed to a classical Lagrangian stochastic model, significantly loosens restriction imposed on the total number of particles in order to obtain "smoothed" solutions. The choice of the numerical value of the parameter p be equal to 0,9 is justified in the article [8].
The integration time step is limited by the maximum value of the flow velocity since the change in the displacement a particle must not exceed the dimensions of an Eulerian cell:
U t< L, where L is the horizontal dimension of the Eulerian cell.
max 2
The effect of large-scale downwelling at the lower border of the mixed layer is introduced in the model. The corresponding vertical current W conveys contamination to the lower cold layers. The decrease of the total activity in a fluid column of the height H and of the surface S is given by the following equation:
dA
— = -pWS (5)
dt
where p = A/ (S- h) is the specific volume activity.
Thus, over the period of the duration t concentration is diminished by a value given by
f Wt~~\
AA = A
1 -e
III. Interaction with the coastline
The way the distribution of concentration is calculated in the vicinity of the coastline is different as opposed to the other parts of computational domain. Hence, all cells are split into computational cells which are occupied by fluid and ghost cells where no motion of fluid occur. Thus, the coastline is represented by a number of computational cells that adjacent to ghost cells. The horizontal dimension of a particle in the vicinity of the coastline is adjusted so as to guarantee intersection of the particle only with computational cells. At a distance from the coastline, which is equal to 3a , the particle dimension is strictly determined by turbulent diffusion. If 3a is larger than a distance from the coastline, then 3a is limited by that distance. Thus, the increase of concentration within a particle approaching to the coastline is modelled as a consequence of the law of conservation of mass.
In Fig. 1 results obtained using a classical Lagrangian approach and a puff Lagrangian stochastic approach are presented. In the latter case one need to filter out artificial values of concentration occurring in the ghost cells, while in the former one particles do not intercept ghost cells due to the fact that their radii are equal to zero. In order to obtain a smooth distribution of concentration in the case of a classical Lagrangian approach a significant increase in number of particles is required. The usage of Gauss particles allows for a smooth distribution of concentration to be obtained while putting a reasonable upper limit on the total number of particles.
Trajectories of Lagrangian particles are calculated regardless of computational mesh. Hence, the motion of particles in the vicinity of the coastline has to be corrected in order to avoid intersection of the center of a particle and the coastline. If the displacement vector of a particle occurring over the time step intersects the coastline, then the method of images is used. Thus, the part of the segment of the displacement vector lying outside the coastline is reflected. Having made a few such iterations, one can force the center of the particle to be in a computational cell. Then, the minimal distance from the center of the particle to the coastline is determined with 3a of the particle being equal to that distance.
1 000e+00 1 347e.00
1 814e<JM
2 443e+00 3290e»00 4 431e+00 5967e+00 8 036e+00 1 082e*01 1 458e+01 1 9S3e»01 2.644e+01 3561e+01 4 795e+01 6458e+01 8 698e*01 1.171 e+02 1 578e»02 2125e+02 2861 e»02
3 854e+02
5190e+02 6980e.02 9.413e+02
1 268e.03 1.707e+03
2 299e*03 3097e.03 4171C+03 5.617e+03 7 565e+03 1 0195.04 1 372e*04
1 848e.04
2 489e»04 3352e+04 4.514e*04 6079e.04 8187e.04 1 103e.05 1 485e.05 2000e*05
Figure 1: Illustration of different approaches to modelling of interaction between the plume of contamination and the
coastline
IV. Input data
In order to perform a calculation, the following parameters are required:
— a two-dimensional array where each element stands either for land or water;
— a relief map for postprocessing;
— a multy-year database of the zonal and meridional components of velocity in the upper layer in the Arctic and the far-East-regions;
— the mixed layer depth distribution in a given region.
The two-dimensional array used to represent the coastline in the model is interpolated from Earth relief data which is quite accessible in the present [9-12]. Databases of ocean currents are generated as a result of simulations of various coupled ocean-atmosphere general circulation models. In the last years, such data have been provided for free access by the largest research centers [13-16]. In this article only free access databases were used.
V. Parallelization method
The accuracy of modelling strongly depends on the total number of particles. In order to obtain statistically significant results, the number of particles is requiered to be large enough. Hence, the efficiency of a method is determined by the parallelization technique applied which in terms of hardware is the way how a whole task is divided into processes. Trajectories of the particles being independent, the solution of the problem can be obtained through performing N independent calculations. The whole release is divided equally between N processes which are carried out independently from each other. Thus, summing results obtained in each process one can obtain the solution of the whole task.
In order to avoid occurrence of identical trajectories in different processes, a unique statistic is ascribed to each process via utilizing different subsequences of a whole sequence of pseudorandom numbers with a long period. In this article all results are obtained using the pseudorandom generator (PNG) [17].
That PNG constitutes a multiplicative congruent sequence with a very long period:
where the modulus M = 2128, the multiplier A = 5100109 modM and the period P = 2126. The values thus obtained conform to the uniform distribution, which in turn can be transformed into values described by the standard normal distribution using the Box-Muller transform. A whole sequence
(6)
ak =uk/M.
of the length P generated by the PNG can be arbitrary divided into a number of subsequences of equal lengths nested into the whole sequence. In order to discern such a subsequence, one needs to know the first element of each subsequence, which is determined for the m-th subsequence of the length n as follows:
u0 =1
u , = u A(n) mod M,
m+1 m W ' (7)
a =u /M, ( )
m m ' '
A(n) = An modM.
Therefore, one does not need to recalculate the values generated before given subsequences in order to skip those subsequences, which is very important in parallel calculations.
Two levels of nesting are used: the first one is a subsequence applied for calculation of a given scenario of emergency situation, the second one constitutes a set of subsequences for running different processes.
The peculiarity of the PNG used in this article lies in its program realization due to the fact that it operates with very large numbers which cannot be processed without overflow using even 64-bit computing. In the model a version of the PNG [18] written in the C language is used.
III. Simulation results
As an example of hypothetical emergency situation various scenarios of releases into the ocean near Razboynik Bay (one the largest repositories of radioactive waste) were considered. The source of releases is located at 132,36°E, 42,89°N. To determine a typical pattern of propagation of contamination in that water area a number of calculation with varying parameters of the release were conducted.
In the first calculation a release of 90Sr and 137Cs was simulated with the total activity of the source as follows: 90Sr -1,87 • 1015 Bq, 137Cs - 6,22 • 1016 Bq. The source rate was kept constant during the calculation which covered a period from 1 January 2001 to 31 December 2004.
In the figures 2 - 4 instant fields of concentration (Bq/l) are presented for some dates within the aforementioned period. The legend used in these figures coincides with that applied for the figure 5.
Knowing a series of values of concentration in every point of the water area obtained for different dates, one can determine a field of statistical characteristics. In figure 5 maximum levels of concentration with a 95% confidence interval are plotted.
This probabilistic estimate is most applicable for the case of a long-term release in the environment. In case of emergency situation when the duration of the release is much less than the characteristic time scale of propagation of contamination, ensemble calculations are used. It is assumed that the release might occur at an arbitrary moment within a given time interval. For every such a moment corresponding simulation of propagation of contamination with fixed parameters of the source is performed. Thus, the empirical distribution function in a given point of the water area is obtained. Given a value of the confidence interval, seldom peaks of concentration can be cut off.
Figure 2: Instant fields of concentration on 10 February 2001 (left) and 2002 (right)
Figure 3 Instant fields of concentration on 20 July 2001 (left) and 2002 (right)
Figure 4: Instant fields of concentration on 25 November 2001 (left) and 2002 (right)
Figure 5: Maximum levels of instant concentration with a 95% confidence interval
The aforementioned approach was tested against the releases in the same water area (Razboynik bay) with the same localization of the source whereas the discharge was supposed to occur instantly. The date of the release varied from 1 to 31 January 2008. Each simulation covered a period of three months so as to allow contamination to be removed from the calculation domain.
A series of time-integrate concentration in every point of the water area is obtained. It represents the impact of the release in an arbitrary point of the water area. In the same way as described before a probabilistic field with a given confidence interval is determined. In the figure 6 maximum levels of time-integrated concentration (Bq*s/l) in every point of the water area with a 95% confidence interval are plotted.
Figure 6: Maximum levels of time-integrated concentration with a 95% confidence interval
IV. Conclusion
In this article a Lagrangian stochastic model for calculation of ocean dispersion and probabilistic approaches to estimations of simulation results are presented. It was demonstrated in the case of a hypothetical emergency situation at Razboynik Bay that stable distribution patterns of contamination may occur in a large water area. Therefore, it is confirmed that an analysis utilizing the methods described in this article can be used for optimization of sampling strategy near coastal nuclear legacy repositories under assumption of the seasonal variability of ocean currents.
Funding
This work was financially supported under a contract between the State Atomic Energy
Corporation ROSATOM and the Nuclear Safety Institute within the framework of the R&D project "A practical methodology for integral safety analysis of nuclear legacy sites, radioactive waste storages, including the development and implementation of a system of computational tools and forecasting software packages. Stage 2017-2019".
References
[1] Problems of the nuclear heritage and ways to solve them. Development of the radioactive waste management system in Russia. - Ed. by Bolshov L.A., Laverov N.P., Linge I.I. - 2013. - 392 p. - V 2.
[2] The Cold War's legacy at the bottom of the Arctic. Radioecological, technical and economic problems of radiation rehabilitation of the seas / Sarkisov A.A., Sivintsev V.L., Vysotsky V.L., Nikitin V.S.; Nuclear Safety Institute (IBRAE) RAS. - Moscow, 2015. - 699 p.: ill. - ISBN 978-59907220-0-2 (bound).
[3] Semenov V.N., Sorokovikov A.V., Sorokovikova O.S. Conservative models for assessing the pollution of surface waters in the event of extreme hypothetical accidents and utilization of nuclear submarines in the Kamchatka Peninsula. Proceedings of IBRAE RAS / Ed. by L.A. Bolshov; Nuclear Safety Institute (IBRAE) RAS. - Moscow: Nauka, 2007 - Issue 9: Modeling of Radionuclide Transport in the Environment [in Russian] / Ed. by R.V. Arutyunyan. - 2008. - 229 p.: ill. - ISBN 978-5-02-036954-2 (bound).
[4] Safety report series N19. Generic models for use in assessing the impact of discharges of radioactive substances to the environment. IAEA, Vienna, 2001.
[5] Sofiev M., Siljamo, P., Valkama, I., Ilvonen, M., Kukkonen, J., A dispersion modelling system SILAM and its evaluation against ETEX data. Atmosph. Environ., Vol. 40, pp. 674-685, 2006, DOI:10.1016/j.atmosenv.2005.09.069.
[6] Monin A.S., Ozmidov R.V., Turbulence in the Ocean, Environmental Fluid Mechanics Vol. 3, 1985, Springer Science & Business Media, 248 p.
[7] Smagorinsky, J., General Circulation Experiments with the Primitive Equations. Monthly Weather Review, 1963. 91(3): p. 99-164.
[8] Belikov V.V., Goloviznin V.M., Katyshkov Yu.V., Semenov V.N., Starodubtseva L.P., Sorokovikova O.S., Fokine A.L. NOSTRADAMUS - computer system for predicting the radiation situation. Verification of the model of contamination transport in the atmosphere. Proceedings of IBRAE RAS / Ed. by L.A. Bolshov; Nuclear Safety Institute (IBRAE) RAS. - Moscow: Nauka, 2007 -Issue 9: Modeling of Radionuclide Transport in the Environment [in Russian] / Ed. by R.V. Arutyunyan. - 2008. - 229 p.: ill. - ISBN 978-5-02-036954-2 (bound).
[9] NOAA. ETOPO1 Global Relief Model [Web resource] https://www.ngdc.noaa.gov/mgg/global/global.html.
[10] Walter H. F. Smith and David T. Sandwell. Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings. Science 26 Sep 1997. - Vol. 277, - Issue 5334, - P. 1956-1962. - DOI: 10.1126/science.277.5334.1956.
[11] Farr T. G. The Shuttle Radar Topography Mission / T.G. Farr [et al.]// Rev. Geophys. -2007. - Vol. 45. - RG2004, doi:10.1029/2005RG000183.
[12] ASTER global digital elevation map announcement [Web resource] https://asterweb.jpl.nasa.gov/gdem.asp.
[13] Geophysical Fluid Dynamics Laboratory [Web resource] https://www.gfdl.noaa.gov/ocean-data-assimilation/.
[14] Climate Data Guide. ORAS4: ECMWF OCEAN REANALYSIS AND DERIVED OCEAN HEAT CONTENT. [Web resource] https://climatedataguide.ucar.edu/climate-data/oras4-ecmwf-ocean-reanalysis-and-derived-ocean-heat-content.
[15] APDRC of the IPRC. The FRA-JCOPE2 17 years NW Pacific Ocean reanalysis data. [Web
resource] http://apdrc.soest.hawaii.edu/datadoc/fra-jcope2.php.
[16] APDRC of the IPRC. NRL NLOM 1/32° Nowcast. [Web resource] http://apdrc.soest.hawaii.edu/datadoc/nlom_30day.php.
[17] Mikhailov G.A., Marchenko M.A. Parallel realization of statistical simulation and random number generators. // Rus. J. Num. Anal. Math. Model- 2002. - Vol. 17. - No. 1. pp.113-124.
[18] Institute of Computational Mathematics and Mathematical Geophysics RAS, Laboratory of Monte Carlo Methods. [Web resource] http://osmf.sscc.ru/~mam/index_rus.html
