Spatiotemporal variation of radionuclide dispersion from nuclear power plant accidents using FLEXPART mini-ensemble modeling

We investigate the spatiotemporal distribution of the radionuclides including iodine-131 (¹³¹I) and cesium-137 (¹³⁷Cs), transported to Qatar from fictitious accidents at the upwind Barakah nuclear power plant (B-NPP) in the United Arab Emirates (UAE). To model the dispersion of radionuclides, we use the Lagrangian particle–air parcel dispersion model FLEXible PARTicle (FLEXPART) and FLEXPART coupled with the Weather Research and Forecasting model (FLEXPART–WRF). A four-member mini-ensemble of meteorological inputs is used to investigate the impact of meteorological inputs on the radionuclide dispersion modeling. The mini-ensemble includes one forecast dataset (Global Forecast System, GFS) and three (re)analysis datasets (native-resolution and downscaled NCEP final analysis – FNL, as well as downscaled ERA5). Additionally, we explore the sensitivity of the radionuclide dispersion simulations to variations in the turbulence schemes, as well as the temporal and vertical emission profiles, and the location of emission sources. According to the simulated age spectrum of the Lagrangian particles, radionuclides enter southern Qatar about 20 to 30 h after release. Most of the radionuclide deposition in the study area occurs within 80 h after release. The most populated areas of Qatar coincide with moderate ¹³¹I concentrations and ¹³⁷Cs deposition, while uninhabited areas in southern Qatar receive the highest amounts. A larger number of long-lived particles is found in the FNL-based simulations, which is interpreted as a greater dispersion of particles at a greater distance from the emission location. The highest simulated ¹³¹I and ¹³⁷Cs deposition shows a pronounced spatiotemporal pattern. The largest impacts are found in the south and southeast of Qatar, during the early daytime development of the boundary layer, and during the cold period of the year. The results show remarkable differences in the spatiotemporal distribution of ¹³¹I and ¹³⁷Cs simulations based on the FNL and GFS datasets, which share a common base meteorological model. As part of a sensitivity analysis involving different model setups, changing the emission point from B-NPP to Bushehr NPP (Bu-NPP) results in a reduced transfer of radioactive materials to Qatar, except in the spring season. Bu-NPP simulations reveal distinct spatial patterns, with peak ¹³¹I concentrations and ¹³⁷Cs deposition observed in northern and eastern Qatar during winter and spring.

1 Introduction

Modeling the spatiotemporal distribution of radioactive compounds using chemical transport models (CTMs), especially after large nuclear accidents, has received widespread attention (Chino et al., 2011; Stohl et al., 2012; Christoudias and Lelieveld, 2013; Evangeliou et al., 2017). Whether explicitly stated or not, these studies aim to determine the extent and transport of radionuclides at various spatial scales. They also aim to improve the performance of CTMs, which are at the core of preparedness programs for potential nuclear accidents (or releases) associated with the growing number of nuclear facilities (Farid et al., 2017). However, case studies of large accidents, which have a typical duration of a few days to a few weeks, are not suitable for studying the effects of diurnal and seasonal (atmospheric) changes on radionuclide dispersion and the sensitivity to the meteorological (re)analyses driving the simulations.

Maurer et al. (2018) found that the performance of dispersion models strongly depends on the successful modeling of boundary layer processes, which vary with season and time of day. Meteorological inputs to the CTMs are generated by atmospheric models, which may exhibit variability when simulating atmospheric conditions over nuclear accident areas (Arnold et al., 2015). Long et al. (2019) studied the effects of the East Asian northeastern monsoon on the transport of radionuclides from the Fukushima nuclear power plant accident to the tropical western Pacific and Southeast Asia. They found that in these regions, radioactivity levels are lower than in other regions of the Northern Hemisphere, which is due to the late arrival of the radionuclide plumes carried by the monsoon circulations: that is, the dispersion of radionuclides from this accident could potentially be different under other atmospheric conditions, which are only captured by the hypothetical, iterative simulation of this event at different times of the day and year.

While the licensing of new nuclear power plants is rare in the Western world and many older Soviet-era stations are approaching decommissioning, the Middle East–North Africa (MENA) region has seen a rise in the planning, proposal, and recent operation of nuclear facilities. The Barakah nuclear power plant (B-NPP) is the latest example, following the Bushehr NPP (Bu-NPP), to become operational in a region with unique climatological conditions. Surprisingly, there is limited literature addressing the potential risks of radionuclide dispersion in this region, unlike the extensive research conducted in Europe, Japan, and the USA.

The main objective of this study is to investigate the spatiotemporal variability of radionuclide transport, surface concentration, and deposition after potential nuclear accidents using a mini-ensemble of meteorological inputs. The gaseous and aerosol radionuclide tracers iodine-131 (¹³¹I) and cesium-137 (¹³⁷Cs) are assumed to be released by fictitious accidents at the B-NPP in the United Arab Emirates (UAE). The amount of ¹³⁷Cs that is deposited on the surface is important due to its relatively long half-life (about 30 years), whereas the ¹³¹I surface-level concentration (about 8 d half-life) is important for human health and the biosphere in the short term (Tsuruta et al., 2019; Pisso et al., 2019; Takagi et al., 2020; Kinase et al., 2020; Wai et al., 2020).

For the dispersion modeling we use the Lagrangian particle model FLEXible PARTicle (FLEXPART) (Stohl et al., 1998, 2005; Brioude et al., 2013). Lagrangian particle dispersion models (LPDMs) model the trajectory of each particle (Nabi et al., 2015) instead of the transport of air pollutants by numerically solving the equation for the conservation of mass (Moussiopoulos, 1997; Zhang and Chen, 2007) notably done with Eulerian CTMs (Brioude et al., 2013). LPDMs are computationally more efficient as they calculate advection and diffusion only for the location of each particle and not for the entire model domain. However, LPDMs also have limitations. They suffer from numerical errors in the interpolation of meteorological fields in space and time. In some cases, the particles may not remain well-mixed during simulation (Brioude et al., 2013). This is mainly due to the treatment of the stochastic motion of the particles and/or the mass balance of vertical velocity with the horizontal winds. Regardless of the formalism used, Girard et al. (2016) have shown that uncertainties in meteorological fields, namely wind speed and direction as well as precipitation, and emission rate can significantly contribute to errors in the dispersion modeling of radionuclides. According to Galmarini et al. (2004), there are three common methods for sensitivity analysis of LPDM simulations with respect to meteorological fields: (i) generating perturbations in the horizontal and vertical position of particles, (ii) using a single meteorological model with perturbations in initial conditions and/or model physics, and (iii) using a number of different meteorological models. By adopting the third approach, we use a four-member mini-ensemble in which the performance of the forecast member, forced by 6-hourly data from the Global Forecast System (GFS), is compared to the members based on (re)analysis datasets. (Re)analysis-based simulations (unavailable in a real-world scenario) are expected to be closer to actual values than forecast-based ones (Leadbetter et al., 2022). Two members of the mini-ensemble are forced by dynamically downscaled meteorological inputs to investigate the effects of downscaling to higher resolution in dispersion modeling of radionuclides. Moreover, we investigated the sensitivity of our simulations to different turbulence schemes under convective conditions, as well as variations in the vertical profile, temporal profile, and point of emission. Specifically, we examined two different vertical profiles, two points of emission, and two temporal profiles of emissions, with one main variant and one secondary variant for each variation. The main variants were used in most of the paper, while the secondary variants were used in a small sensitivity test to evaluate the impact of these variations on our results.

Finally, to study the potential risks to the local population, the radionuclide simulations are examined in relation to the population density of the catchment region of interest (Qatar). To our knowledge, this is the first study to investigate potential radionuclide releases in the region, and our results could provide important information for developing risk assessment and preparedness plans. Specifically, our findings, which account for variations in meteorological conditions and model inputs, could help identify areas and times of the year that are more susceptible to higher concentrations of radionuclides and thus aid in developing appropriate mitigation strategies. The paper is organized as follows: Sect. 2 describes the LPDM, meteorological inputs, study area, and source term of radionuclides used in this study. The results are presented and discussed in Sect. 3. The conclusions of the study are presented in Sect. 4.

2 Mini-ensemble model configuration

A brief description of the FLEXPART modeling structure follows.

2.1 FLEXPART and FLEXPART–WRF dispersion modeling

