Models for predicting the transport of radionuclides in the Red Sea

Models for predicting the transport of radionuclides in the Red Sea R. Periáñez Dpt Física Aplicada I, ETSIA, Universidad de Sevilla, Ctra Utrera km 1, 41013 Sevilla, Spain A B S T R A C T Two radionuclide transport models for the Red Sea are described: a Lagrangian model to deal with emergency situations and a Eulerian model better suited to longer term simulations, as for instance required in case of chronic radionuclide releases. Baroclinic circulation is obtained for both transport models from HYCOM ocean model. The Lagrangian model also includes tides, which are obtained from a standard tidal model customized to the Red Sea, and currents induced by local winds. Both models describe exchanges of radionuclides between water and sediments. A number of simulations were carried out to illustrate capabilities of the models. Additionally, flushing times over the Red Sea were evaluated with the Eulerian model, as another example of model use. 1. Introduction The Red Sea suffers intense shipping activities since it is located between the Indian Ocean and the Mediterranean Sea, connected to the last through the Canal of Suez. Bab el-Mandeb (BeM) Strait connects it with the Indian Ocean (Fig. 1). Several authors (Mamoney and Khater,2004;El-Taher et al.,2018a,b) have recognized the potential radioactive impact along the Egyptian Red Sea coast due to phosphate and oil industries, as well as due to ship traffic. Some of these ships are nuclear powered or may be carrying radioactive materials. Thus they are sources of potential accidental radioactive contamination. This contamination risk increased during the first and second Gulf wars (Saharty and Dar,2010), when the Red Sea was used to transfer military ships, weapons and submarines from the Mediterranean Sea to the Arabian Gulf. In present days, ships carrying nuclear wastes sail along the Red Sea (Saharty and Dar,2010). The IAEA (International Atomic Energy Agency) has launched MODARIA and MODARIA-II (Modeling and Data for Radiological Impact Assessments) programs in the last years with the purpose, among others, of improving numerical tools to support decision-making after an accidental or deliberate radionuclide release in the sea. A significant conclusion from these programs (Periáñez et al.,2016a,2019a) was the need of having site specific tools which are carefully developed, adapted to the site and tested; and finally made available for any marine area potentially exposed to a radionuclide release. A rapidresponse, as required by decision-makers, could then be achieved in case of an accident. The purpose of this work is to present a number of tools developed for the Red Sea (for which radionuclide transport models cannot be found in literature); which illustrates the procedure by which the needs described in the IAEA programs could be satisfied. As described in Periáñez et al. (2016a,2019a), three stages after a nuclear accident may be defined: emergency phase, post-emergency phase and long term phase. These phases are characterized by increasing spatio-temporal scales, and each one requires a specific kind of model to give response to decision makers. Two radionuclide transport models are consequently described in this work: 1. Lagrangian model: It is useful to make predictions in the initial stages of the accident. 2. Eulerian model: It is used to make predictions at longer temporal scales and particularly useful in the case of releases which are persistent in time, since the number of particles required in a Lagrangian model would be very large in such situations. A detailed description of the relative advantages of Eulerian and Lagrangian models is presented in Periáñez et al. (2019a) and it is not repeated here. Both models are fully three-dimensional, due to the stratified nature of the Red Sea, and may be used both for conservative (remain dissolved) and for reactive radionuclides. Consequently, the contamination of the bed sediments can be predicted as well. The Eulerian model is forced by baroclinic circulation in the sea (due to density differences and atmospheric circulation). The Lagrangian model, in addition, includes tidal currents and local wind forecasts, in order to have more accurate predictions in the area near the accident. Actually, it was found (Periáñez,2020) that baroclinic currents are the dominant transport forcing within most of the Red Sea, both for typical winter and summer conditions; tides and tidal residuals being essentially negligible in most of the sea (tides are only significant in the region of BeM Strait). Thus, we provide the opportunity of using more detailed circulation fields (including tidal effects) in the case of an accident leading to an acute release (a Lagrangian model is more E-mail address: [email protected]. 2 Fig. 1. Model domain. Water depths are measured in m below the mean sea level. appropriate in this case). But the use of only baroclinic circulation may be enough in most cases (see details in Periáñez,2020). However, of course it is possible to force the Eulerian model with tides and tidal residuals as well. The models are based upon previous ones developed by the author, thus they are only briefly described in Section 2, where the hydrodynamic and radionuclide transport components are presented in separate paragraphs and appropriate references are provided. Results are described in Section 3. Results from the transport models are given in Section 3.2, where some application examples are described. Finally, the Eulerian model was applied to a conservative radionuclide to evaluate e-folding times over several regions of the Red Sea and in different seasons (Section 3.2.4). 2. Methods 2.1. Hydrodynamic methods Water currents in the sea result from the addition of tidal currents, baroclinic circulation and local wind conditions. The procedures to describe these currents are presented in this section. Tides were computed by means of a two dimensional depthaveraged model. Elevations and currents predicted by the model were treated through standard tidal analysis (Pugh,1987) and tidal constants (amplitudes and phases) were then calculated and stored. Five constituents were considered: three semidiurnal (𝑀2,𝑆2and 𝑁2) and two diurnal (𝐾1and 𝑂1). A residual transport cannot be produced by the pure harmonic currents which result from tidal analysis. Consequently, the tidal model also calculated the Eulerian residual transport for each tidal constituent. Details of this procedure may be seen in Periáñez (2012) for instance. Since the tidal model is two-dimensional, it provides averaged currents over the water column. A standard current profile is then used to generate a vertical structure in tidal currents, which decrease from sea surface to the bottom because of friction with the seabed. A 3D tidal current field is produced in this way, which can used in the threedimensional transport model. Details may be seen in Pugh (1987) and Periáñez and Pascual-Granged (2008). The tidal model was successfully applied in the past to a number of different areas (Periáñez,2007,2009,2012;Periáñez and Pascual- Granged,2008;Periáñez and Abril,2014;Periáñez et al.,2013). Equations and details may be seen in these references. In addition, tides calculated in the Red Sea were already used in an oil-spill transport model developed for this region (Periáñez,2020). A comparison between measured and calculated tides for the five considered constituents is presented in such paper and not repeated here. Baroclinic circulation in the Red Sea was obtained from HYCOM (Hybrid Coordinate Ocean Model, Bleck,2001) model. It is a 3D primitive equation global circulation model with 40 vertical layers and its horizontal resolution is 0.08◦. A number of HYCOM model applications can be seen in the model web page (https://www.hycom.org/). Daily currents for the area of interest were obtained from HYCOM data server. Local winds can be obtained from weather forecasts. Following Proctor et al. (1994) the wind-induced current decreases, below the wind-driven surface layer, logarithmically to zero at a depth which typically is about 20 m. The mathematical form of this current profile may be seen in Pugh (1987), for instance. The surface current is a percentage of the wind speed, usually in the range 2%–3%. The winddriven layer speed was increased to 3.5% of the wind speed to account for Stokes drift (Proctor et al.,1994). Although the surface current should be deflected to the right/left (northern/southern hemisphere), this angle is small and may be neglected (Proctor et al.,1994). The procedure to include the local wind forecast in the transport model is presented in Section 2.4. 2.2. Lagrangian model A pollutant release into the sea is simulated in a Lagrangian model by means of a number of particles, each one equivalent to a number of units (for instance Bq), whose trajectories are calculated during the simulated period (Periáñez et al.,2016b). The transport model considers physical transport (advection due to water currents and mixing due to turbulence) plus some specific processes for radionuclides: radioactive decay and interactions with bed sediments. Turbulent mixing, radioactive decay and exchanges of radionuclides between water and sediment are described through a stochastic method whose details can be seen in a number of papers (Periáñez and Elliott, 2002;Kobayashi et al.,2007;Periáñez et al.,2019a). A summary of the main equations can be seen in Appendix A. In particular, a dynamic method is applied to describe water/sediment interactions. A kinetic coefficient 𝑘1describes the transfer of radionuclides from water to sediment and a coefficient 𝑘2governs the inverse process. As in other works, 𝑘1is derived from the radionuclide equilibrium distribution coefficient 𝑘𝑑(IAEA,2004) and a standard experimental value for 𝑘2 (Periáñez,2012;Periáñez et al.,2013,2016b). Typical values are used for the diffusion coefficients, which are considered to be constant: 1.0 × 10−5 m2/s and 10 m2/s for the vertical and horizontal coefficient respectively (Elliott et al.,2001;Periáñez et al., 2019b). As commented in Appendix A this is just an approximation to speed up calculations, but more complex descriptions are possible. The horizontal diffusion coefficient is appropriate for the present grid resolution, according to the standard relation (Periáñez,2005): 𝐾ℎ= 0.2055 × 10−3𝛥𝑥1.15 (1) This equation leads to a diffusion coefficient of 7,4 m2/s for the present resolution of HYCOM grid, which has been rounded to 10 m2/s. Time step was fixed as 600 s. 3 0 Fig. 2. Baroclinic circulation in the Red Sea at the end of January and July. Only one of each 16 vectors is drawn for more clarity. Fig. 3. Concentrations in surface water (Bq/m3) five days after the release calculated with the Eulerian and Lagrangian models. The black dot is the release point. 2.3. Eulerian model The Eulerian model is based upon the numerical solution of the three dimensional advection/diffusion equation for the radionuclide concentrations over the model domain. Additional terms are added to this equation in order to account for radioactive decay and water/sediment interactions. These interactions are described through kinetic coefficients, exactly as in the Lagrangian model. Equations and details may be seen for instance in Periáñez (2009,2012) and Periáñez et al. (2013), and a summary is given in Appendix B. Second order accuracy finite difference schemes were used for the advection and diffusion terms. Time step was fixed as 600 s to satisfy stability conditions. The same diffusion coefficients as in the Lagrangian model were used in the Eulerian one. 2.4. Computational scheme The model domain extends from 10◦to 30.5◦N in latitude and from 32◦to 45◦E in latitude. Bathymetry was obtained from NOAA ETOPO1 data, using its on-line tool to design a grid1at a resolution of one minute of arc in both longitude and latitude. It is presented in Fig. 1. As mentioned before, all the hydrodynamic data is obtained in advance and stored in files which are read by the radionuclide transport codes. Additional files are required by these codes, which specify the release characteristics (date, time, position, depth, magnitude), radionuclide properties (radioactive decay constant and equilibrium 1https://maps.ngdc.noaa.gov/viewers/wcs-client/. 4 Fig. 4. Calculated concentrations in surface water (Bq/m3, left side) and bed sediments (Bg/kg, right side) after 60 and 182 days of continuous release at the black dot location. Concentrations are drawn in logarithmic scale. Table 1 Summary of required inputs for each model. The file input.dat is described in the text. Lagrangian Eulerian HYCOM currents Yes Yes Tide amplitudes and phases Yes No Tidal residuals Yes No File input.dat Yes No Equilibrium arguments and nodal factors Yes No Bathymetry Yes Yes Release data and simulation time Yes Yes Radionuclide properties Yes Yes Local wind data Yes No distribution coefficient — which may be obtained from IAEA (2004)), simulation time and, depending on the transport model which is used, wind forecast and components of the currents to be used, as explained in what follows: tidal currents and residuals may be individually switched on and off (to 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. A summary of the required input for each model is given in Table 1. Baroclinic currents were previously downloaded from the HYCOM website as already mentioned. Tidal amplitudes and phases are provided for the five included constituents, as well as Eulerian residual transports for each of them. The file input.dat provides the switches, as zeros or ones, for each constituent tidal current and residual in order that each one is included (1) or not (0) in the simulation as mentioned before. Equilibrium arguments (at Greenwich) and nodal factors for each constituent and for the year of simulation are parameters required to predict the exact tidal state at the day and time of the simulation (Pugh,1987;Parker,2007;Boon,2011). Wind conditions (speed and direction) are not uniform over the model domain. However, the local wind 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), each one characterized by a wind speed, direction and start and end of the episode in hours after the beginning of the radionuclide release. This time-evolving wind conditions may be obtained from appropriate weather forecasts. It must be noted that atmospheric forcing is already included in the HYCOM calculations of water circulation. However, with the present definition of ‘‘wind episodes’’, we have the opportunity of describing radionuclide transport in case that the accident occurs during a local storm, for instance, which is not described in HYCOM. The need of adding this local wind in some cases was clearly shown in Periáñez (2020). The models provide radionuclide concentrations over the model domain in two water layers: a surface layer whose thickness may be defined by the user and a deep layer which extends from the bottom of the surface layer to the seabed. In addition, concentrations in bed sediments are also provided. The thickness of the surface sediment layer over which radionuclide concentration is calculated is defined as 5 cm. 3. Results 3.1. Hydrodynamics As already mentioned, the tidal model was recently applied to the Red Sea (Periáñez,2020). Computed tidal charts and comparisons between calculated and measured tides may be seen in such paper and thus are not repeated here. 5 Baroclinic circulation for the end of January and July may be seen in Fig. 2 as a couple of examples. These winter and summer conditions are different: in winter the flow is mainly directed to the north along the east side of the sea. In summer this circulation is not apparent. Now the flow reverses and is directed to the south in BeM Strait. This effect is known since a long time (Patzert, 1974) and is attributed to monsoon reversal. 3.2. Radionuclide transport 3.2.1. Comparison of the Lagrangian and Eulerian models The relative advantages of Eulerian and Lagrangian models have been discussed in detail in the review paper by Periáñez et al. (2019a), and the stochastic method used to describe water/sediment interactions was compared with the Eulerian method in Periáñez and Elliott (2002). However, it is illustrative to present results for the same simulation with the two models described in this paper. The transport of a hypothetical instantaneous release of 1012 Bq of 137Cs occurring at the surface in the point with coordinates 40.5◦E, 18◦N was simulated. The release was supposed to occur on January 31 and a 5 day long simulation was carried out. The recommended (IAEA, 2004) 𝑘𝑑 for Cs in coastal waters (4.0 × 103) was used. Calculated concentrations in the sea surface are presented in Fig. 3. The Eulerian model is clearly more diffusive than the Lagrangian one, as already known, in spite that a second-order accuracy advective scheme is used. Although the same color bar is used in both maps, peak concentrations are 210 and 115 Bq/m3 in the Lagrangian and Eulerian models respectively. Water currents at this time of year and area of the sea are directed to the north (Fig. 2). The Lagrangian model gives an elongated radionuclide patch in the current direction, but the Eulerian one produces a more intense mixing in the direction transverse to the current. This is numerical diffusion since the same diffusion coefficients were used in both models. Given the spatial resolution of the Eulerian model (0.08◦) numerical diffusion is significant (it is proportional to the grid cell size and to the water velocity (Periáñez, 2005)). Actually, the same experiment as in Fig. 3 but with horizontal diffusion set to zero was repeated and results were virtually the same as with the used nominal diffusion equal to 10 m2/s. Thus, results of the Eulerian model cannot get closer to the Lagrangian model ones in regions with high currents. Only in regions with very weak currents numerical diffusion will tend to zero. Presented results confirm the fact that Lagrangian models are better suited to emergency situations since they can handle the high radionuclide concentration gradients between clean and contaminated water, which appear after an acute release in the sea, better than Eulerian models. Note that after 5 days, peaks concentration provided by the Eulerian model is a factor two lower than the Lagrangian model peak. An additional factor contributing to this effect, as discussed in Periáñez et al. (2019a), is the fact that a true point source can be defined in a Lagrangian model. In contrast, radionuclides are instantaneously mixed into the release grid cell in an Eulerian model. 3.2.2. A chronic release The Eulerian model was applied to simulate a hypothetical continuous release occurring in the surface at a coastal location (35.2◦E, 25.0◦N). It was supposed to start on April 30 and a half year long simulation was carried out with a constant radionuclide release rate equal to 106 Bq/s. The radionuclide was 137Cs. Calculated concentrations in surface water and sediments (5 cm thick surface layer) after 60 and 182 days are presented in Fig. 4 as examples, where such concentrations are drawn in logarithmic scale. It may be seen that currents transport radionuclides towards the south, although a branch of marked water is directed to the north and some radionuclides are present in an eddy in the northern Red Sea. Only sediments along the coast are contaminated; water depth increases Fig. 5. Calculated concentrations in surface water (Bq/m3, up) and bed sediments (Bg/kg, down) after 10 days of an instantaneous release at the black dot location. rather quickly from the coast, reaching more than 1000 m (Fig. 1) in the central Red Sea at these latitudes. Radionuclides do not reach such deep sediments, although the band of contaminated sediment is increasing width with time, as can be seen in the figure. 3.2.3. Acute releases Two hypothetical accidents have been simulated to illustrate the model performance, one occurring in shallow water near the coast, and the other in the central Red Sea, in deep water. The first accident was supposed to occur at position 40.5◦E, 15.2◦N, where water depth is only 41 m, on February 1st as an example. A total amount of 1012 Bq of 137Cs was instantaneously released at the surface, and a 10 day long simulation was carried out. Resulting concentrations in surface water, defined as a 20 m thick layer, and sediments (5 cm thick surface layer) are presented in Fig. 5. The radionuclide patch is transported to the north by water currents, passing between the many small islands in the area. Bed sediments are contaminated as water containing radionuclides moves over them. Since water/sediment interactions are described by means of a dynamic model, the shapes of the radionuclide distributions in water and sediments are not the same. Once the sediment is contaminated, it buffers radionuclides 6 Fig. 6. Calculated concentrations in surface water (Bq/m3, up) and bed sediments (Bg/kg, down) after 10 days of an instantaneous release at the black dot location. Release depth is 1 m. and consequently the patch of marked sediments extends from the release point with similar values of concentrations. In contrast, the radionuclide peak in water is well located at about 15.5◦N. The same accident was now supposed to occur at position 39.0◦E, 20.0◦N, where water depth is 619 m. Four simulations were carried out changing the release depth: 1 m, 100 m, 300 m and 600 m. The purposes were to illustrate the functioning of the model and also to show how the release depth affect the fate of radionuclides since water currents have a significant vertical structure. Some results are presented below. In the first case (release depth 1 m) the resulting radionuclide concentrations in the surface water layer, again defined as 20 m thick, and in bed sediments are presented in Fig. 6. Surface currents transport dissolved radionuclides towards the north, reaching shallower areas where bed sediments are contaminated. Actually two radionuclide patches in sediments are observed in Fig. 6 (down). Results are very different if releases are deeper. In all other experiments radionuclides do not reach the surface water layer. In order to summarize results, only the positions of particles in water (all of Fig. 7. Positions of particles in water and sediments after 10 days of an instantaneous release at the black dot location for different release depths. them in the deep layer) and bed sediments are presented in Fig. 7. Resulting concentrations are not shown. If the release occurs at 100 m, radionuclides travel to the west, reaching even deeper waters and consequently they do not reach the seabed. The magnitude of water 7 currents decreases with depth, thus when releases are deeper more radionuclides stay closer to the release point, as can be observed in the cases of 300 and 600 m. Some small spots of contaminated sediments can be also observed in these cases. Finally, the different shapes of the radionuclide patches indicate the water current variability with depth. 3.2.4. Flushing times As a final example, the Eulerian model was applied to estimate water flushing (or turn-over) times in the Red Sea, for surface water. These are defined (see for instance Prandle, 1984) as the time in which the tracer (in this case radionuclide) inventory in the water column within a bounded region decreases in a factor 𝑒−1, i.e. 0.37. These are relevant parameters for the water quality of a system since they provide the time scale for a radionuclide discharged into a water body to be transported out of it (Shen and Haas, 2004). The Red Sea was divided into 18 latitude bands 0.8◦ wide each one, from the north of BeM Strait to the entrance of the Gulf of Suez (see Fig. 1 for locations). One simulation was carried out for each latitude band. In these simulations an arbitrary radionuclide concentration was defined in the band surface layer (2 m thick) and the inventory within was calculated along one year. From a numerical fitting of the inventory time evolution to an exponential decay function, the flushing time was obtained. Radioactive decay and water/sediment interactions were not included in the calculations to estimate flushing times due to water circulation only. It must be commented that the Eulerian model does not include tidal currents and tidal residuals. However, as described in Periáñez (2020) tide effects are generally negligible. Moreover, the starting time of the simulation could affect the flushing time, since baroclinic currents change along the year. Thus, flushing times were evaluated for spring and fall as examples (this means that simulations for each band were carried out twice: with starting times at the beginning of each season) to investigate if there are significant differences. Results of these numerical experiments are presented in Fig. 8. An example of the numerical fitting of the radionuclide inventory within one band to a function of the form 𝐼(𝑡) = 𝐼0𝑒−𝑡∕𝑇𝑓(2) is also shown, where 𝐼0is the initial inventory and 𝑇𝑓is flushing time of the considered region. The spatial distribution of flushing times in fall and spring may be seen in the maps. The white zones in the northern Red Sea indicate that numerical fitting to the exponential curve could not be carried out. Eddies retain radionuclides in this area and a inventory decrease in the simulated time is not apparent. In general, flushing times are below 100 days in most of the sea and for both seasons. The most significant difference is that the highest flushing times are located in the central Red Sea in fall, while they are displaced to the north in spring. 4. Conclusions Two radionuclide transport models (Eulerian and Lagrangian) have been developed for the Red Sea, a region for which these models have not been described before. Both models include baroclinic circulation data obtained from HYCOM ocean model; the Lagrangian model can also include tidal currents, derived from a tidal model previously customized to the Red Sea. Both transport models also include water/sediment interactions, described in a dynamic way by means of kinetic transfer coefficients. A number of numerical experiments were carried out in order to illustrate the performances of models. The Eulerian model was applied to a relatively long simulation of a chronic release; showing that radionuclides are trapped in gyres and also that sediments along the coast (where water is shallower) are contaminated. These bands of marked sediments slowly increase their width with time. The Lagrangian model was applied to simulate hypothetical accidental releases occurring near the coast (in shallow water) and in the central Red Sea (where water is deeper) and also considering different release depths. The fate of radionuclides clearly depends on the release depth due to the significant vertical structure of baroclinic currents. All results presented in this paper consider a radionuclide source which is fixed at a given location. However, a source moving in the marine space could be considered after an appropriate modification of the code. There is not any practical or computational limitation to this. Finally, as another example, the Eulerian model was used to estimate flushing times over the Red Sea, which is a relevant parameter for water quality assessments. In general, it was found that they are below 100 days over most of the sea, presenting some seasonal variability. This work is an example on how specific radionuclide transport models could be customized to regions potentially exposed to an accidental release, in line with conclusions from the marine working group in IAEA MODARIA and MODARIA-II programs. Models are quickly and easily prepared for a given simulation, requiring few parameters and just the editing of some input files. Model running times are short, and they may be run in a laptop PC. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Appendix A. Lagrangian model Advection in a Lagrangian model is computed solving the following equation for each particle: 𝛥𝑥 =𝑢 𝛥𝑡 +𝜕𝐾ℎ 𝜕𝑥 𝛥𝑡 (3) 𝛥𝑦 =𝑣 𝛥𝑡 +𝜕𝐾ℎ 𝜕𝑦 𝛥𝑡 (4) where 𝛥𝑥 and 𝛥𝑦 are the changes in particle position (𝑥, 𝑦);𝑢and 𝑣are water velocity components in the west–east and south–north directions, respectively, at the particle position and depth; and for the moment when the calculation is done, since currents are changing in time. Derivatives of the horizontal diffusion coefficient (𝐾ℎ) prevent the artificial accumulation of particles in regions of low diffusivity (Proehl et al.,2005;Lynch et al.,2015). Constant diffusion coefficients are used in the present work, thus these terms are not relevant here. The use of constant diffusivities is just a simplification to speed up calculations. Of course more complex descriptions could be implemented, as the Smagorinsky scheme (Cushman-Roisin and Beckers,2011). Actually, Lagrangian transport models were applied to the Pacific Ocean (Periáñez et al.,2019b), some with constant and some with Smagorinsky diffusivities, providing comparable results. It must be commented that a first order accuracy equation is used to describe advection. However, Elliott and Clarke (1998) did not find improvements in results when a second order accuracy scheme was used to simulate the movement of surface drifters, with respect to the first order equation. Moreover, in marine transport problems, the effects of turbulence will mask any small errors in the advection scheme. Vertical currents in the sea are small, thus vertical advection is masked by vertical mixing due to turbulence. As a consequence, it is a common approach to neglect vertical advection in marine transport modeling. Water currents, 𝑢(𝑥, 𝑦, 𝑧, 𝑡)and 𝑣(𝑥, 𝑦, 𝑧, 𝑡), are the addition of tidal current, tidal residual, baroclinic current and local wind induced current. The maximum size of the horizontal step given by the particle due to turbulence, 𝐷ℎ, is (Proctor et al.,1994;Hunter,1987;Periáñez and Elliott,2002): 𝐷ℎ=√12𝐾ℎ𝛥𝑡 (5) 8 Fig. 8. Calculated flushing times (days) for fall and spring. An example of the time evolution of the radionuclide inventory within a latitude band and the numerical fitting (red line) to an exponential decay is shown in the top. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) in the direction 𝜃= 2𝜋𝑅𝐴𝑁, where 𝑅𝐴𝑁 is an uniform random number between 0 and 1. This equation gives the maximum size of the step. In practice, it is multiplied by another independent random number to obtain the real size at a given time and for a given particle. This is required to ensure a simulation of a Fickian diffusion process (Proctor et al.,1994;Hunter,1987). Similarly, the size of the vertical step is (Proctor et al.,1994;Hunter, 1987;Periáñez and Elliott,2002): 𝐷𝑣=√2𝐾𝑣𝛥𝑡 (6) given either towards the sea surface or the sea bottom; 𝐾𝑣is the vertical diffusion coefficients and 𝛥𝑡 is time step. Radioactive decay is treated in a stochastic way as well (Periáñez and Elliott,2002). Thus, a decay probability is defined as: 𝑝𝑑= 1 − 𝑒−𝜆𝛥𝑡 (7) where 𝜆is the radioactive decay constant. A new random number is generated. If 𝑅𝐴𝑁 ≤𝑝𝑑the particle has decayed and is removed from the computation. A similar stochastic method is applied to describe interactions between dissolved radionuclides and the bed sediments. Thus, if a particle is dissolved, the probability that it is adsorbed by the sediment is: 𝑝𝑎= 1 − 𝑒−𝑘1𝛥𝑡 (8) If a new random number 𝑅𝐴𝑁 ≤𝑝𝑎the particle is adsorbed by the sediment. For a particle which is fixed to the sediment, the probability that it is redissolved is: 𝑝𝑟= 1 − 𝑒−𝑘2𝜙𝛥𝑡 (9) and the same procedure follows. 𝜙is a correction factor that takes into account that part of the sediment particle surface may be hidden by other sediment particles. Thus, this part is not interacting with water. Full details may be seen in the references cited in the main body of the paper. Appendix B. Eulerian model The equation which gives the time evolution of radionuclide concentration in the dissolved phase, 𝐶𝑑(Bq/m3), is: 𝜕𝐶𝑑 𝜕𝑡 +𝜕𝑢𝐶𝑑 𝜕𝑥 +𝜕𝑣𝐶𝑑 𝜕𝑦 =𝐾ℎ(𝜕2𝐶𝑑 𝜕𝑥2+𝜕2𝐶𝑑 𝜕𝑦2)+𝐾𝑣(𝜕2𝐶𝑑 𝜕𝑧2) +𝛿𝑏(𝑘2 𝐿𝜌𝑠𝜙𝐴𝑠 𝜓−𝑘1𝐶𝑑) − 𝜆𝐶𝑑(10) 9 where 𝐴𝑠 is the radionuclide concentration in bottom sediments (Bq/kg), 𝑢 and 𝑣 are water velocities along the 𝑥 and 𝑦 axis and 𝐾ℎ and 𝐾𝑣 are, respectively, the horizontal and vertical diffusion coefficients. 𝐿 is the sediment mixing depth (the distance to which the dissolved phase penetrates the sediment), 𝜌𝑠 is the sediment bulk density (dry mass divided by wet volume) and 𝜙 is the correction factor introduced in the Lagrangian model description. Finally, 𝜓 is the thickness of the deepest water layer (as configured in HYCOM model). Standard values are used for these parameters since specific values for the Red Sea are not available (and in any case they would not be uniform over the sea). The same comments with respect to diffusion coefficients and neglecting vertical water velocities made above for the Lagrangian model are equally valid for the Eulerian one. The delta function is introduced to take into account that only the deepest water layer interacts with the bed sediment. Thus 𝛿𝑏= 1 for the deepest layer and 𝛿𝑏= 0 elsewhere. The equation for the temporal evolution of radionuclide concentration in the bottom sediment mixed layer is: 𝜕𝐴𝑠 𝜕𝑡 =𝑘1 𝐶𝑑(𝑏)𝜓 𝐿𝜌𝑠 −𝑘2𝐴𝑠𝜙−𝜆𝐴𝑠(11) where 𝐶𝑑(𝑏)means the radionuclide concentration in water at the deepest water layer, which is in contact with the bed sediment. Conservation of mass was carefully checked and errors resulted negligible. Details on the model can be seen in the references cited in the paper. References Bleck, R., 2001. An oceanic general circulation model framed in hybrid isopycnic-Cartesian coordinates. Ocean Model. 4, 55–88. Boon, J.D., 2011. Secrets of the Tide. Woodhead Publishing, USA. Cushman-Roisin, B., Beckers, J.M., 2011. Introduction to Geophysical Fluid Dynamics. Elsevier. El-Taher, A., Alashrah, S., Madkour, H.A., Al-Sayed, A., El-Erian, T.M., 2018a. Radionuclides distribution in marine sediment from Abu Soma Bay, Egyptian Red Sea coast. J. Environ. Sci. Technol. 11, 95–103. El-Taher, A., Zakaly, H.M.H., Elsaman, R., 2018b. Environmental implications and spatial distribution of natural radionuclides and heavy metals in sediments from four harbours in the Egyptian Red Sea coast. Appl. Radiat. Isot. 131, 13–22. Elliott, A.J., Clarke, S., 1998. Shallow water tides in the Firth of Forth. Hydrogr. J. 87, 19–24. Elliott, A.J., Wilkins, B.T., Mansfield, P., 2001. On the disposal of contaminated milk in coastal waters. Mar. Pollut. Bull. 42, 927–934. Hunter, J.R., 1987. The application of Lagrangian particle tracking techniques to modelling of dispersion in the sea. In: Noye, J. (Ed.), Numerical Modelling, Applications to Marine Systems. Elsevier, North-Holland, pp. 257–269. IAEA, 2004. Sediment Distribution Coefficients and Concentration Factors for Biota in the Marine Environment. Technical Reports Series 422, Vienna. Kobayashi, T., Otosaka, S., Togawa, O., Hayashi, K., 2007. Development of a nonconservative radionuclide dispersion model in the ocean and its application to surface cesium-137 dispersion in the Irish Sea. J. Nucl. Sci. Technol. 44, 238–247. Lynch, D.R., Greenberg, D.A., Bilgili, A., Mcgillicuddy, D.J., Manning, J.P., Aretxabaleta, A.L., 2015. Particles in the Coastal Ocean. In: Theory and Applications, Cambridge University Press, NY. Mamoney, M.H.E., Khater, A.E.M., 2004. Environmental characterization and radioecological impacts of non-nuclear industries on the Red Sea coast. J. Environ. Radioact. 73, 151–168. Parker, B.B., 2007. Tidal Analysis and Prediction. NOAA Special Publication NOS CO-OPS 3. Patzert, W.C., 1974. Wind-induced reversal in Red Sea circulation. Deep Sea Res. Oceanogr. Abstr. 21, 109–121. Periáñez, R., 2005. Modelling the Dispersion of Radionuclides in the Marine Environment. an Introduction. Springer. Periáñez, R, 2007. Chemical and oil spill rapid response modelling in the Strait of Gibraltar-Alborán Sea. Ecol. Model. 207, 210–222. Periáñez, R., 2009. Environmental modelling in the Gulf of Cadiz: heavy metal distributions in water and sediments. Sci. Total Environ. 407, 3392–3406. Periáñez, R., 2012. Modelling the environmental behavior of pollutants in Algeciras Bay (south Spain). Mar. Pollut. Bull. 64, 221–232. Periáñez, R., 2020. A lagrangian oil spill transport model for the red sea. Ocean Eng. http://dx.doi.org/10.1016/j.oceaneng.2020.107953. Periáñez, R., Abril, J.M., 2014. A numerical modelling study on oceanographic conditions in the former Gulf of Tartessos (SW Iberia): tides and tsunami propagation. J. Mar. Syst. 139, 68–78. Periáñez, R., Bezhenar, R., Brovchenko, I., Duffa, C., Iosjpe, M., Jung, K.T., Kobayashi, T., Liptak, L., Little, A., Maderich, V., Min, B.I., Nies, H., Osvath, I., Suh, K.S., de With, G., 2019a. Marine radionuclide transport modelling: Recent developments, problems and challenges. Environ. Model. Softw. 122, 104523. Periáñez, R., Bezhenar, R., Brovchenko, I., Duffa, C., Jung, K.T., Kobayashi, T., Lamego, F., Maderich, V., Min, B.I., Nies, H., Osvath, I., Outola, I., Psaltaki, M., Suh, K.S., de With, G., 2016a. Modelling of marine radionuclide dispersion in IAEA MODARIA program: lessons learnt from the Baltic Sea and Fukushima scenarios. Sci. Total Environ. 569/570, 594–602. Periáñez, R., Bezhenar, R., Brovchenko, I., Jung, K.T., Kamidara, Y., Kim, K.O., Kobayashi, T., Liptak, L., Maderich, V., Min, B.I., Suh, K.S., 2019b. Fukushima137Cs releases dispersion modelling over the Pacific Ocean. Comparisons of models with water, sediment and biota data. J. Environ. Radioact. 198, 50–63. Periáñez, R., Casas-Ruíz, M., Bolívar, J.P., 2013. Tidal circulation, sediment and pollutant transport in Cádiz Bay (SW Spain): a modelling study. Ocean Eng. 69, 60–69. Periáñez, R., Elliott, A.J., 2002. A particle tracking method for simulating the dispersion of non conservative radionuclides in coastal waters. J. Environ. Radioact. 58, 13–33. Periáñez, R., Pascual-Granged, A., 2008. Modelling surface radioactive, chemical and oil spills in the Strait of Gibraltar. Comput. Geosci. 34, 163–180. Periáñez, R., Suh, K.S., Min, B.I., 2016b. The behaviour of137Cs in the North Atlantic Ocean assessed from numerical modelling: Releases from nuclear fuel reprocessing factories, redissolution from contaminated sediments and leakage from dumped nuclear wastes. Mar. Pollut. Bull. 113, 343–361. Prandle, D., 1984. A modelling study of the mixing of137Cs in the seas of the European Continental Shelf. Philos. Trans. R. Soc. Lond. Ser. A 310, 407–436. Proctor, R., Flather, R.A., Elliott, A.J., 1994. Modelling tides and surface drift in the Arabian Gulf: application to the Gulf oil spill. Cont. Shelf Res. 14, 531–545. Proehl, J.A., Lynch, D.R., McGuillicuddy, D.J., Ledwell, J.R., 2005. Modeling turbulent dispersion on the North Flank of Georges Bank using Lagrangian particle methods. Cont. Shelf Res. 25, 875–900. Pugh, D.T., 1987. Tides, Surges and Mean Sea Level. Wiley, Chichester, p. 472. Saharty, E.A.A., Dar, M.A., 2010. The concentration levels of some isotopic radionuclides in the coastal sediments of the Red Sea, Egypt. Isot. Radiat. Res. 42 (1), 11–27. Shen, J., Haas, L., 2004. Calculating age and residence time in the tidal York River using three dimensional model experiments. Estuar. Coast. Shelf Sci. 61, 449–461.