Simulating Long Range Transport of Radioactive Aerosols Using a Global Aerosol Transport Model

Study of long range transport of radioactive gases and aerosols is necessary to estimate radiological impact to the members of public and environment during nuclear reactor accidents. In the present work, an attempt has been made to employ an atmospheric circulation model that predicts meteorological parameters at the regional and global scales, and a transport model that utilizes these meteorological parameters for the radioactive aerosol dispersion. Non-hydrostatic Icosahedral Atmospheric Model (NICAM) coupled with Spectral Radiation-Transport Model for Aerosol Species (SPRINTARS) is adapted to fulfil this objective, which is used for simulating effects of conventional aerosols on atmospheric pollution and climate system. As an illustrative case study, global simulation is carried out for a horizontal resolution of ~110 km to model the dispersion of radioactive aerosol (S, I and Cs) releases from the Fukushima Daiichi Nuclear Power Plant (FDNPP) accident. The results obtained from this case study are compared with the available literature data and other simulation results. The statistical analyses show that the comparisons are better for locations far away (> 150 km) from the emission location, and the results are further discussed. Continuous run of this system will help in predicting the activity concentration in forecast mode, and it may be used for decision support, particularly in the case of long range transport (> 100 km).


INTRODUCTION
Very low probable accidents in radiological facilities such as Fukushima Daiichi Nuclear Power Plant (FDNPP) may lead to the release of radioactive materials in the environment.The radionuclides in the form of gases and aerosols are transported in the atmosphere from local to global scale based on the energy and quantity of release from these facilities.The study of long range transport of these radioactive materials assist in estimating their radiological consequences and nuclear emergency management.Atmospheric concentration of these radionuclides and their ground deposition constitute the important parameters required to quantify the radiological impact.Near-and farfield simulation of these parameters from the emission/release point are carried out using air pollution and regional/global models respectively.Various dispersion models based on Lagrangian and Eulerian framework coupled with meteorological models have been used widely to simulate the long range transport.Apart from the study of radiological consequences, these simulations will also help in estimating aerosol residence time and transport pathways in the atmosphere (Papastefanou, 2010;Huh et al., 2013).
Lagrangian models such as Flexible Particle dispersion model (FLEXPART) and Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) coupled with meteorological models were used to simulate the dispersion of airborne radionuclides in the environment (Srinivas and Venkatesan, 2005;Srinivas et al., 2012;Rakesh et al., 2015).The numerical accuracy of Lagrangian based models depends on the number of air parcels or particles released and the resolution of meteorological input data.Eulerian models are also used frequently to simulate regional and global radioactive aerosol transport (e.g., Takemura et al., 2011;Danielache et al., 2012;Christoudias et al., 2013;Draxler et al., 2015).Draxler et al. (2015) conducted a regional scale simulation using five different atmospheric transport and deposition models and compared their ensemble mean results ( 137 Cs and 131 I air concentration and their deposition) with the measurements at a single location about 110 km from the FDNPP.Similarly, Kristiansen et al. (2016) performed a global scale study to estimate aerosol lifetimes using 137 Cs and 133 Xe as radiotracers released from Fukushima accident using an ensemble of 19 global models.This study shows that the ratio of modelled to observed surface aerosol concentration varies many orders of magnitude depending on the model and location.
In the present work, a similar attempt has been made to implement a non-hydrostatic atmospheric circulation model (NICAM -Nonhydrostatic ICosahedral Atmospheric Model) coupled with an aerosol transport model (SPRINTARS -Spectral Radiation-Transport Model for Aerosol Species) to simulate the global atmospheric dispersion of radioactive aerosols.NICAM-SPRINTARS was previously used for simulating effects of aerosols on atmospheric pollution and climate system at global as well as regional scale (Dai et al., 2014;Goto, 2014;Goto et al., 2015).The important features and advantages of this model for aerosol simulations are highlighted by Goto et al. (2015).
In the present study, an emission module for radioactive aerosols is introduced in NICAM-SPRINTARS model by considering their particle size distribution, release rate and correction for radioactive decay during their transport.The model is then applied to simulate transport of three specific radionuclides ( 35 S, 131 I and 137 Cs) attached to aerosols, released during FDNPP accident.Furthermore, simulated atmospheric radionuclide concentrations are compared with the measurements and other model data from the literature, and statistical analyses are carried out.