We refer the reader to Stohl et al. (2005), Pisso et al. (2019), and Brioude et al. (2013) for detailed descriptions of the FLEXPART and FLEXPART–WRF models. Here, only the principles of FLEXPART Lagrangian modeling are discussed to facilitate the presentation of the results in the next section. FLEXPART was developed as an LPDM based on the zero acceleration scheme (Eq. 1). It solves a Langevin equation (Eq. 2) to model the trajectories of Lagrangian particles. The new position of the particles is influenced by large-scale winds, local turbulence (stochastic component), and mesoscale motions:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathit{X}\left(t+\mathrm{\Delta }t\right)=\mathit{X}\left(t\right)+v\left(\mathit{X}t\right)\mathrm{\Delta }t\\ & {\displaystyle \frac{\mathrm{d}\mathit{X}}{\mathrm{d}t}}=v\left[\mathit{X}\left(t\right)\right],\end{array}\end{array}$$

where t is time, Δt is the time increment, X is the position vector, and $\mathit{v}=\stackrel{\mathrm{‾}}{v}+v_{\mathrm{t}}+v_{\mathrm{m}}$ is the wind vector composed of the grid-scale wind $\stackrel{\mathrm{‾}}{v}$, the turbulent wind fluctuations v_t, and the mesoscale wind fluctuations v_m. FLEXPART also quantifies changes in the mass, or mixing ratio, of transported particles carried away by calculating various removal processes (Stohl et al., 2005; Grythe et al., 2017). Turbulent motions $v_{{\mathrm{t}}_i}$ for wind components i are parameterized assuming a Markov process based on the Langevin equation (Eq. 2):

$$\begin{array}{}\text{(2)}& \mathrm{d}{v}_{{\mathrm{t}}_i}={\mathit{\alpha }}_i\left(x{v}_{\mathrm{t}}t\right)\mathrm{d}t+b_{ij}\left(x{v}_{\mathrm{t}}t\right)\mathrm{d}{W}_j,\end{array}$$

where α is the drift term, b is the diffusion term, and dW_j represents incremental components of a Wiener process with mean zero and variance dt that are uncorrelated in time (Stohl et al., 2010). The minimum value of time step Δt_i is 1 s. Δt_i is used only for the horizontal turbulent wind components of Eq. (3).

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{t}_i={\displaystyle \frac{\mathrm{1}}{\mathrm{ctl}}}\min \left(\mathrm{\Delta }{\mathit{\tau }}_{L_{\mathit{\omega }}}{\displaystyle \frac{h}{\mathrm{2}\mathit{\omega }}}{\displaystyle \frac{\frac{\mathrm{0.5}}{\partial {\mathit{\sigma }}_{\mathit{\omega }}}}{\partial z}}\right)\end{array}$$

In Eq. (3), h represents the height of the atmospheric boundary layer and z denotes the height of the model level. The constant ctl represents a predefined characteristic timescale. ${\mathit{\tau }}_{L_{\mathit{\omega }}}$ is the Lagrangian timescale for the vertical velocity autocorrelation. Additionally, ω represents the turbulent vertical wind component, while σ_ω represents its standard deviation. To solve the Langevin equation for the vertical wind component, a shorter time step $\mathrm{\Delta }{t}_{\mathit{\omega }}=\frac{\mathrm{\Delta }{t}_i}{i_{\mathrm{fine}}}$ is used. Under convective conditions, when the turbulence is skewed, larger areas are occupied by downdrafts (instead of updrafts). This can lead to higher surface concentrations (deposition) in areas near pollution sources (Pisso et al., 2019). To investigate the sensitivity of radionuclide quantities to the turbulence scheme used, one sensitivity run by FLEXPART–WRF is performed after replacing the standard Gaussian turbulence model (GTM) with the skewed turbulence model (STM) (Luhar et al., 1996; Cassiani et al., 2015). The results of either scheme are used to calculate the vertical velocity component of the drift term in Eq. (2). The implementation of STM requires shorter time steps, dt in Eq. (2), to better resolve turbulence in the convective planetary boundary layer. Therefore, we use ctl=10 and i_fine=10 as recommended by Brioude et al. (2013) and Pisso et al. (2019) in this sensitivity run. For the computation of ${\mathit{\sigma }}_{v_i}$ and $\mathrm{\Delta }{\mathit{\tau }}_{L_i}$, FLEXPART uses the parameterization scheme proposed by Hanna (1982). To account for mesoscale motions, a method similar to that of Maryon (1998) is used.

In radionuclide dispersion modeling, particle mass reduction occurs primarily through three processes: radioactive decay, dry deposition, and wet deposition. The following exponential equations characterize radioactive decay (Eq. 4) and wet deposition (Eq. 5):

$$\begin{array}{}\text{(4)}& m(t+\mathrm{\Delta }t)=m\left(t\right)\exp (-\mathrm{\Delta }t/\mathit{\beta }),\end{array}$$

where m is the particle mass and the time constant $\mathit{\beta }=\frac{T_{\frac{\mathrm{1}}{\mathrm{2}}}}{\ln \left(\mathrm{2}\right)}$ is determined from the radionuclide half-life T_½.

$$\begin{array}{}\text{(5)}& m(t+\mathrm{\Delta }t)=m\left(t\right)\exp (-\mathrm{\Lambda }\mathrm{\Delta }t)\end{array}$$

The scavenging coefficient Λ is calculated differently depending on whether the particles are in the aerosol or gas phase and whether the scavenging takes place inside or below the clouds (Stohl et al., 2005).

In this study, we use FLEXPART 10.4 and a modified version of FLEXPART 9.0.2 to input meteorological simulations from the Weather Research and Forecasting (WRF) model (Grell et al., 2005) (hereafter referred to as FLEXPART and FLEXPART–WRF, respectively). Compared to previous versions, FLEXPART has undergone a significant revision of the wet deposition calculation. The latest updates incorporate the dependence of wet deposition on aerosol particle size and precipitation type (Grythe et al., 2017). FLEXPART calculates wet deposition only for cloudy grid cells wherein the precipitation rate exceeds 0.01 mm h⁻¹. Therefore, the accuracy of cloud pixel detection plays a critical role in the accuracy of the location and amount of wet deposition simulated by FLEXPART. In previous versions, including version 9.0.2 used in the development of FLEXPART–WRF 3.2, in-cloud grid cells are defined as those with precipitation and relative humidity above 90 % (Brioude et al., 2013). The grid cells below the in-cloud grid cells up to the surface are defined as below-cloud grid cells (Seibert and Arnold, 2013; Pisso et al., 2019). In recent updates to the FLEXPART's source code, the above threshold has been modified using the 3D cloud water mixing ratio (q_c) fields. The threshold of q_c>0 (q_c=0) now identifies grid cells within the cloud (below the cloud) (Pisso et al., 2019). Therefore, the differences in approaches in calculating scavenging between FLEXPART and FLEXPART–WRF may partially explain discrepancies between simulations for the wet season. Pisso et al. (2019) have discussed this and other updates in more detail in the latest version of FLEXPART. For dry deposition, for all particles below twice the reference height (h_ref), the reduction in particle mass is calculated using Eq. (6):

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }m\left(t\right)=m\left(t\right)\left[\mathrm{1}-\exp \left({\displaystyle \frac{-v_{\mathrm{d}}\left(h_{\mathrm{ref}}\right)\mathrm{\Delta }t}{\mathrm{2}{h}_{\mathrm{ref}}}}\right)\right],\end{array}$$

where v_d is the dry deposition velocity calculated as the ratio of $v_{\mathrm{d}}\left(z\right)=\frac{F_{\mathrm{C}}}{C\left(z\right)}$ for h_ref of 15 m. F_C and C(z) are the flux and concentration of a species at height z within the layer at constant flux. If the necessary information to parameterize dry deposition of gases and particles is not available, v_d can be assumed to be constant. The concentration of a given gas or aerosol species in a grid cell is equal to the weighted average of the total particle mass within the grid cell divided by the volume of the grid cell (Pisso et al., 2019), as defined in Eq. (7):

$$\begin{array}{}\text{(7)}& C\left(z\right)={\displaystyle \frac{\mathrm{1}}{v}}{\sum }_i^N\left(m_i{f}_i\right),\end{array}$$

