APERTRACK: A particle-tracking model to simulate radionuclide transport in the Arabian/Persian Gulf

Progress in Nuclear Energy 142 (2021) 103998 0149-1970/© 2021 The Author. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Contents lists available at ScienceDirect Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene APERTRACK: A particle-tracking model to simulate radionuclide transport in the Arabian/Persian Gulf R. Periáñez Dpt Física Aplicada I, ETSIA, Universidad de Sevilla, Ctra Utrera km 1, 41013 Sevilla, Spain ARTICLE INFO Keywords: Arabian/Persian Gulf Tide Baroclinic circulation Radionuclide Transport Sediment ABSTRACT A Lagrangian rapid-response model for simulating the transport of radionuclides in the Arabian (or Persian) Gulf is described. The model is based on a tide model including five constituents, which was solved in advance, and baroclinic circulation was obtained from HYCOM operational ocean model. The radionuclide model includes physical transport (advection and diffusion), radioactive decay and geochemical processes (interactions of radionuclides between water and sediments, described in a dynamic way). The model can lead with instantaneous or continuous releases. Some hypothetical releases from a coastal nuclear power plant were simulated. Results show that the moment of release affects the fate of radionuclides due to the temporal variability of baroclinic currents. Also, comparing results for releases of Cs and Pu, it was seen how the geochemical behaviour of the radionuclide clearly affects the further radionuclide distributions. It is easy to setup the model for a particular release and it provides a fast response; thus the present model is an appropriate tool to support decision-making after a nuclear accident. 1. Introduction The Arabian (or Persian) Gulf, from now on APG, is a shallow water body with a mean depth of 36 m (Alosairi and Pokavanich,2017). It is connected to the Gulf of Oman (Indian Ocean) through the Strait of Hormuz, thus it is a semi-enclosed marginal sea (Fig. 1). Countries which surround the APG are the United Arab Emirates, Saudi Arabia, Qatar, Bahrain (which consists of more than 30 islands in the APG), Kuwait, Iraq and Iran (this last in the eastern side). Circulation in the APG is forced by both winds and thermohaline (density driven) forcing. Given the excess of evaporation over precipitation and river inflow, a inverse estuarine circulation results; with the high salinity waters leaving the APG through a deep layer of the Strait of Hormuz and being replaced by a fresher surface inflow from the Indian Ocean (Kämpf and Sadrinasab,2006). This inflow occurs along the Iranian coast (Johns et al.,2003). Tides in the Gulf form standing waves, being dominant the semidiurnal and diurnal tides. The dimensions of the Gulf lead to a resonance of both tides, with one amphidromic point in the case of the diurnal and two in the case of the semidiurnal ones (Hyder et al.,2013). Desalination plants are the main freshwater source to the APG countries (Alosairi and Pokavanich,2017). For instance, in Abu Dhabi in 2007 desalination plants produced more than 2.3 million cubic metre of fresh water per day, which accounted for 36% of the total water production (Environmental Agency Abu Dhabi,2009) in such country. E-mail address: [email protected]. In addition, commercial and subsistence fisheries provide a living for a large sector of the coastal population (Abdi et al.,2006). The coastal environment of the APG has been exposed to various sources of radioactive pollution (Al-Ghamdi et al.,2016), including desalination plants (which are the main source of radium in the brine discharged to the sea) and phosphate industry (radium in phosphogypsum waste). Oil spills are relatively common in the APG, in addition to the massive oil releases during the 1991 Gulf War, which have both added natural radionuclides into the local marine environment. A review on radioactivity levels in the APG may be seen in Uddin et al. (2020). In addition to what it is commented above, the APG, Strait of Hormuz and Gulf of Oman are one of the most important waterways in the world, thus exposed to pollution incidents due to shipping activities (mainly potential oil spills). But recently, there has been concern about the nuclear power plants which are now operating along the APG coasts (Kamyab et al.,2018). There are two operational NPPs in the region, Bushehr in Iran and Barakah in UAE, whose unit 1 was connected to the power grid in summer 2020. About seventeen more are planned in the Kingdom of Saudi Arabia, with the intention that they are operational by 2030 (Uddin et al.,2020). Consequently, it is relevant to have a numerical model able to assess the effects of radioactive releases into the APG from such NPPs (or from a ship transporting nuclear wastes for instance). Discharges could https://doi.org/10.1016/j.pnucene.2021.103998 Received 17 February 2021; Received in revised form 5 October 2021; Accepted 8 October 2021 Progress in Nuclear Energy 142 (2021) 103998 2 R. Periáñez Table 1 Model availability. Program name APERTRACK Developer R. Periáñez, University of Sevilla Contact [email protected] Hardware Desktop PC Program code Fortran Cost Free Availability https://personal.us.es/rperianez/ be due to the normal operation of the plants or to acute accidental releases. A significant conclusion from IAEA (International Atomic Energy Agency) MODARIA and MODARIA-II (Modelling and Data for Radiological Impact Assessments) programmes (Periáñez et al.,2019a, 2016a;IAEA,2019) was the need to have site specific models which are carefully adapted to the region and made available for any marine area potentially exposed to a radionuclide release. This would help the decision-making process after an accident. Recent studies describing marine radionuclide transport models applied to other areas potentially exposed to nuclear accidents are, for instance, those of Periáñez et al. (2021) for the northern Indian Ocean and Tsabaris et al. (2021) for the eastern Mediterranean Sea. A review of models applied to simulate Fukushima releases in the Pacific Ocean may be seen in Periáñez et al. (2019a). Some models are described in literature concerning the dispersion of oil spills in the APG (Proctor et al.,1994;Faghihifard and Badri,2016; Al-Rabeh et al.,2000); however this is not the case with radionuclides. A radionuclide transport modelling work for the APG which could be found is that of Kamyab et al. (2018). These authors applied CROM1 model to simulate a hypothetical accident at Bushehr NPP; but CROM is essentially a Gaussian model based on the generic models described in IAEA (2001) suitable for steady conditions at a local scale, not able to deal with spatio-temporal variations of currents due to tidal oscillations and thermohaline forcing, thus its applicability in this case is questionable. Hassanvand and Mirnejad (2019) calculate tides in the northern APG and describe their effects in transporting radionuclides released from Bushehr in a qualitative way (without applying a transport model). They again use CROM to estimate transport and doses. The purpose of this paper is to fill such gap, presenting a radionuclide transport model for the APG which could be used for both chronic and accidental releases, including realistic descriptions of tidal and baroclinic currents, and finally including interactions of radionuclides between water and sediments; in line with recommendations in IAEA (2019). Moreover, the model is able to provide a fast response, thus it would be useful to support the decision-making process after an accident. Availability of the model is summarized in Table 1. The model is described in Section 2, where hydrodynamic methods (for tides and baroclinic circulation) and radionuclide transport description are presented separately. Results are presented in Section 3; first results of the tidal and baroclinic models are described (Section 3.1). Next some examples of simulations of radionuclide releases in the APG are presented (Section 3.2). 2. Model description 2.1. Tidal modelling A two dimensional depth-averaged model was used to simulate tides in the APG. Calculated elevations and currents are treated through standard tidal analysis (Pugh,1987, Chapter 4) and tidal constants (amplitudes and phases) are then calculated and stored for each grid cell in the computational domain. Five constituents were considered: three semidiurnal (𝑀2,𝑆2and 𝑁2) and two diurnal (𝐾1and 𝑂1). Tidal model 1ftp://ftp.ciemat.es/pub/CROM. equations (see for instance Periáñez,2012; a summary is presented in Appendix A) are solved for each constituent and tidal analysis is carried out for each constituent as well. The Eulerian residual transport is calculated, according to the procedure described in Periáñez (2012) and summarized in Appendix A, to obtain tidal residual currents. Boundary conditions to solve the equations consist of specifying water surface elevations and phases, from measured tidal constants, along the open boundaries of the domain. Measurements were obtained from Pous et al. (2012). The model domain extends from 47◦E to 57◦E in longitude and from 23◦N to 31◦N in latitude (Fig. 1). Resolution is the same as HYCOM model, 0.08◦(see Section 2.2). Once that amplitudes and phases (adapted phase, i.e., for the local time meridian) for each grid cell and constituent (calculated from the tidal analysis) are known, the tidal prediction equation is used to evaluate the exact tidal state during each time step of the radionuclide simulation and location in the APG. The procedure is described in Parker (2007) and Boon (2011) and summarized in Appendix A. The tidal model is two-dimensional, thus it provides averaged currents over the water column. A three-dimensional current field is generated using a standard current profile, since currents decrease from sea surface to the bottom because of friction. Details may be seen in Pugh (1987) and Periáñez and Pascual-Granged (2008). The present tidal model was successfully tested for several regions at quite different spatial scales (Periáñez,2007,2009,2012;Periáñez et al.,2013;Periáñez and Abril,2014;Periáñez,2020a). 2.2. Baroclinic circulation HYCOM (Hybrid Coordinate Ocean Model, (Bleck,2001)) model was used to obtain baroclinic circulation in the APG. HYCOM is a primitive equation general circulation model with 40 vertical layers increasing in thickness from the surface to the sea bottom and 0.08◦ horizontal resolution in both latitude and longitude. Examples of HYCOM model applications over the world are presented in the model web page (https://www.hycom.org/). Actually, this model has already been used to study circulation in the APG (Yao and Johns,2010a,b). Daily currents were downloaded from HYCOM data server for the APG (the same domain specified above for the tidal model). Note that the tidal model is required since tides are not included in HYCOM. 2.3. Radionuclide transport The model is Lagrangian as commented before, thus the radionuclide release into the sea is simulated by means of a number of particles. Each particle is equivalent to a number of units (for instance Bq), and trajectories are calculated during the simulated period. The transport model considers physical transport (advection due to water currents and mixing due to turbulence) plus radioactive decay and interactions of radionuclides with bed sediments (adsorption/desorption reactions). Radionuclide concentrations are obtained from the number of particles within each grid cell and compartment (surface water, deep water and sediment as explained in Appendix B) and the number of units (Bq) which corresponds to each particle. Turbulent mixing, radioactive decay and exchanges of radionuclides between water and sediment are described through a stochastic method (Periáñez and Elliott,2002;Kobayashi et al.,2007;Periáñez et al., 2019a). A dynamic method is applied to describe water/sediment interactions, thus a kinetic coefficient 𝑘1describes the transfer of radionuclides from water to sediment and a coefficient 𝑘2governs the inverse process. A summary of the involved equations may be seen in Appendix B. As in other works, 𝑘1is derived from the radionuclide equilibrium distribution coefficient 𝑘𝑑(provided for instance in (IAEA, 2004)) and a standard experimental value for 𝑘2(Periáñez,2009; Periáñez et al.,2013,2016b). Equations are summarized in Appendix B. Progress in Nuclear Energy 142 (2021) 103998 3 R. Periáñez Fig. 1. Map of the APG, which corresponds to the present model domain. Isobaths of 20, 40, 60, 80 and 100 m are drawn. The locations of Bushehr and Barakah NPPs are also shown. Fig. 2. General scheme of the modelling procedure. The user must specify only release data and other. Equilibrium arguments and nodal factors for year 2021 are set as default option. 2.4. Model input A number of files specify the release characteristics (date, time, position in geographic coordinates, depth, magnitude and duration) and simulation time, radionuclide properties (decay constant and equilibrium distribution coefficient (which may be obtained from IAEA (2004) as mentioned above), and, finally, an optional wind forecast (see next paragraph) and components of the currents to be used: tidal currents and residuals may be individually switched on and off (to Progress in Nuclear Energy 142 (2021) 103998 4 R. Periáñez Fig. 3. Comparison between calculated and observed amplitudes for the semidiurnal (left) and diurnal (right) constituents considered in the model. The map in the top shows points where tidal constants were measured (black dots). allow comparisons if they are included or not in the simulations or to speed them up by removing tides in the calculations). These switches are provided in a specific file named input.dat. The file tidedata.dat contains equilibrium arguments and nodal factors of the 5 tidal constituents for year 2021, set as default, as explained in Appendix A. Thus, this file should be modified only if a simulation for a different year is to be carried out. Equilibrium arguments and nodal factors for the corresponding year should then be used. A list of the input files which should be modified for a particular simulation is given in Table 2. In the case of a simulation to assess the effects of an acute release due to an accident, for instance, it may be relevant to include a local wind, which is considered uniform in the release area. Wind data are provided in a file as a number of different ‘‘wind episodes’’ (any number can be used with a maximum of 100), each one characterized by a wind speed, direction and start and end times measured in hours after the pollutant release beginning. This time-evolving wind conditions may be obtained from weather forecasts. It should be commented that HYCOM calculations already include atmospheric forcing. However, the present Table 2 Input files which must be modified for each specific simulation. It is required to modify tide-data.dat only if a simulation for other year than 2021 is to be carried out. tide-data.dat Equilibrium arguments and nodal factors release.dat Release data and simulation time RN.dat Contaminant properties (decay constant and 𝑘𝑑) input.dat Switches to include or not tidal circulation wind.dat Local wind data definition of ‘‘wind episodes’’ gives the opportunity of describing transport in case that an accident occurs, for instance, during a local storm which is not described in HYCOM. The need of adding this local wind in some oil spill simulations in the Red Sea was clearly shown in Periáñez (2020a) and was also used in a radionuclide transport model for the same sea (Periáñez,2020b). The wind-induced current is considered to decrease logarithmically to zero from the surface. The mathematical form of this profile may be seen in Pugh (1987), for instance. It should be clearly pointed out that using this ‘‘local wind’’ is Progress in Nuclear Energy 142 (2021) 103998 5 R. Periáñez Fig. 4. Calculated chart for the 𝑀2tide. Phases are given with respect to the local time meridian (adapted phases). Fig. 5. Calculated chart for the 𝑂1tide. Phases are given with respect to the local time meridian (adapted phases). optional and should be included only if a wind forecast is known and it includes unusual weather conditions. Otherwise atmospheric forcing already included in HYCOM calculations is enough for the transport calculations. 2.5. Model output The model output consists of radionuclide concentrations over the model domain in two water layers: a surface layer whose thickness is defined as 10 m, but can be changed by the user in the code, and a deep layer which extends from the bottom of the surface layer to the seabed. Actually, the model provides the radionuclide inventory in units/m2in the deep layer. Concentrations in bed sediments are provided in a 5 cm thick sediment layer. In addition, the model provides the position of particles (both in the water column and in sediments) at the end of the simulation. All this information may be drawn with the Octave scripts which are provided with the model. A general scheme of the modelling procedure is presented in Fig. 2. All required inputs are in blue boxes. The marine data is pre-computed and does not require any action by the user, which only needs to modify the release data and other. Once input is defined, the transport code (pink) performs the calculations and provides output (green). The number of particles used in the model is 200000. A simulation over three months takes about 10 min on a desktop PC working over Ubuntu 18.04 operating system. All the required codes were written in Fortran. Progress in Nuclear Energy 142 (2021) 103998 6 R. Periáñez Fig. 6. Water circulation as downloaded from HYCOM model at the end of four months of the year for the sea surface. Only one of each 16 vectors is drawn for more clarity. 3. Results 3.1. Hydrodynamics The tidal model was calibrated changing the bed friction coefficient until the best agreement between calculated and observed (from Pous et al. (2012)) tidal elevations was achieved. Such agreement was measured as 𝜒2, according to the equation (Glover et al.,2011): 𝜒2=1 𝑁 𝑁 ∑ 𝑖=1 (𝑍𝑜𝑏𝑠 𝑖−𝑍𝑐𝑎𝑙 𝑖)2 𝜎2 𝑖 (1) where 𝑁is the number of observations, 𝑍𝑜𝑏𝑠 𝑖and 𝑍𝑐𝑎𝑙 𝑖are observed and calculated elevations respectively, and finally 𝜎𝑖is the uncertainty in each measurement, taken as 0.01 m according to the observations presented in Pous et al. (2012). A comparison between observed and calculated elevations for the five constituents may be seen in Fig. 3. Although agreement is generally good, there are stations where higher discrepancies appear. Most likely it is due to the relatively coarse resolution of the model: tides where simulated using the same grid as HYCOM, which is 0.08◦. Using a finer grid would improve results, but this would be overcome by errors and difficulties in interpolating currents from one grid to the other in order to deal simultaneously with tidal and baroclinic currents. As a couple of examples, tidal charts for one semidiurnal (𝑀2) and one diurnal (𝑂1) tide are respectively presented in Figs. 4 and 5. These charts are in good agreement with earlier calculations made for the APG (Pous et al.,2012;Hyder et al.,2013;Akbari et al.,2016). Thus, in the case of the 𝑀2tide there are two amphidromes, at (50◦E, 28◦N) and (53◦E, 25◦N) approximately. In contrast, diurnal tides show a single one. In case of the 𝑂1tide it is located approximately at (52◦E, 27◦N). These locations are in agreement with those presented in Akbari et al. (2016). Fig. 6 presents a few examples of surface water circulation as calculated by HYCOM model at the end of the indicated months. Circulation is essentially cyclonic in January and anticyclonic in September, showing the well-known surface inflow of Indian Ocean waters along the Iranian coast (Johns et al.,2003). A cyclonic eddy is also apparent in the northern part in July. Actually, it was found, through numerical simulations, that in summer the north-westward coastal current flowing along Iran evolves into a series of mesoscale anticyclonic eddies with typical diameter about 120 km. One of these eddies was apparent in the northern region of the Gulf (Thoppil and Hogan,2010). 3.2. Radionuclide dispersion The model can be applied to any radionuclide, simply using its specific distribution coefficient and radioactive decay constant. Here we present some examples with 137Cs and 239,240Pu, which have very different geochemical behaviours: the first is quite conservative while plutonium presents a high affinity to be fixed to sediment particles. A summary of model runs which were carried out is presented in Table 3. An hypothetical accident occurring at Bushehr NPP (coordinates 50.88◦E, 28.82◦N) was simulated. A 137Cs release was supposed to last 90 days, with a total activity released equal to 1 PBq. This is just an example, but it is the same order of magnitude as the direct release from Fukushima into the Pacific Ocean during the first three months after the Progress in Nuclear Energy 142 (2021) 103998 7 R. Periáñez Fig. 7. 137Cs concentrations (logarithmic scale) in surface water (Bq/m3) and bed sediments (Bq/kg) after 90 days of a release starting in March 21 in Bushehr NPP (run 1). Details of the hypothetical accident are given in the text. Table 3 Summary of model runs. Starting time of the releases was 12:00 h local time in all cases (year 2021). Local wind was not included in any case and all tidal constituents and residuals were considered. Release magnitude was 1 PBq during 90 days in all runs. Radionuclide Location Simulated time Starting time Run 1 137Cs Bushehr NPP 90 days March 21 Run 2 137Cs Bushehr NPP 90 days June 21 Run 3 137Cs Bushehr NPP 90 days September 21 Run 4 137Cs Bushehr NPP 90 days December 21 Run 5 137Cs Bushehr NPP 1 year March 21 Run 6 239,240Pu Bushehr NPP 90 days March 21 Run 7 137Cs Barakah NPP 90 days March 21 2011 tsunami (Kobayashi et al.,2013). The 137Cs 𝑘𝑑was fixed as 4.0 m3/kg, which is the established value for coastal waters by IAEA (2004) and radioactive decay constant for this radionuclide is 7.29 × 10−10 Fig. 8. 137Cs concentrations (logarithmic scale) in surface water (Bq/m3) and bed sediments (Bq/kg) after 90 days of a release starting in June 21 (run 2). Details of the hypothetical accident are given in the text. s−1(half life of 30.17 year). The release was supposed to occur at the sea surface and simulation time was 90 days. Four simulations were carried out with different starting times: 21 March, 21 June, 21 September and 21 December. All releases were finally supported to start at 12:00 h local time. The optional local winds are not included in these calculations, but only the atmospheric forcing already described within HYCOM model. The five tidal constituents and their residuals were included (all switches in file input.dat set to 1). The simulations shown as examples are relatively long (90 days) simply to illustrate general transport patterns in the APG. In the case of an accident it may be relevant to carry out short term (few days) simulations to support decision-making and undertake preventing actions in the region around the accident. Maps of 137Cs in surface water, taken as a 10 m thick layer, and bed sediments after the simulations were obtained, which are presented in Figs. 7 to 10. It seems clear that the starting time of the release Progress in Nuclear Energy 142 (2021) 103998 8 R. Periáñez Fig. 9. 137Cs concentrations (logarithmic scale) in surface water (Bq/m3) and bed sediments (Bq/kg) after 90 days of a release starting in September 21 (run 3). Details of the hypothetical accident are given in the text. affects the subsequent radionuclide distributions due to the temporal variability of baroclinic circulation. Thus, if the release starts with spring (Fig. 7) radionuclides move to the north and to the central APG, then travelling to the south along the western side. Sediments are contaminated as waters containing 137Cs move over them. Since the water/sediment interaction model is dynamic, sediments buffer radionuclides which are later released as water above them is cleaned. Thus, the concentration map for surface water is an instantaneous picture of the radionuclide distribution at exactly that time; but the map for sediments integrate the whole path followed by the release. If the release starts with summer (Fig. 8), the sediment map indicates that transport has been predominantly directed to the north, while it is directed to the south if the release starts with fall (Fig. 9). In this case there is also some transport to the south along the western side, as in Fig. 7. Finally, if the release starts with winter (Fig. 10) radionuclides remain close to Bushehr NPP; transport is mainly directed to the south Fig. 10. 137Cs concentrations (logarithmic scale) in surface water (Bq/m3) and bed sediments (Bq/kg) after 90 days of a release starting in December 21 (run 4). Details of the hypothetical accident are given in the text. along the Iranian coast and radionuclides do not reach the western coast of the APG in the simulated temporal frame. As a conclusion, it seems evident that the moment when an accident occurs determines the fate of the released radionuclides and the portion of the APG coast which is potentially contaminated. However, the four simulations show that radionuclides do not reach the north extreme of the APG. It can be probably attributed to the freshwater input from Shatt Al Arab river (Tigris and Eufrates) at the Gulf head, although it should be noted that the present day inflow is much smaller than it once was because of dam projects in Turkey (Hyder et al.,2013). The accident starting in March (Fig. 7) has been simulated during one year and results are presented in Fig. 11; where 137Cs concentrations in surface water, bed sediment and inventory of radionuclides in the bottom water layer (from 10 m depth to the seabed), in Bq/m2, may be seen. If the simulation time is extended, a significant amount of radionuclides reach the bottom water layer and are able to contaminate Progress in Nuclear Energy 142 (2021) 103998 9 R. Periáñez Fig. 11. Same as Fig. 7 but for a one year long simulation (run 5). Inventory (Bq/m2) of 137Cs in the bottom water layer is also shown. Note the different colour scales for waters. Details of the hypothetical accident are given in the text. the bed sediments. Actually, virtually all the sediments of the APG contain 137Cs (Fig. 11). Radionuclides in the bottom water layer reach the Strait of Hormuz, travelling with the deep outflow water, and will leave the APG entering the Gulf of Oman. The geochemical behaviour of the radionuclide affects the fate of the release. For instance, the experiment shown in Fig. 7 was repeated Fig. 12. 239,240Pu concentrations (logarithmic scale) in surface water (Bq/m3) and bed sediments (Bq/kg) after 90 days of a release starting in March 21 (run 6). Details of the hypothetical accident are given in the text. but supposing that the released radionuclide was 239,240Pu, whose recommended 𝑘𝑑value is 100 m3/kg according to IAEA (2004). Thus, it is much more reactive than 137Cs, presenting a higher affinity to be fixed to the sediment. This can be clearly seen comparing Fig. 12, which shows the plutonium results, with the previous Fig. 7:239,240Pu is quickly fixed to the sediments in the release area, thus presents low mobility in a shallow marine environment like the APG is. As a final example, exactly the same accident as shown in Fig. 7 was simulated for 137Cs but occurring in Barakah NPP (coordinates 52.23◦ E, 23.97◦N) in UAE. Thus, details on the release are presented above. Concentrations resulting from this Barakah NPP release can be seen in Fig. 13. In this case the released 137Cs moves towards the Strait of Hormuz, but currents in this region of the APG are weaker and the extension of the contaminated area is much smaller than for the previous simulation.