Description of the Model
NICAM is a Global Cloud Resolving Model (GCRM) developed by Tomita and Satoh (2004), having fully compressible non-hydrostatic dynamical core and icosahedral horizontal grids.It is a seamless global-to-regional model which does not involve nesting technique and boundary conditions from global model (Goto et al., 2015).More details of the code development and its capabilities are found in various literature (Satoh, 2002;Tomita et al., 2002;Tomita et al., 2005;Satoh et al., 2008).Tomita (2008) developed a new grid transformation method which forms stretched grid that is isotropic and inhomogeneous.In the present study, NICAM uses Large Scale Condensation (LSC) for the cloud microphysics scheme described by Le Treut and Li (1991) with cumulus parameterization by Chikira and Sugiyama (2010) and Chikira (2010).To model the subgrid scale turbulence of the planetary boundary layer, Mellor-Yamada Nakanishi-Ninno Scheme, level 2 (MYNN2) is used (Satoh et al., 2014).Takemura et al. (2000) using the framework of atmospheric general circulation model, MIROC (Model for Interdisciplinary Research on Climate).However, it can be coupled to any other general circulation model, like NICAM, to provide the meteorological parameters.SPRINTARS is a global aerosol transport model which simulates the dispersion of conventional aerosols like soil dust, sea salt, black carbon, organic carbon, sulfate aerosols in the atmosphere.Takemura et al. (2002) simulated long range transport of Asian dust and anthropogenic aerosols from East Asia, and studied aerosol optical properties in the global scale using SPRINTARS.Furthermore, SPRINTARS is used to study changes in the meteorological field due to direct and indirect effects of atmospheric aerosols (Takemura et al., 2005;Suzuki et al., 2013).This aerosol transport model mainly includes the process of emission, advection, diffusion and deposition.The deposition process includes dry deposition, wet deposition and gravitational settling (Takemura et al., 2000).The dry deposition flux is due to turbulent mixing and the wet deposition flux includes subcloud scavenging (due to gravitational coagulation) and incloud scavenging (removal of aerosol in cloud water by precipitation).In this model, particle size distribution is split into 10 different size bins equally spaced in log scale, ranging from 0.1 to 10 µm.Aerosol micro-physical processes such as nucleation, coagulation,and condensation are not considered in this model.

Code Modifications
Since the present study involves simulation of radioactive aerosol dispersal from a single release location, the emission module of the code needs to be modified.The introduction of point source emission for radioactive aerosols in the code system requires the information about grid structures and their formation.The emission location for a given latitude and longitude is fixed using three parameters, viz., processes handling regions, region division number (r-level, it is a key parameter for deciding how many regions are made to distribute all the grid points), and grid point locations in the code (Sarkar, 2015).In order to account for the radioactive decay of the radionuclide, the aerosol emission module is modified by introducing the radioactive decay term in the governing equation of SPRINTARS code.The modified governing equation for radioactive aerosol transport is then given by, where, ρ air is the density of air, q a is the aerosol mixing ratio, defined as the ratio of mass of radiaoctive aerosols to the mass of dry air in a given volume., w are horizontal and vertical wind components respectively, F E , F D , and F S are the emission, diffusion, and deposition fluxes of aerosols, respectively, and λ is the radioactive decay constant.

Input Parameters
Source term denotes information about the actual or potential releases of radioactive material from a given source which includes total quantity, release rate, chemical composition, physical form, mode of release, and energy accompanying release.Most of the radiologically important radionuclides get released to the atmosphere either in the form of gas or aerosols.The transport behaviour of radioactive aerosol particles is governed by many parameters such as height of release, energy of release, source emission, radioactive half-life of each radionuclide, particle size distribution, etc.In the present work, we study the transport of aerosols bearing three radionuclides ( 35 S, 131 I and 137 Cs) released from the FDNPP accident. 35S is an activation product generated from neutron activation of 35 Cl and 34 S present in the cooling sea water which is injected into the damaged reactor.The neutrons for the activation of 35 Cl and 34 S are produced from the probable re-criticality of the melted fuel during the early phase of the FDNPP accident.The sea water cooling was supplied up to 12 days from March 14 till the fresh water supply was restored.Hence, 35 S emission in the numerical simulations is carried for ~12 days starting from 14 th March.However, 131 I and 137 Cs are the fission products produced in the reactor core, their releases from the damaged fuel started much earlier (11 th March) compared to 35 S attached aerosol emission, and they are significant upto 30 days (11 th April).Since the formation and release mechanisms are different for 35 S and fission product aerosols ( 131 I and 137 Cs), the total release period and starting time of release are different in their respective simulations.