where V is the volume of the grid cell, m_i is the particle mass, N is the total number of particles, and f_i is the fraction (weight) of the mass of particle i assigned to the particular grid cell. The amount of dry and wet deposition over the given grid cell is accumulated by default over the time dimension of the output, unless the age composition of the air parcels is required, as is the case in this study. A unique feature of FLEXPART and FLEXPART–WRF is the grouping of simulations based on the age of the Lagrangian particles. This means that simulated concentrations and deposition at each time step are obtained from aggregating the simulations over the additional dimension of particle age. In this study, we use the age spectrum of Lagrangian particles to investigate the age composition of radioactive materials from the source to the receptors. Air parcel ages are studied at an hourly resolution. As a result, the output grid, which has a horizontal resolution of 10 km and 14 vertical levels from 5 to 5000 m a.g.l., gains two dimensions: time and age spectra with the same length of 96 (hours). Thickness-weighted averages of simulated concentrations near the surface, hereafter referred to as near-surface concentrations, are calculated from concentrations within the bottom four model layers between 5 and 100 m a.g.l. (with layer thicknesses of 5, 5, 40, and 50 m). To statistically compare the age distributions of the mini-ensemble members, we use the maximum normalized difference according to Eq. (8). Assuming that a and b are two air parcel age distributions that have been smoothed with the Gaussian kernel (Chung, 2020), their maximum normalized difference is calculated as the maximum value of the absolute differences between a and b divided by the maximum value of a and b. Larger variations in the distributions are indicated by higher values of this indicator (Jin and Kozhevnikov, 2011).

$$\begin{array}{}\text{(8)}& \mathrm{maximum}\mathrm{normalized}\mathrm{difference}={\displaystyle \frac{max\left(\mathrm{abs}\left(a-b\right)\right)}{max\left(max\left(a\right),\max \left(b\right)\right)}}\end{array}$$

2.2 Meteorological data

FLEXPART and FLEXPART–WRF are driven offline by gridded meteorological fields from numerical weather prediction models. In this study, we obtain meteorological products from the GFS (NCEP, 2015a), the NCEP final analysis (FNL) (NCEP, 2015b), and the ECMWF reanalysis fifth generation (ERA5) (Hersbach et al., 2020). GFS is a fully coupled model that represents the interactions between oceans, land, and atmosphere. The GFS dataset, available since 2015, provides 6-hourly forecast inputs at 03:00, 09:00, 15:00, and 21:00 with a spatial resolution of 0.25^∘. The FNL dataset is produced using the same base meteorological model that produces the GFS dataset. The former provides the 3-hourly combination of analysis (at 00:00, 06:00, 12:00, and 18:00) and forecast meteorological fields (at 03:00, 09:00, 15:00, and 21:00) at a spatial resolution of 0.25^∘ starting in 2015. ERA5 is the latest generation of ECMWF reanalysis data, covering the period from 1 January 1950 to the present. It is generated in hourly time steps with a spatial resolution of about 0.25^∘. GFS and FNL are used to force FLEXPART directly. In addition, WRF 4.2 is applied to dynamically downscale the ERA5 and FNL to be used as the input to FLEXPART. They are used in FLEXPART–WRF with higher downscaled spatiotemporal resolutions of 10 km and hourly. Though we devoted efforts to incorporating ERA5 reanalysis data directly as inputs for FLEXPART, we were hindered by technical difficulties encountered with the data preprocessor tool, flex_extract, thereby impeding our ability to do so. While we were working to resolve the technical issues with the developers, a solution was not readily available during the writing and revision of this paper. We have indirectly employed ERA5 reanalysis data from the Research Data Archive (RDA) by incorporating it into the WRF model and subsequently into the FLEXPART–WRF. This approach has enabled us to assess the influence of ERA5 on our downscaled simulations and to compare it with FNL. To the best of our knowledge, this is the first time that the ERA5 reanalysis has been used with the FLEXPART–WRF model setup. Table 1 provides an overview of the WRF model configuration.

Table 1 Overview of the WRF model configuration.