Sulphur-35
35 S is a radionuclide with half-life of 87 days naturally produced in the atmosphere by cosmic ray spallation of 40 Ar present in air (1%).High levels of 35 S activity at different sampling sites around the world was observed during the progress of FDNPP accident.It was then attributed to the 35 S production by slow neutron capture of 35 Cl and 34 S present in the coolant seawater (Priyadarshi et al., 2011).
In the present study, we have used an estimated average emission flux of ~10 16 atoms s -1 of 35 S in two particle size categories of effective radii, 0.13 µm and 0.2 µm with equal weightage (Danielache et al., 2012).The release period started from 12 UTC, March 14, 2011 and continued for ~12 days (Danielache et al., 2012).

Iodine-131
131 I is one of the major fission product with a fission yield of 2.83 percent and it is easily absorbed in the thyroid gland causing short-term exposure to the public in the case of accidental releases.It is a significant contributor to the internal dose from beta (maximum 606 keV with 89% abundance) and gamma (364 keV with 81% abundance) emissions.Approximately 1.5 × 10 17 Bq of total 131 I (both vapor and aerosol forms) was released to the atmosphere during FDNPP accident (Meszaros et al., 2016).We have considered 131 I aerosol emission rates (Fig. 1) given by Katata et al. (2015) in the present simulation with a characteristic size distribution as shown in Table 1.The aerosol size distribution is lognormal with 0.4 µm as mode diameter and Geometric Standard Deviation (GSD) of 1.3 (Baklanov and Sorensen, 2001).

Model Setup/Initialization
NICAM is initialized with NCEP Final Analysis (FNL) data (http://rda.ucar.edu/datasets/ds083.2) having a spatial resolution of 1° × 1° available at a temporal resolution of 6 hours.In the present NICAM simulation, the horizontal resolution of 110 km corresponds to g-level 6 (globally), and 40 vertical levels are considered.The upper layer of the first level from the Earth surface is approximately located at 80 m and the top most vertical layer is situated at ~40 km height from mean sea level.Horizontal grids are distributed across 160 regions (r-level 2), and these regions are assigned to 16 processes for parallel computation in the present study.The simulation is performed with nudging of meteorological parameters (horizontal wind components -u & v and temperature -T) extracted from FNL data for all vertical layers, and the necessary output variables are averaged over one hour.Global NICAM-SPRINTARS model simulations are carried out in the ANUPAM-AGGRA parallel computing system (Nuclear India, 2014).The emission location is fixed at a grid point (37.4°N,141.0°E) as this is the closest to the FDNPP.The simulation period is 12 days (288 h) for 35 S dispersion since the neutron emission and 35 S production rates are negligible beyond 12 days as shown by Danielache et al. (2012).However, significant amount of 131 I and 137 Cs emissions continued upto 30 days and hence the simulations are carried out for one month with a step size of 240 s.The release rates of these radionulides have been determined using inverse estimation methods by Priyadarshi et al. (2011) for 35 S and Katata et al. (2015) for 131 I and 137 Cs.The aerosols are released from the lowermost vertical layer of model which has upper bounday at ~80 m elevation.Chemical reactions during the transport are not considered and hence they behave as passive tracers.

Global Dispersion of 35 S Simulations
The dispersion of 35 S attached aerosols is simulated with constant emission flux of ~10  (Giulianelli and Willard, 1974;Priyadarshi et al., 2011) near the emission location.Further transport of the sulfate aerosols in the atmosphere is considered as a passive tracer throughout the simulation period.The results are then compared with the literature data to examine the model performance.We have used number concentration (atoms m -3 ) instead of activity concentration (Bq m -3 ) for comparing the results with available literature data directly.Fig. 2 shows 35 S attached aerosol concentration (1 h average values at the surface) after 3 and 11 days of transport.The results show that most of the activity is transported towards the Pacific Ocean during the initial release period (11 days) except for a few days when the plume is directed towards East and moving over land surface (Fig. 2).The simulation results show that the spread of activity is mainly within the northern hemisphere, and the 35 S concentration varies between 10 3 -10 7 atoms m -3 .
Present simulation results are also compared with i) measurements and ii) model predictions by Danielache et al. (2012) at two particular locations viz., Atmosphere and Ocean Research Institute (AORI), Japan [35.9°N, 140.0°E] and Scripps Institute of Oceanography (SIO), USA [32.8°N, 117.3°W]where 35 S air concentration was continuously measured and reported during the FDNPP accident (Priyadarshi et al., 2011).Present simulation shows that the temporal variation of 35 S concentration is similar to that of observations and are able to identify the peak values for both the locations.Comparison of some of the peak values at these locations is as follows: at AORI, 35 S concentration reached a maximum value of 1.7 × 10 4 atoms m -3 on 24-25 March, 2011, present simulation and Danielache et al. (2012) predict ~10 6 and ~10 7 atoms m -3 respectively during this period and identified the peaking behaviour; at SIO, measurements show that the concentration reaches a peak value of 1050 atoms m -3 during 24-28 March, 2011.In this case, models (present simulation and Danielache et al., 2012) predict a peak concentration of ~10 3 atoms m -3 and they are in close agreement with the observations.These results show that the model predictions are better for faraway loactions from the source emission, and the estimates for the near-source location is poor due to the global grid  (Anand andMayya, 2011, 2015), aerosol particle properties such as wet scavenging rates, size distribution, and meteorological parameters.