Table 2 provides an overview of the input data and the corresponding simulation codes. For comparison of the above meteorological inputs with observations, daily total precipitation and daily averages of wind speed and temperature data are obtained from 157 climate stations within the model domain. These observations are freely available from the Global Surface Summary of the Day (GSOD) database, which is maintained by the National Oceanic and Atmospheric Administration (NOAA). The GSOD data are subject to a rigorous two-tier quality control (QC) process (https://www.ncei.noaa.gov/data/global-summary-of-the-day/doc/readme.txt, last access: 5 May 2023). The first level entails the application of meticulous automated QC procedures that effectively purge the raw data of random errors. The second level of QC is performed to create daily summaries. However, there is still a slight likelihood of unknown errors being present within the GSOD dataset.

Table 2 Summary of meteorological inputs used to run FLEXPART and FLEXPART–WRF.

2.3 Emission scenario and study area

We simulate fictitious nuclear accidents at the B-NPP (Fig. 1a) in which 6.9 kg (22 PBq) of ¹³⁷Cs and 0.042 kg (192 PBq) of ¹³¹I are released during the first 24 h of each 96 h simulation period. These amounts are the upper bounds of emissions estimated by Babukhina et al. (2016) for the Fukushima Daiichi nuclear accident in March 2011. Although we simulate a fictitious release of radioactivity at a level comparable to the Fukushima nuclear accident, our study does not intend to replicate past accidents or simulate a specific real-world case. These source terms are used to provide a real-world comparison that gives the reader a tangible point of reference.

The number of Lagrangian particles required is dictated by the specific problem at hand (Papagiannopoulos et al., 2020; Thompson et al., 2015; Fast and Easter, 2006). The computational timescales increase linearly with the increasing number of particles, while the statistical error of the simulations decreases with the square root of the particle density (Pisso et al., 2019). In this study, a total of 10⁴ Lagrangian particles are released uniformly during the 24 h emission period. To assess if this rate is sufficient, we study the dry deposition process, which is directly influenced by the boundary layer conditions. We performed a preliminary test run with GFS data with a 10-fold increase in particles. As shown in Fig. S1, the simulated dry deposition of ¹³⁷Cs and ¹³¹I and the wet deposition of ¹³⁷Cs do not undergo significant changes with the increase of an order of magnitude over Qatar.

Regarding the height of the releases, the particles are emitted at model levels between 100 and 300 m above ground level (a.g.l.) above the emission point. Each simulation starts at the beginning of each day and lasts 96 h (Fig. 1b). We run this scenario for every day of the year 2019, resulting in a total of 1460 simulation days (365 d of the year × 4 forward simulation days). In the following analysis, when reporting particle ages associated with the simulations, we always refer to the time when a particle is released.

As expected, the temporal and spatial characteristics of the release of radioactive materials (like other pollutants) significantly impact the simulations. In this study, we perform three sensitivity analyses on the release conditions, in addition to one sensitivity test on the skewed turbulence model (STM), using FLEXPART–WRF based on ERA5 reanalysis. Specifically, we conduct the continuous release sensitivity test (CRST), whereby pollutants are released continuously for 96 h, as opposed to only the first 24 h of the accident. Following Morino et al. (2011), we also increase the upper limit of release height from 300 to 700 m (referred to as the release height sensitivity test, RHST). The third sensitivity test, designated as the release location sensitivity test (RLST), is conducted by releasing radionuclides from the Bushehr NPP (Bu-NPP) and investigating the impact of source–receptor position on the transport of radionuclides, particularly in relation to seasonal variations in atmospheric patterns. The simulation domain extends between 17–33^∘ N and 40–60^∘ E, and for post-processing we analyze simulations over Qatar between 24.25–26.35^∘ N and 50.65–51.75^∘ E.

Figure 1 (a) The study area includes the Barakah nuclear power plant (B-NPP) and the Bushehr new nuclear power plant (Bu-NPP), both represented by red squares, as well as the country of Qatar. The base map and overlay information, excluding the names of the NPP facilities, are sourced from Google Earth. The LPDM simulation cycle is schematically depicted in (b), which is originally retrieved from https://www.janicke.de/en/lasat.html (last accessed: 5 May 2023).

3 Results and discussion

3.1 Meteorological inputs vs. observations

We compare the relevant meteorological input fields including daily total precipitation, daily mean surface wind speed, and temperature with observations from 157 stations spanning the model domain (Table 3). A relatively large discrepancy between observations and model inputs (as well as among different model inputs, as depicted in Figs. S2 and S3) is obtained for precipitation. The average Spearman correlation coefficient (r) of ∼ 0.41 and root mean square error (RMSE) of 3.2 mm indicate a moderate relationship between simulated precipitation inputs and observations in all members. The ERA5–WRF and FNL datasets have the highest correlation (0.44) and lowest RMSE (3 mm) compared to observations, respectively. All datasets, especially the downscaled ones, underestimate precipitation amounts, resulting in an average mean bias error (MBE) of −0.22. The seasonal box plots of precipitation datasets (Fig. S4a) reveal that the underestimation of precipitation is most pronounced in extreme precipitation events, particularly in the downscaled simulations. Specifically, the upper quartile of precipitation in observations during the rainy seasons of winter and fall is approximately 7.5 and 5.5 mm, respectively. In contrast, the corresponding values in the FNL, GFS, ERA5–WRF, and FNL–WRF simulations are markedly lower, amounting to approximately 6 and 3.1 mm, 5.8 and 3.3 mm, 4.8 and 2.3 mm, and 2 and 1 mm, respectively.

The underestimation of precipitation in atmospheric modeling datasets is a well-known challenge, attributed to the complex interplay of subgrid-scale processes that are often parameterized in the models. These parameterizations rely on assumptions and approximations, introducing uncertainties in precipitation simulations. Additionally, the intricate topography of the study area can significantly influence precipitation patterns (Tapiador et al., 2019).

The more severe underestimation observed in downscaled simulations can be partly attributed to the greater sensitivity of these simulations to the choice of physical parameterizations. The use of higher spatial resolution in downscaled simulations means that physical parameterizations are applied over smaller areas, which can magnify the impact of uncertainties in these parameterizations. For instance, parameterizations of cloud microphysics and convection can be particularly important in determining the intensity and duration of extreme precipitation events (Kain, 2004). While this study did not specifically investigate model resolutions finer than 10 km, increasing the model resolution to 1 km or less provides the capability to resolve finer subgrid phenomena. This allows for the deactivation of convective schemes designed for coarser resolutions, as the resolved deep convection becomes more prominent.

All meteorological datasets are closer to the observations when evaluating the simulations of wind speed and temperature. The ERA5–WRF dataset shows the best agreement with the wind speed observations in terms of correlation (0.65) and RMSE (1.5 m s⁻¹), but overestimates wind speed (MBE > 0.5) to a greater extent than the GFS and FNL datasets (MBE < −0.13). Based on the box plot of wind speeds shown in Fig. S4b, the GFS and FNL datasets underestimate wind speed mainly in low intensities, while the ERA5–WRF and FNL–WRF overestimate wind speed in high intensities. All datasets are significantly correlated with the temperature observations and with each other, as shown by the distribution of data points along the identity line (Figs. S2 and S3). The distributions of the simulated and modeled temperature exhibit the least dissimilarity (Fig. S4c) in comparison to precipitation and wind speed. In order to deepen our understanding of the model performance at the primary emission source and receptors in Qatar, we conducted a similar analysis on 12 stations within the area encompassing the B-NPP and the state of Qatar (Figs. S5, S6, and S7). Although some variations are noted, such as an improvement in the correlations between simulated and observed wind speed, the overall patterns of performance for the input datasets over this region remain consistent with those identified earlier. Based on these findings, we infer that the FNL–WRF and ERA5–WRF (downscaled) datasets exhibit a better correlation with observations. However, this comes at the expense of increased bias values when compared to the FNL and GFS datasets.

The intercomparison of the input datasets shows that the FNL and GFS simulations have the best agreement, as can be expected since they are produced using the same data assimilation and forecasting system. The better agreement found between the ERA5–WRF and FNL–WRF datasets than between the latter and the FNL dataset shows the impact of the downscaling on the homogeneity of the meteorological inputs. It should be noted that a systematic comparison of meteorological datasets with a representative sample of surface observations is needed to determine the optimal choice of meteorological inputs for forcing a CTM for different regions. Here we have access to daily surface meteorological data, while transport modeling is performed in several layers of the atmosphere at hourly time steps or shorter.

Table 3 Comparison of modeled and observed daily precipitation, wind speed, and temperature. The bold font is utilized to highlight the best statistical value contributed by each mini-ensemble member for every meteorological variable.

3.2 Age composition of radionuclide plumes

The seasonal distributions of the air parcel ages for ¹³¹I and ¹³⁷Cs concentrations are shown in Figs. 2 and S8, respectively. The high degree of similarity between the age distributions of ¹³¹I and ¹³⁷Cs is due to the fact that the removal process in LPDMs only affects mass concentrations and not particle positions. Therefore, the following discussion of the ages of the ¹³¹I particles is also valid for the ages of the ¹³⁷Cs particles. All mini-ensemble members, in particular GFS and FNL, predominantly simulate air parcels with lifetimes of less than 20 to 30 h. Other than fall, the age distribution of the particles is relatively similar in all other seasons. Seasonal atmospheric patterns with a low frequency of occurrence in the GFS and FNL datasets, discussed in Sect. 3.3, seem to be responsible for the difference in particle age distributions in fall. Compared to the other mini-ensemble members, the FNL-based simulations show a delayed appearance of the second peak in air particle ages. This could be due to a decrease in transport speed to the same receptors and/or a higher number of air parcels reaching areas away from the source (discussed later).

To understand the relationship between radionuclide concentrations and air parcel ages, the latter are divided into three groups corresponding to low (the top row in Fig. S9), moderate (the bottom row in Fig. S9), and high (the bottom row in Fig. 2) near-surface concentrations of ¹³¹I (see Figs. S8 and S10 for ¹³⁷Cs). These categories are determined for each member based on concentrations below the 33rd percentile, between the 33rd and 66th percentiles, and above the 66th percentile, respectively. The rapid transport of dense radionuclide clouds within the region of interest is demonstrated by the strongly positively skewed age distributions of particles corresponding to moderate and high radionuclide concentrations in all members and all seasons. The inconsistency of age distributions for high and moderate concentrations in summer can be attributed to an unusual transport pattern in this season.

It is worth noting that all members of the mini-ensemble show closer agreement for air parcel ages corresponding to the moderate and high levels of ¹³¹I and ¹³⁷Cs concentrations than for those corresponding to low concentrations. This is consistent with the principles of FLEXPART dispersion modeling. According to Pisso et al. (2019), the error rate of the simulations decreases with the square root of the particle density. As a result, all four members simulate the age of Lagrangian particles corresponding to moderate and high concentrations in better agreement and likely more accurately near the source. Using the maximum normalized distance, a quantitative comparison of the particle age distributions, corresponding to all concentrations, is shown in Fig. 3. The age distribution in the ERA5–WRF- and FNL–WRF-based simulations shows a greater similarity than to that of the FNL-based simulations in all seasons, except fall. This is due to the dynamic downscaled meteorological inputs with the same spatiotemporal resolution and to a common simulation code in ERA5- and FNL–WRF-based simulations. The age distributions based on the FNL and GFS inputs, which are generated by a similar base model, do not show the same degree of similarity in all seasons. They have the smallest difference in fall and spring and a larger difference in winter and summer.

Figure 2 Distributions of the air parcels corresponding to all near-surface ¹³¹I concentrations (top row) and of those above the 66th percentile (bottom row). The y and x axes show density values and air parcel ages, respectively.

Figure 3 Seasonal maximum normalized difference of air parcel age distributions.

Due to the relatively long half-life of ¹³⁷Cs, its deposition rate in the affected areas is of great importance. To analyze the relationship between the age composition of air parcels and the amount of ¹³⁷Cs deposition, the deposition values cumulatively aggregated across time steps (j) and age spectra (i) are normalized to the total amount of ¹³⁷Cs deposition in each grid cell (k) at the end of each simulation run (l).

$$\begin{array}{}\text{(9)}& {}^{\mathrm{137}}{\mathrm{Cs}}_{{\mathrm{klna}}_{\mathrm{norm}\mathrm{\_}\mathrm{depso}}}=\left\{\begin{array}{l}\frac{{}^{\mathrm{137}}{\mathrm{Cs}}_{kl(n-\mathrm{1})n}}{{\sum }_{j=\mathrm{1}}^{\mathrm{96}}{\sum }_{i=\mathrm{1}}^{\mathrm{95}}{}^{\mathrm{137}}{\mathrm{Cs}}_{klij}}\mathrm{if}n=\mathrm{2}\\ \frac{{\sum }_{j=\mathrm{1}}^{n-\mathrm{1}}{\sum }_{i=\mathrm{1}}^{j-\mathrm{1}}{}^{\mathrm{137}}{\mathrm{Cs}}_{klij}+{\sum }_{i=\mathrm{1}}^a{}^{\mathrm{137}}{\mathrm{Cs}}_{klin}}{{\sum }_{j=\mathrm{1}}^{\mathrm{96}}{\sum }_{i=\mathrm{1}}^{\mathrm{95}}{}^{\mathrm{137}}{\mathrm{Cs}}_{klij}}\mathrm{if}n>\mathrm{2}\end{array}\right.\end{array}$$

Here, n is the time step with a maximum of 96 (the last time step) and a is the given particle age with a maximum of n−1. ¹³⁷Cs_klij=0 in two conditions: (1) n>=26 if iϵ [1, n−25] and (2) i>=j.

Figure 4 shows the normalized deposition amounts (${}^{\mathrm{137}}{\mathrm{Cs}}_{kln{a}_{\mathrm{norm}\mathrm{\_}\mathrm{depso}}}$) in winter, when both dry deposition and wet deposition occur in the study area. Similar deposition patterns are obtained for other seasons (Fig. S11). As shown in Fig. S12, the main reason for the small difference between the seasonal deposition patterns is the lack of precipitation and subsequent wet deposition in the region. Although the spatial pattern of the deposition varies considerably, as indicated by the range of quartiles, the median of the normalized deposition shows that about 80 % of the deposition occurs within 80 h after an accident. The cumulative deposition at the end of the simulation period is mostly in the areas farthest from the source, and the total deposition reaches 100 % as it approaches the end of the 96 h simulation period. This is evident in Fig. S13, where the ages of the deposited particles peak around 20–30 h, with a rather small number after 80 h.

Figure 4 The normalized cumulative ¹³⁷Cs deposition compared to the total ¹³⁷Cs deposition of each 96 h simulation period at each grid cell in winter. The x axis shows the age of the Lagrangian particles, and the y axis is the normalized deposition. The boxes show the quartiles of the normalized deposition. The whiskers represent the range between 1.5 times the interquartile range above the upper quartile and below the lower quartile.

3.3 Spatiotemporal distribution of radionuclides

In this section, we explore the seasonal and diurnal changes in the concentrations and deposition of radionuclides within the study area. To analyze the radionuclide concentrations, we used the average of simulations in the lowest four layers of the model between 5 and 100 m. Given that ¹³¹I is highly radioactive and that ¹³⁷Cs has a high solubility and a relatively long half-life, we will focus on their concentrations and deposition, respectively. Figures 5 and 6 illustrate the seasonal median of the 96 h integrated ¹³¹I concentrations (¹³¹I^{intg_conc_seas}, in units of ng m⁻³) and of the maximum (total) ¹³⁷Cs deposition (¹³⁷Cs^{tot_depos_seas}, in units of ng m⁻²) within the study area.

Our mini-ensemble simulations show that the ¹³¹I^{intg_conc_seas} values above $\mathrm{1.43}\times {\mathrm{10}}^{-\mathrm{3}}$ ng m⁻³ occur frequently in the cold period of the year, especially in fall, in the simulations of all members. This may be due in part to the lower (higher) planetary boundary layer height (PBLH) in the cold (warm) seasons and the synoptic conditions. An exception is the distribution of ¹³¹I^{intg_conc_seas} in the FNL-based simulations during the summer. This may be caused by a rare atmospheric circulation in summer, discussed later, and it may change as the modeling period is extended. In terms of spatial distribution, ¹³¹I^{intg_conc_seas} above $\mathrm{1.43}\times {\mathrm{10}}^{-\mathrm{3}}$ ng m⁻³ occurs close to the source in the southeastern part of the domain, and the intensity of the ¹³¹I^{intg_conc_seas} decreases with distance to the north. In other words, southeastern Qatar is the first area to be affected by dense ¹³¹I clouds in the event of a nuclear accident, especially during the cold period of the year. The advance of ¹³¹I^{intg_conc_seas} above $\mathrm{1.43}\times {\mathrm{10}}^{-\mathrm{3}}$ ng m⁻³ in the south of Qatar in simulations based on FNL inputs occurs in both fall and winter, but is found only in winter in the GFS-based simulations. ¹³¹I^{intg_conc_seas} values peak in the fall in both downscaled runs, with inputs from the ERA5–WRF and FNL–WRF datasets, but the extent of high ¹³¹I^{intg_conc_seas} values in the southeast of Qatar in fall is much larger in the ERA5–WRF run. This is also the case when comparing the FNL–WRF and FNL-based simulations. The differences have resulted in ¹³¹I^{intg_conc_seas} values varying by a factor of 2 to 10 in the south of the area of interest between mini-ensemble members.

To investigate the influence of the particle release time on the radionuclide dispersion, the seasonal median of the particle release time (hours in local time, LT) coinciding with the maximum concentration of ¹³¹I and the completion of ¹³⁷Cs deposition is considered (contours in Figs. 5 and 6). The results show that the highest ¹³¹I concentrations coincide with particles released between 09:00 and 14:00 LT in most parts of the study area. This is the time of day when the development of the planetary boundary layer, the intensification of the land–sea thermal gradient, and the resulting daytime onshore winds coincide in the region. Among the members, the earliest particle release times (between 06:00 and 09:00 LT) leading to the highest ¹³¹I concentrations are observed in simulations based on ERA5–WRF data.

All members have simulated ¹³⁷C deposition above 10 ng m⁻² over a significant part of the study area during the cold period of the year. The highest ¹³⁷Cs^{tot_depos_seas} values are observed in the south and southeast, near the emission point. Compared to the FNL- and GFS-based simulations, the simulations based on the FNL–WRF and ERA5–WRF datasets show the much greater extent of the ¹³⁷C deposition over almost the whole of Qatar in the winter and, to a lesser extent, in the fall. Among the mini-ensemble members, the highest ¹³⁷Cs^{tot_depos_seas} in the southeast is seen in the simulations based on ERA5–WRF inputs, followed by the FNL–WRF-based simulations, in fall. While the ERA5–WRF-based simulations of ¹³⁷Cs^{tot_depos_seas} exceed 30 ng m⁻², the deposition in GFS-based simulations is up to 10 times lower in some points of the same area and same period. The ¹³⁷Cs^{tot_depos_seas} in the warm period of the year, except for the simulations based on the FNL dataset in summer, is mostly either close to or below 10 ng m⁻². The higher ¹³⁷Cs^{tot_depos_seas} in the cold period of the year, especially in winter, can be largely attributed to the seasonal increase in ¹³⁷Cs transport by southerly winds, discussed later, and to the relative increase in wet deposition. The above results indicate that if a nuclear accident occurs during the cold period of the year, the magnitude of extreme ¹³¹I concentrations and ¹³⁷Cs deposition in the south of Qatar may be up to 10 times greater than during the warm seasons. A similar order of magnitude in simulated variability is seen across the mini-ensemble members. With respect to the particle release times, ¹³⁷Cs^{tot_depos_seas} occurs mainly when particles are released between 10:00 and 17:00 LT in the presence of a turbulent boundary layer.

The analysis of the age spectra in Sect. 3.2 shows the higher frequency of long-lived air parcels in the FNL-based simulations (see top row in Fig. 2). To further examine these results, the spatial distributions of the full-year median of the 96 h integrated ¹³¹I concentrations (¹³¹I^{intg_conc_full}) and the full-year median of the air parcel ages coinciding with the maximum concentration of ¹³¹I found in each 96 h run are examined (Fig. 7). As expected, the age of the Lagrangian particles decreases southward with proximity to the source. All mini-ensemble members simulate particle ages lower than 25 h along with relatively high levels of ¹³¹I^{intg_conc_full} at the southeastern corner of the study area. The longer-lived particles (above 35 h) are found in higher latitudes. In simulations based on the FNL dataset, particle ages in northern Qatar exceed 40 to 50 h, while they mostly do not exceed 40 h in other mini-ensemble members. Furthermore, we obtain the larger extent of relatively higher ¹³¹I^{intg_conc_full} to central Qatar in these simulations (considering both ¹³¹I^{intg_conc_seas} and ¹³¹I^{intg_conc_full} shown in Figs. 5 and 7). Our analysis suggests that the longer-lived particles in FNL-based simulations can be attributed to a higher number of air parcels reaching areas away from the source, rather than a lower transport speed.

Figure 5 Seasonal median of 96 h integrated ¹³¹I concentrations (¹³¹I^{intg_conc_seas}). The contour lines (in local time, hours of the day) depict the seasonal median of the Lagrangian particle release time (local time) coinciding with the maximum ¹³¹I concentrations found in each 96 h run.

Figure 6 Same as Fig. 5, but for the seasonal median of total ¹³⁷Cs deposition (¹³⁷Cs^{tot_depos_seas}) (units of ng m⁻²). The contour lines are the seasonal median of the Lagrangian particle release time (local time) coinciding with the completion of ¹³⁷Cs deposition found in each 96 h run.

¹³¹I^{intg_conc_seas} and ¹³⁷Cs^{tot_depos_seas} are analyzed against the population density of Qatar (Fig. 8). The desert areas of southern and southeastern Qatar, over which radionuclides enter the country, host a small population or are almost uninhabited (Fig. 8a). Figure 8b and c show that extremely high levels of ¹³¹I^{intg_conc_seas} (greater than $\mathrm{0.6}\times {\mathrm{10}}^{-\mathrm{2}}$ ng m⁻³) and ¹³⁷Cs^{tot_depos_seas} (greater than 60 ng m⁻²) occur mostly in areas with a population density of fewer than five persons per arc second. In the populated areas (with a density of more than 8 persons per arc second) ¹³¹I^{intg_conc_seas} and ¹³⁷Cs^{tot_depos_seas} do not exceed $\mathrm{0.4}\times {\mathrm{10}}^{-\mathrm{2}}$ ng m⁻³ and 40 ng m⁻² in most cases. Due to the exceptional weather pattern that occurs in the simulations based on the FNL dataset in the eastern part of Qatar (where the most densely populated areas are located) in summer, these simulations cause the highest values of ¹³¹I^{intg_conc_seas} in densely populated areas. Otherwise, all mini-ensemble members simulate the highest ¹³¹I^{intg_conc_seas} during the cold seasons. The highest ¹³⁷Cs^{tot_depos_seas} in the same areas is found based on ERA5 and FNL–WRF in the cold period of the year.

Figure 7 The full-year median of 96 h integrated ¹³¹I concentrations (¹³¹I^{intg_conc_full}). The contour lines are the full-year median of age spectra coinciding with the maximum ¹³¹I concentrations found in each 96 h run.

Figure 8 (a) The gridded population density of Qatar at 1 arcsec resolution in 2020. The relationship between ¹³¹I^{intg_conc_seas} (ng m⁻³) and ¹³⁷Cs^{tot_depos_seas} (ng m⁻²) as well as the population density of Qatar are shown in (b) and (c), respectively. The shapes and colors of the markers represent seasons and mini-ensemble members.

In addition to the magnitude of the transported radioactive materials, the temporal distribution of extreme events in the affected area is also of importance for preparedness programs. Figure 9a shows the frequency of occurrence (FoO) of ¹³¹I concentrations above the 66th percentile. In all members of the mini-ensemble, more than half of the events, especially in the northern part, take place in the winter. While the FoO of high ¹³¹I varies between around 15 % and 30 % in spring, the lowest FoO is observed in the summer when less than about 10 % of the extreme events occur. The seasonal FoO of ¹³⁷Cs concentrations also shows very similar results (Fig. S14). A case in which ERA5–WRF simulated the northwestward movement of the near-surface ¹³¹I concentrations (ng m⁻³) in January at noon is shown in the upper panel of Fig. 9b. We find a similar pattern in a large number of events in which high levels of radionuclides are transported to Qatar. This synoptic pattern is related to the juxtaposition of low- and high-pressure cells located to the west and east of the region. The resulting pressure gradient causes strong south-southeast winds to develop between two cyclonic and anticyclonic cells, bringing the dense ¹³¹I clouds into the study area. This pattern mainly occurs in the late winter and early spring, coinciding with the southward movement of the westerlies and the eastward movement of the Saudi Arabian subtropical high-pressure system (De Vries et al., 2016). The lower panel of Fig. 9b shows the summertime mean of near-surface ¹³¹I concentrations (ng m⁻³), obtained from simulations based on the ERA5–WRF database. The near-surface atmospheric circulation is superimposed on these simulations. The seasonal pattern found here illustrates why very few extreme events are observed in the study area during the summer. The northwest–southeast winds, known as the shamal winds (Yu et al., 2016), cause the simulated radionuclides to move away from the study area. Therefore, the highest (lowest) contamination risk to the population of Qatar from the occurrence of a nuclear accident in B-NPP and the resulting release of radionuclides is expected to occur in winter (summer).

Figure 9 (a) Frequency of occurrence (%) of ¹³¹I concentrations above the respective 66th percentile. (b) The panel at the top shows the simulation of near-surface¹³¹I concentrations (ng m⁻³) based on the ERA5–WRF dataset on 14 January 2019. The panel at the bottom is the summer average of near-surface ¹³¹I concentrations. The wind field (m s⁻¹) at 12:00 on the first day of the 96 h simulation period is used to represent the atmospheric circulation during the radionuclide transport.

3.4 Sensitivity of simulations to turbulence scheme, emission profiles, and emission location

This section focuses on investigating the sensitivity of radionuclide simulations using ERA5–WRF inputs through four distinct sensitivity runs. Each of these runs explores specific variations in model setup, such as turbulence scheme, emission duration, emission height, and emission location. To establish a baseline, all sensitivity runs are compared against the control run that utilizes the GTM as the turbulence scheme and emission characteristics that are explained in Sect. 2.3. Results show that the use of the skewed turbulence model (STM) leads to a more frequent occurrence of high ¹³¹I concentrations and ¹³⁷Cs deposition within Qatar (Fig. 10). The upper quartiles of the STM-based ¹³¹I concentrations are around 10 % and 16 % higher than those of the GTM-based control simulations in winter and year-round, respectively. The upper quartiles of STM-based deposition simulations also increase by 23 % and 12 % over the same periods in comparison to GTM-based control deposition simulations. As noted previously, the elevation in ¹³¹I concentrations and ¹³⁷Cs deposition based on the STM is not surprising, given that the STM promotes more downdrafts than updrafts. This pattern leads to higher surface concentrations and deposition in the proximity of pollution sources (Pisso et al., 2019).

In the CRST run with the extension of the release duration from 24 to 96 h, a decrease in ¹³¹I concentrations and the amount of ¹³⁷Cs deposition compared to the control run is observed, except for ¹³¹I concentrations in spring. This result is to be expected, as the particles released at the end of the 96 h simulation period have less time to reach the study area, resulting in a considerable decrease in the concentration and deposition of radioactive materials. For instance, this implementation leads to a reduction of approximately 75 % and 45 % (170 % and 140 %) of the upper quartiles of ¹³¹I concentrations (¹³⁷Cs deposition) in winter and fall compared to the control run.

Figure 10 (a) These box plots depict the 96 h integrated simulations of ¹³¹I concentrations (ng m⁻³) based on ERA5–WRF inputs. The simulations include control runs shown in cyan, sensitivity runs with the skewed turbulence model (STM) in brown, increased release duration from 24 to 96 h (CRST) in blue, modified release height from 100–300 to the 100–700 m a.g.l. (RHST) in pink, and a change in the release location from Barakah to Bushehr nuclear power plant (RLST) in green. (b) The box plots show the same for the total deposition of ¹³⁷Cs (ng m⁻²). The quartiles of the simulations are depicted by the box borders, while the whiskers illustrate the range between 1.5 times the interquartile range above the upper quartile and below the lower quartile.

Morino et al. (2011) used the Eulerian chemical transport model, CMAQ, to investigate the sensitivity of ¹³¹I and ¹³⁷Cs deposition to the vertical profile of emissions from the Fukushima nuclear accident. They concluded that the fractions of ¹³¹I and ¹³⁷Cs deposition were insensitive to the change in emission height. In contrast, our simulations indicate that an increase in emission height leads to a significant decrease in ¹³¹I concentrations and ¹³⁷Cs deposition, except in the summer. Specifically, the upper quartile of ¹³¹I concentrations in the RHST run decreases by approximately 25 % and 130 % in the winter and fall, respectively, compared to the control run. Similarly, the increase in emission height is associated with a decrease in ¹³⁷Cs deposition, with the upper quartile of simulations showing a considerable decrease of up to 65 % in RHST simulations compared to the control implementation during the fall. The severity of deposition reduction lessens during the winter and spring seasons and is instead accompanied by a slight increase during the summer. It seems that increasing the release height results in particle entry into the upper layers of the atmosphere, especially during the cold months when regional and larger-scale atmospheric patterns are more prevalent. This results in the transportation of the majority of radionuclides to more distant locations, leading to a decrease in ¹³¹I concentrations and ¹³⁷Cs deposition.

The RLST run provides insights into the impact of shifting the emission point from B-NPP to Bu-NPP. The results demonstrate notable disparities in the transportation of radioactive materials to Qatar. Specifically, the results reveal that in the event of a nuclear accident at Bu-NPP, there is a relatively reduced transfer of radioactive materials to Qatar compared to emissions originating from B-NPP. This trend holds true for most of the year. This outcome can indeed be attributed to the greater distance between Bushehr NPP and Qatar, as anticipated. The upper quartile of the full-year median of ¹³¹I concentrations and ¹³⁷Cs deposition in RLST simulations decrease by a factor of approximately 2.2 and more than 3, respectively. The seasonal variability of the transport and deposition of radioactive materials from Bu-NPP is a noteworthy aspect to consider in our analysis. The simulations of the RLST run show that ¹³⁷Cs deposition is more prominent in winter and spring compared to other seasons. Moreover, the upper quartile value of ¹³¹I concentrations from Bu-NPP exceeds that of B-NPP by more than 80 % in spring, which is in contrast to the simulations of other seasons.

In Fig. 11a, a distinct spatial pattern emerges in the Bu-NPP simulations when compared to the B-NPP simulations (Figs. 5, 6, and 7). The simulated ¹³¹I^{intg_conc_full} exhibits relatively high values in the northern and eastern regions, including Doha, particularly during the winter and spring seasons. This occurrence can be attributed to the downward movement of westerlies, facilitating the transport of air masses enriched with ¹³¹I towards these areas. The full-year median of the total ¹³⁷Cs deposition (¹³⁷Cs^{tot_depos_full}) (Fig. 11b) closely follows this temporal pattern, showing prominent peaks during the spring and winter seasons within the northern and eastern regions of Qatar. Conversely, during the fall and summer seasons, the levels of ¹³¹I^{intg_conc_full} and ¹³⁷Cs^{tot_depos_full} are greatly reduced. The age of particles corresponding to high levels of ¹³¹I and ¹³⁷Cs simulations indicates their earlier entry into the northern half of Qatar during these two seasons, reaching the region within approximately 40–50 and 60 h after release, respectively.

Figure 11 Simulated ¹³¹I^{intg_conc_full} (a) and ¹³⁷Cs^{tot_depos_full} (b) based on ERA5–WRF inputs originating from Bushehr NPP. The contour lines are the full-year median of age spectra coinciding with the maximum ¹³¹I concentrations and the completion of ¹³⁷Cs deposition found in each 96 h run.

3.5 Intercomparison of mini-ensemble members

Radionuclide simulations based on (re)analysis datasets are assumed to be a better approximation of the actual atmospheric conditions than those based on the forecasting dataset (GFS in this study) (Leadbetter et al., 2022). The statistical evaluation metrics used here are recommended by Maurer et al. (2018), who assessed transport models simulating Xe-133 using measurements from six International Monitoring System (IMS) stations. In addition to other common evaluation metrics, they used the fractional bias (FB) and the fraction within a factor of 5 (F5). FB, in the range −2 and 2, is the bias of the simulated mean values normalized by the sum of the simulation and measurement means multiplied by 2. F5 is the fraction of simulations that are at most 1 factor larger (5) or smaller (0.2) than the reference values. The Spearman correlation coefficient (r) calculated between the simulations of ¹³¹I concentrations shows that the GFS-based simulations are closely associated with those of the (re)analysis-based members (red circles in the bottom row in Fig. 12). The GFS-based simulations obtain the highest correlation with the FNL-based simulations (0.85) followed by FNL–WRF-based simulations (0.73). While the GFS-based simulations have (on average) a small FB (0.07) compared to the FNL-based simulations, they are positively biased compared to the simulations based on the FNL–WRF (−0.99) and ERA5–WRF (−0.98) datasets. According to F5, 74.06 % and 55.69 % of the GFS-based simulations are within a factor of 5 of the FNL and FNL–WRF-based simulations, respectively, whereas F5 decreases to 54.72 % between the GFS- and ERA5–WRF-based simulations. The RMSE between simulations based on the GFS and FNL datasets ($\mathrm{1.49}\times {\mathrm{10}}^{-\mathrm{2}}$ ng m⁻³) is smaller than that found between the former and the FNL–WRF- (0.0235 ng m⁻³) and ERA5–WRF-based ($\mathrm{2.41}\times {\mathrm{10}}^{-\mathrm{2}}$ ng m⁻³) simulations. As for the other metrics, the GFS-based simulations produce the lowest and highest NMSE against the FNL (1.24) and ERA5–WRF (8.78) simulations, respectively. In short, the GFS- and FNL-based simulations have the highest agreement due to the large similarities between their meteorological inputs (Sect. 2.2). Among all mini-ensemble members, the largest difference occurs between simulations based on ERA5–WRF and FNL datasets (r=0.55, FB = −0.92, F5 = 50.6 %, RMSE = $\mathrm{2.62}\times {\mathrm{10}}^{-\mathrm{2}}$ ng m⁻³, and NMSE = 11.18). On the other hand, the downscaling of the inputs and the application of the same simulation code yield higher agreement of the simulations based on the FNL–WRF and ERA5–WRF datasets than with those based on the FNL and FNL–WRF datasets (r=0.7 vs. 0.67, FB = 0.01 vs. −0.93, F5 = 69.83 % vs. 53.3 %, RMSE = $\mathrm{0.78}\times {\mathrm{10}}^{-\mathrm{2}}$ ng m⁻³ vs. $\mathrm{2.48}\times {\mathrm{10}}^{-\mathrm{2}}$ ng m⁻³, and NMSE = 2.74 vs. 10.08). These results point out that particle dispersion modeling is primarily influenced by both the meteorological inputs and the dispersion model of choice, which is consistent with the results of Karion et al. (2019).

Figure 12 The intercomparison of the 96 h integrated simulations of near-surface ¹³¹I (red circles) and ¹³⁷Cs (blue squares) concentrations (ng m⁻³) simulated at each grid point in each day of 2019. Solid and dotted lines show regression and identity lines, respectively. The asterisk next to the Spearman correlation coefficient (r) indicates a statistically significant correlation at p<0.05.

The evaluation metrics of ¹³⁷Cs concentrations (blue squares in the top row in Fig. 12) lead to the same results as discussed above. For example, the simulations based on the GFS dataset show the highest agreement with the FNL-based simulations (r=0.85, FB = 0.07, F5 = 74.98 %, RMSE = 2.64 ng m⁻³, and NMSE = 1.19). Consistent with our findings, the simulations derived from FNL and ERA5–WRF datasets exhibit the most robust agreement, while the simulations based on FNL–WRF and ERA5–WRF datasets demonstrate the least agreement among the ensemble members relying on the (re)analysis inputs. The GFS-based simulations mostly lead to closer agreement with those based on the (re)analysis datasets when comparing the ¹³⁷Cs concentration simulations than the ¹³¹I concentration simulations.

The assessment of the simulated ¹³⁷Cs deposition reveals a notable lack of agreement between all mini-ensemble members, as depicted in Fig. 13. In particular, the disparities in wet deposition rates are considerably more pronounced than those in dry deposition rates. This can be largely explained by the high degree of uncertainty associated with the inputs required to calculate wet deposition, such as precipitation occurrence and rate as well as cloud water content. These uncertainties contribute to the overall uncertainty in the deposition parameterization (Gudiksen et al., 1988), thereby leading to inconsistencies in the wet deposition rates, which are more pronounced than the dry deposition rates. Despite the crucial role of precipitation in modeling, our mini-ensemble members exhibited poor performance in simulating precipitation values throughout the region, specifically in Qatar (Sect. 3.1). Moreover, FLEXPART and FLEXPART–WRF use different scavenging schemes (see Sect. 2.1). However, similar to the simulated concentrations, GFS-based simulations of ¹³⁷Cs (dry–wet) deposition show the closest agreement with FNL-based simulations. We also find that the largest difference between the simulated dry deposition of ¹³⁷Cs from all mini-ensemble members occurs between those forced by the FNL and ERA5–WRF datasets. Better agreement between the simulations based on the FNL–WRF and ERA5–WRF datasets than between those based on the FNL and FNL–WRF datasets is also obtained for the simulations of ¹³⁷Cs dry deposition, but not for the simulations of ¹³⁷Cs wet deposition. The above results also apply to the simulations of ¹³¹I deposition (Fig. S15), except that in this study ¹³¹I is assumed to be insoluble, remaining in the gas phase and not subject to wet deposition.

Figure 13 Same as Fig. 12 but for ¹³⁷Cs dry (in red circles) and wet (in blue squares) deposition (ng m⁻²).

4 Conclusions

In this study, we examined the spatiotemporal dispersion of radionuclides, including ¹³¹I and ¹³⁷Cs, in the fictitious event of nuclear accidents at the Barakah nuclear power plant (B-NPP) in the United Arab Emirates (UAE). The resulting concentrations and deposition of the studied radionuclides were simulated using the Lagrangian particle dispersion model FLEXible PARTicle (FLEXPART) and FLEXPART coupled with the Weather Research and Forecasting model (FLEXPART–WRF). To investigate the diurnal and seasonal variations of the radionuclide dispersion, the particles (air parcels) were released within the first 24 h of a 96 h simulation period at B-NPP, between 100 and 300 m above ground level, iterated daily over the year 2019. The source term is scaled to the maximum estimates of the radioactivity emissions from the Fukushima accident (6.9 kg or 22 PBq of ¹³⁷Cs and 0.042 kg or 192 PBq) for the purpose of providing a realistic benchmark for the reference values. We found differences in the simulations with respect to the meteorological inputs. We investigated the meteorological uncertainties by constructing a mini-ensemble with three members based on (re)analysis datasets including the final analysis (FNL) at native resolution, FNL and the ECMWF 5th Generation Reanalysis (ERA5) downscaled by WRF (ERA5–WRF and FNL–WRF), and one member based on the Global Forecast System (GFS) run by NCEP. We compared the daily means of surface wind speed and temperature as well as daily total precipitation from the above datasets with observations measured from 157 monitoring stations within the model domain. This mini-ensemble provided the basis for comparing the simulations based on the forecast and (re)analysis datasets. The FNL and GFS simulations were compared with the simulations based on the FNL–WRF and ERA5–WRF datasets to determine the impact of downscaling and of using different or the same model simulation codes (FLEXPART vs. FLEXPART–WRF) on the modeled dispersion. We conducted a sensitivity analysis on ERA5–WRF-based simulations, exploring the effects of variations in the turbulence scheme under convective conditions, emission duration, emission height, and emission location. We examined the simulations of all four members in relation to Qatar's population density, with the goal of identifying possible risks to populated areas from a nuclear accident in the region. A summary of the results of the study is presented below.

1. Transport of radionuclides. To analyze the time interval between emission and exposure to radionuclides, we investigated the age composition of the radionuclide plumes. Analysis of air parcel ages indicates that dense radionuclide clouds arrive in the south of the study area approximately 20 to 30 h after the emission. A significant amount of the ¹³¹I released is transported to the most distant parts of the study area up to 40 to 50 h after the accidents. All members simulated a large fraction of the ¹³⁷Cs deposition occurring within the first 75 to 80 h after the emission. The largest deposition of longer-lived particles was found in FNL-based simulations for all seasons, except in the fall. We attribute this to the more distant transport of air parcels from the emission point in FNL-based simulations compared to other members. The two members which are forced by the downscaled datasets (FNL–WRF and ERA5–WRF) simulated a more similar distribution of air parcel ages than those simulated by the other two members.

2. Distribution of extremely high concentrations and deposition of radionuclides. We calculated the seasonal median of 96 h integrated ¹³¹I concentrations (¹³¹I^{intg_conc_seas}, in units of ng m⁻³) and the total ¹³⁷Cs deposition (¹³⁷Cs^{tot_depos_seas}) across the study area. As expected, all members simulated much higher ¹³¹I^{intg_conc_seas} over the south-southeast of Qatar, which is the closest point of the study area to the emission point. The inter-seasonal comparison of the simulations showed that the mini-ensemble members simulated the largest advance of ¹³¹I^{intg_conc_seas} from the source to the south-southeast of Qatar in the cold period of the year. The simulations of ¹³¹I^{intg_conc_seas} between the members of the mini-ensemble differed by a factor of 10. The differences in the ¹³¹I^{intg_conc_seas} simulations based on FNL and FNL–WRF datasets demonstrated how the use of different model simulation codes and the downscaling of meteorological inputs can affect the FLEXPART modeling and, consequently, the decisions made based on its simulations after nuclear accidents. Similarly, remarkable differences were found in the spatiotemporal distribution of ¹³¹I^{intg_conc_seas} simulations based on the FNL and GFS datasets. This is the case even though these datasets are produced by the same base meteorological model. We attribute this to the fact that differences, however small, in the meteorological inputs that lead to cumulative deviations in the transport and concentration calculations of atmospheric pollutants. As with ¹³¹I^{intg_conc_seas}, ¹³⁷Cs^{tot_depos_seas} simulations from all the mini-ensemble members peaked in the southeastern part of the study area in the cold period of the year when both wet and dry deposition occur. The largest expansion of ¹³⁷Cs^{tot_depos_seas} above 10 ng m⁻² was found in the simulations based on the ERA5–WRF and FNL–WRF datasets in winter. The highest levels of ¹³⁷Cs^{tot_depos_seas} (above 30 ng m⁻²) were found in the simulations based on the ERA5–WRF dataset in the southeastern corner of Qatar in the fall. This region received far less ¹³⁷Cs^{tot_depos_seas} (around 10 ng m⁻² and less) in the simulations based on the GFS and FNL datasets in the same period. The examination of the release time of the air parcels resulting in extreme ¹³¹I concentrations and ¹³⁷Cs^{tot_depos_seas} showed that the corresponding particles are mostly released between 09:00 and 14:00 LT. The development of the boundary layer height, the intensification of the thermal gradient between the land and sea, and the resulting onshore winds increase the transport of radionuclides to the study area during this time of day. The analysis of the frequency with which ¹³¹I and ¹³⁷Cs concentrations above the 66th percentile are transported to the populated areas of eastern Qatar showed that, for all members, more than 50 % of the extreme cases occur in winter and between 15 % and 30 % in spring. The above results indicate that any nuclear accident in the winter will more likely be accompanied by the highest radionuclide concentrations and deposition within the study area. This pronounced intra-annual variation is attributed to a seasonal atmospheric pattern in which south-southeast winds transport dense radionuclide clouds. The collocation of population density showed that the populated areas (with more than eight persons per arc second) receive moderate to low levels of ¹³¹I^{intg_conc_seas} and ¹³⁷Cs^{tot_depos_seas} (below $\mathrm{0.4}\times {\mathrm{10}}^{-\mathrm{2}}$ ng m⁻³ and 40 ng m⁻²). Uninhabited areas in southern Qatar receive the highest levels of ¹³¹I^{intg_conc_seas} and ¹³⁷Cs^{tot_depos_seas} (above $\mathrm{0.6}\times {\mathrm{10}}^{-\mathrm{2}}$ ng m⁻³ and 60 ng m⁻²).

3. Sensitivity of simulations to model parameters and emission characteristics.

   We analyzed the sensitivity of radionuclide simulations based on the ERA5–WRF dataset to four variations of model setup. Results show that the use of the skewed turbulence model (STM) leads to a more frequent occurrence of high ¹³¹I concentrations and ¹³⁷Cs deposition within Qatar. Increasing emission height leads to a significant decrease in ¹³¹I concentrations and ¹³⁷Cs deposition, particularly during the cold period. Extending the emission duration causes a decrease in ¹³¹I concentrations and the amount of ¹³⁷Cs deposition. As anticipated, due to the longer distance, the change in the emission point from B-NPP to Bushehr NPP (Bu-NPP) results in reduced transfer of radioactive materials to Qatar, except during the spring season. The simulations from Bu-NPP exhibit distinct spatial patterns compared to those from B-NPP. Concentrations of ¹³¹I peak in the northern and eastern regions during winter and spring, which can be attributed to a southward shift of the westerlies. The deposition of ¹³⁷Cs follows a similar pattern. The entry of particles inducing high intensities of ¹³¹I and ¹³⁷Cs over the northern region of Qatar occurs within approximately 40–50 h and 60 after their release, respectively. The significant changes observed in the simulations during the sensitivity analysis can be attributed to the promotion of more downdrafts than updrafts with the STM, particle entry into the upper layers of the atmosphere with increasing emission height, less time for particles to reach the study area with extended emission duration, and the greater distance between the emission source and catchment area along with the varying impact of atmospheric patterns on the transport of radionuclides with a change in emission point.

4. Intercomparison of mini-ensemble members. The simulations of ¹³⁷Cs and ¹³¹I concentrations based on the GFS dataset were found to have the closest agreement with the simulations based on the FNL dataset because they share a common meteorological base model. The comparison of wind speed, precipitation, and temperature from the GFS and FNL datasets showed the best agreement compared to the other two input datasets. However, we also found important differences in the spatiotemporal distribution of GFS- and FNL-based simulations. This may be due to the cumulative effect of differences in the meteorological inputs on particle dispersion. The comparison of simulations based on ERA5–WRF and FNL–WRF indicated that the downscaling of the inputs and the application of the same simulation increase the agreement of resulting simulations to an extent that exceeds the degree of similarity between simulations based on FNL and FNL–WRF, with the same source of meteorological inputs. The deposition simulations of all members showed relatively large inconsistency for both radionuclides. This was more pronounced for the simulations of ¹³⁷Cs wet deposition. This is in part because wet deposition forcing factors are among the most challenging meteorological parameters to model accurately. Moreover, the recently updated wet deposition scheme implemented in the FLEXPART uses different methods to determine the occurrence of wet deposition than FLEXPART–WRF (Girard et al., 2016; Gudiksen et al., 1988; Evangeliou et al., 2017).