I and 137 Cs Releases into the Atmosphere and their Global Dispersion
131 I and 137 Cs aerosol releases from FDNPP accident are studied in detail by Katata et al. (2015) compared to 35 S release scenario and their time-dependent release rates are shown in Fig. 1.Present simulations are carried out for 30 days with these release rates and input parameters described in Section 2. Figs. 3 and 4 show 1-h average surface concentration (Bq m -3 ) of 137 Cs and 131 I respectively from the NICAM-SPRINTARS simulations.Simulation results show that the wind flow pattern is in the South-East direction over the Pacific Ocean initially, and after two days (March 13, 18 UTC) the wind direction is towards East.As a result, the radioactive plume with trace amounts of radioactivity reaches North America and European continent after 7 and 10 days respectively.This global dispersion pattern is in close agreement with the published work of Meszaros et al. (2016).Comparison of Figs. 3 and 4 show that short half-life (also shorter residence time) of iodine isotope restricts the area of dispersion and concentration in the atmosphere within 30-days of simulation, although the radioactive inventory of 131 I released from FDNPP is approximately one order higher than 137 Cs (Fig. 1).Our model results are also compared with the online version of HYSPLIT model (the results are not shown here) and found to be in good agreement.

Comparison of Model Results with Measurements
NICAM-SPRINTARS simulation results are compared with some of the radioactivity measurements ( 131 I and 137 Cs air concentrations) that were carried out at different locations around the globe during the progress and after FDNPP accident.Furthermore, statistical analysis are carried out between the measurements and simulation datasets which will indicate the model performance.The following parameters are estimated in this analysis: i) Normalised Mean Square Error (NMSE), ii) Fractional Bias (FB) and iii) Correlation Coefficient (CC).Correlation coefficient is one of the most robust parameter (Stohl et al., 1998), however multiple statistical parameters are used to arrive at the better model performance (Draxler, 2006).
Measurement and model data of 131 I and 137 Cs airborne concentrations at 140.5°E], about 110 km away from FDNPP (Draxler et al., 2015;Lebel et al., 2016, for   The correlation coefficient between measurement dataset and NICAM-SPRINTARS model data for the Tokai location is estimated as 0.56 for 131 I and -0.1 for 137 Cs; FB is -1.63 and -1.59, and NMSE is 101.1 and 123.88 for 131 I and 137 Cs respectively.These values are compared with that of Draxler et al. (2015) for many model performances as shown in Table 2.
The differences in the model comparisons (CC for 137 Cs, FB and NMSE for 131 I) are due to i) Cs and I source terms used in the models; and ii) horizontal resolution in the present study (~100 km).The source term in the present study is from Katata et al. (2015) whereas Draxler et al. (2015) considered the source term from Terada et al. (2012).Although the horizontal resolution of ~100 km is considered as finer for global level simulation, it is coarser for the regional model and hence the deviation in the simulation results at this location (~100 km to the source emission point).The other reasons for the disagreement with the measurement data are due to the conversion of Iodine aerosols to vapours and vice versa during transport which could have changed the deposition rate of aerosols and hence the air concentration.The comparisons with other models by Draxler et al. (2015) shows that NICAM-SPRINTARS is acceptable for simulating radioactive aerosols, although some biases are found.
Since, Tokai-mura is nearer to the emission (FDNPP) in the global scale, we compare some of the simulation data with measurements at far locations such as Ukraine (Chaisan et al., 2013) and other CTBT locations published in the literature (CTBTO, 2011).The temporal variation of 131 I and 137 Cs air concentrations for a particular site (Ukraine, 51.0°N, 30.0°E) is shown in Fig. 5 along with model results.It shows that the present simulation captures the arrival of radioactive plume at this location after ~10 days, and the results are in agreement with the measurements (Chaisan et al., 2013)     model performance at 5 out of 9 locations, RN70, RN75, RN53, RN63 and RN61; however few sites (RN38, RN34) are poorly correlated (Table 3).Furthermore, a scatter plot analysis of 131 I and 137 Cs air concentrations from NICAM-SPRINTARS simulations and measurement data is carried out.Fig. 6 shows nine CTBTO locations and Ukraine data are pooled together in this plot.Overall, it is noted that this scatter plot are within the factor of 10 (as the band indicated in Fig. 6) for 131 I, and similarly 55% for 137 Cs.All these statistical parameters (CC, NMSE and FB) and scatter plot analysis indicate good predictability of the present model.However, some of the differences between present model results and measurements/other simulations can be attributed to input parameters such as source term and particle characteristics (size), and aerosol microphysical processes expended in the numerical model.

CONCLUSIONS
The present study demonstrates the applicability of NICAM-SPRINTARS for simulating long range transport of radioactive aerosols.Comparison with other models and measurements shows higher CC, and smaller NMSE and FB values that denotes good performance of the NICAM-SPRINTARS model, although some biases are found in few cases.This study also shows that source term, aerosol microphysical processes, and particle size characteristics are the main factors that affect the radioactive aerosol concentration in the atmosphere.Hence, inclusion of aerosol microphysical processes and appropriate input parameters in the model will enhance the estimation of aerosol metrics.The present NICAM-SPRINTARS system in forecast mode will help in predicting the activity concentration, and hence it may be used for nuclear emergency management, particularly for long range, beyond 100 km.Also, it will be used to evaluate the validity of emission scenarios.The existing system will be further strengthened by the addition of all relevant radionuclides and their properties associated with aerosols.

Fig. 2 .
Fig. 2. Simulated global atmospheric concentration (atoms m -3 ) of 35 S after (a) 3-d, (b) 11-d 131 I only) are compared with the present simulation results.The individual concentration values from NICAM-SPRINTARS are closer to 10-ensemble model mean fromDraxler et al. (2015) at many temporal data points.

Fig. 6 .
Fig. 6.Scatter plot of 131 I and 137 Cs air concentration simulated by NICAM-SPRINTARS model vs. measurement data.

Table 1 .
Size distribution of 131 I and 137 Cs aerosols.

Table 3 .
Similarly for 131 I, CC varies from -0.51 to 0.49, NMSE from 1.83 to 20.27, and FB from -1.17 to 1.18.These estimated statistical parameters show better

Table 2 .
Comparison of statistical parameters at Tokai, Japan.

Table 3 .
Statistical parameters -comparison ofmodel and measurement data.