1. Introduction

In the last years there has been an increasing interest in developing and improving numerical models aimed to perform simulations of the transport of radioactive releases in the marine environment, specially after Fukushima Dai– ichi nuclear power plant accident in 2011 [1].

Model domains may be small (typically up to a few hundred km) if these models are specifically designed as predictive tools to assess the effects of an accidental radioactive release in a region which is potentially exposed to a nuclear accident, as for instance a coastal region where one or several nuclear facilities are located or a marine area with intense shipping activities including transport of radioactive material or nuclear powered vessel transit. Model predictions would then be used to support decision–making after the accident. Recent examples to this type of model may be seen, for instance, in [2] for the Red Sea and [3] for the Arabian/Persian Gulf. These could be denoted as emergency models.

Other models have much larger spatial domains, being able to simulate the dispersion of radionuclides at full oceanic scales. These models can be used to assess the effects of a nuclear accident at long temporal scales after the incident. This was for instance the case with Fukushima releases in the Pacific Ocean after the 2011 tsunami (see for instance [4,5,6,7] among many others). A detailed review of works aimed at simulating the transport of Fukushima releases in the Pacific Ocean may be seen in [1]. Nevertheless, these models can also be used to assess the effects of chronic releases to the sea at large spatio–temporal scales, for instance from nuclear fuel reprocessing plants. This is the case with the releases to the Atlantic Ocean carried out from nuclear fuel reprocessing plants of Sellafield in the UK and La Hague in France [8,9,10,11].

The purpose of this paper is to summarize several modelling works which were carried out to simulate the transport of radionuclides at oceanic scales for the Atlantic, the Indian and the Pacific Ocean. In addition, it may be useful as a summary of the techniques which are required to carry out this kind of study.

In the case of the Atlantic the focus will be on the potential leakages of nuclear wastes dumped in the ocean. About 85 PBq of radioactive waste was deliberately dumped at more than 80 locations over the world [12]. About 46 PBq of activity, mainly in the Northeast Atlantic, were dumped between 1956 and 1982. A voluntary moratorium on radioactive waste dumping in the sea started in 1985 by the contracting parties of the London Convention. Dumping was finally banned in 1993. The International Atomic Energy Agency (IAEA) published a detailed report on these activities [13] in year 2015.

Several nuclear power plants (NPP) are operating along the shores of the northern Indian Ocean. These are: Kanupp NPP in Pakistan and Kakparar, Tarapar, Kaiga, Kudankulam and PFBR NPPs in India. A radionuclide transport model which could be used to simulate radionuclide transport from a hypothetical accident, as was done in the case of Fukushima accident, is described in [14]. This model can also provide insight into dispersion patterns resulting from different release scenarios (different release location and season). Thus, the model is a tool which may help to assess the consequences of a nuclear accident, but in addition it can also be used in a more process–orientated investigation of transport in the ocean. The main results of the work presented in [14] are summarized in this review.

As already mentioned, a number of papers have investigated the oceanic transport of radionuclides released after Fukushima accident in the Pacific Ocean (see review in [1]). Nevertheless, there are plans to release the water used to cool the reactors to the sea, since there is not more space to store it. This water is contaminated by radionuclides [15], which therefore may contaminate the marine environment. These planned releases have been recently simulated for ³H [16] and for different release scenarios. Results were that most of tritium would be rapidly diluted in the coastal waters off Japan, and then transported eastward along the Kuroshio Current extension. The radionuclide patch would arrive to the north-American Pacific coast in about 4.5 years. The case of oceanic transport of ¹³⁷Cs planned releases was studied in [17]; main results of this work are presented here.

The paper is organized as follows: theoretical foundations are given in section 2. A brief description of the general oceanic circulation is presented in 2.1 and the numerical models used to describe it are introduced in section 2.2. The Lagrangian radionuclide transport model used in these works is presented in section 2.3. Main results are finally summarized in section 3 for the Atlantic, Indian and Pacific Ocean separately.

2. Theory

2.1 Ocean Circulation

The large–scale transport of tracers in the ocean is dominated by the ocean circulation patterns, which are essentially defined by atmospheric forcing in addition to Coriolis force due to rotation of the Earth. The general features of ocean circulation are [18]:

These features may be seen in the map in Fig. 1. Circulation in the northern Indian Ocean is different to the northern Atlantic and Pacific due to two reasons: a significant part of the ocean is occupied by land and also it is affected by monsoons. Thus large scale circulation is seasonal. Essentially, during winter in the northern hemisphere a high pressure system is formed in central Asia, which results in northeast winds over the ocean. This is the so–called winter (or northeast) monsoon. This high pressure system weakens and in summer a low pressure system of thermal origin results which produces a reversal in wind direction. This is the so–called summer (or southwest) monsoon. The winter monsoon extends from November to February and the summer one from May to September. The remaining months are transition ones from one monsoon to the other. Circulation in the Indian Ocean, which for instance is analyzed in detail in [19,20] is affected by this atmospheric forcing.

Typical surface water circulation schemes for both monsoons are presented in Fig. 2. The summer monsoon map corresponds to August 1st and the winter monsoon map is obtained in December 31th. The names of the main currents are included in the figure. Since currents are seasonally reversing, in winter the overall current direction in the northern ocean is to the west: the East India Coastal Current (see Fig. 2) flows to the south, the Winter Monsoon Current to the west in the south of India and the West India Coastal Current to the north along the western India coast. The Somali Current flows to the south along east Africa coast. This pattern reverses during summer, with the Somali Current flowing northwards and a overall current direction to the east. Circulation in the Arabian Sea is anticyclonic in summer and cyclonic in winter [19].

Coastal upwelling in the Indian Ocean is different to the Atlantic and Pacific oceans, where upwelling occurs along the eastern coasts. The strongest upwelling area in the Indian Ocean is placed along Somalia and Arabia coasts, when the SW monsoon leads to Ekman transport away from the coasts [21].

2.2 Ocean Circulation Models

Three ocean circulation models were used to obtain water currents for each considered application, depending on the availability of data for each region and in the required temporal frame. Also, the fact that ocean model results were adequately compared with circulation data was taken into account. Such three–dimensional water currents were downloaded from each model data server in each corresponding web page, as described below.

2.2.1 OFES Model

OFES (Ocean global circulation model For the Earth Simulator) model (http://www.jamstec.go.jp/esc/research/ AtmOcn/product/ofes.html cite note-1) was used to obtain circulation in the Atlantic Ocean. This model is operated by JAMSTEC (Japan Agency for Marine-Earth Science and Technology).

The considered domain in the north Atlantic extends from 50°W to 25°E in longitude and from 45.1°N to 75.1°N in latitude. A comparison of model performances in several regions of the ocean (including the North Atlantic) is given in [22]. Horizontal resolution is 0.1° with 54 vertical levels with increasing thickness from the surface towards the seabed. Monthly mean circulation has been used. Surface water circulation for January 2008 is presented in Fig. 3 as an example. The Norwegian Coastal Current (NCC) flows along the Norwegian coast to the north. The East Greenland Current (EGC) flows southwards along the east Greenland coast. Finally, the North Atlantic Current (NAC), deriving from the Gulf Stream, is seen in the central Atlantic with several eddies and meanders.

2.2.2 HYCOM model

HYCOM (Hybrid Coordinate Ocean Model) model [23] was used to obtain circulation in the Indian Ocean. It is a primitive equation general circulation model with horizontal resolution 0.08^◦ and 40 vertical layers. HYCOM web page (https://www.hycom.org/) gives several examples of applications over the world, as well as a detailed description of the model. Particularly, the model is forced by winds, short and long wave radiation and freshwater inflow through the main rivers. The modelling system uses the Navy Coupled Ocean Data Assimilation (NCODA) method [24]. Water circulation in the Indian Ocean shown in Fig. 2 was actually derived from this model. Only one of each 100 current vectors were drawn for more clarity. The considered domain extends 3–30° in latitude and 48–100° in longitude.

2.2.3 FORA Model

FORA model [Four-dimensional Variational Ocean Re- Analysis for the Western North Pacific (FORA- WNP30)] was used for the Pacific Ocean. This constitutes the first dataset covering the western North Pacific over the last three decades (1982-2014) at eddy-resolving resolution. The model is a cooperative work of JAMSTEC and the Meteorological Research Institute, Japan Meteorological Agency (JMA/MRI) using the Earth Simulator [25]. The domain extends 117°E–160°W and 15°N–65°N in longitude and latitude, respectively. Horizontal resolution is 0.1° with 54 vertical levels, whose thickness increase from the surface to the seabed. Monthly climatological data from 2011 to 2014 were used.

An example of surface water circulation (averaged value for March 2011) is provided in Fig. 4. Large scale circulation in the western Pacific Ocean is dominated by the interaction between the Kuroshio and Oyashio currents. The Kuroshio Current is the western boundary current in the north Pacific, flowing along the coast of Japan towards the north and which curves to the central Pacific Ocean, then forming the so–called Kuroshio Extension. The Oyashio Current is a cold current flowing from the north. These currents converge in the coastal waters off Japan. Such convergence leads to the generation of unsteady eddies in the area.

2.3 Lagrangian Transport Modelling

The transport model used is Lagrangian as commented before. Eulerian models, in which a differential equation for the temporal evolution of radionuclide concentrations over the domain is solved, could be used as well. A detailed description of the relative advantages of each approach may be seen in the review in [1].

In a Lagrangian model 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 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 model grid cell and compartment (surface water, deep water and sediment) and the number of units (Bq) which corresponds to each particle.

Advection in a Lagrangian model is computed solving the following equation for each particle:

where L_x and L_y are the changes in particle position (x, y); u and v 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 (K_h) prevent the artificial accumulation of particles in regions of low diffusivity [26,27]. Different schemes which may be used to fix values for the diffusion coefficients are described below.

A first order accuracy scheme was applied to solve advection. However, Elliott and Clarke [28] did not find improvements in results if a second order accuracy scheme were used. Moreover, turbulence masks small errors in the advection scheme [28].

Turbulent mixing is simulated using a stochastic method. The maximum size of the horizontal step given by the particle due to such process, D_h, is [29,30]:

in the direction θ = 2πRAN , where RAN is an uniform random number between 0 and 1 and Δt is time step in the Lagrangian model. This equation gives the maximum size of the step. The real size at a given time and for a given particle is obtained multiplying the equation by another independent random number. This procedure is required to ensure that a Fickian diffusion process [29] is simulated.

Similarly, the size of the vertical step given by each particle is [29,30]:

given either upward or downward. In this equation K_v is the vertical diffusion coefficient.

The Smagorinsky scheme (see for instance [31]), used in the models for the three oceans, describes the horizontal diffusion coefficient as a function of the grid spacing, Δx and Δy, and the horizontal shear of currents:

The vertical diffusion coefficient is set as a constant typical value K_v = 1.0 × 10⁻⁵ m²·s⁻¹ [32,33] in the cases of Atlantic and Pacific ocean models. However, a variable vertical diffusivity was tested in the Indian Ocean, following the description in [34] which was used to describe radionuclide transport in the Pacific Ocean immediately after Fukushima accident:

where z is water depth below the sea surface. This way increased turbulence in the mixed layer of the ocean is described, while it is smaller in the pycnocline and deep ocean. Using a constant mixed layer thickness is just an approximation to speed up calculations; of course this thickness is seasonally variable. For instance, its climatological mean ranges from approximately 30 to 100 m in the Indian Ocean [35].

Radioactive decay and exchanges of radionuclides between water and sediment are described through stochastic methods as well. A dynamic approach is applied to describe water/sediment interactions, thus a kinetic coefficient k₁ describes the transfer of radionuclides from water to sediment and a coefficient k₂ governs radionuclide release from the sediment.

In the case of radioactive decay, a decay probability during a time step Δt is defined as [30]:

where λ is the radioactive decay constant. A new independent random number is generated. If RAN ≤ p_d then the particle decays and it is removed from the computation.

The probability that a dissolved particle is adsorbed by the sediment during Δt is defined as:

If a new independent random number results RAN ≤ p_a, then the particle is adsorbed by the sediment. The probability that a particle which was fixed to the sediment is redissolved during Δt is written as:

and the same procedure follows. ϕ is a correction factor that takes into account that part of the sediment surface is hidden by surrounding sediments. Thus, this part is not interacting with water [36].

The number of units corresponding to each particle, R is deduced from the number of particles released each time step and of the radionuclide release rate from each nuclear facility (which may be constant in time or not). Typically, several millions particles are used in each simulation.

The concentration of radionuclides in a given water layer, k, of each grid cell (i, j), C_k(i, j), is:

where ΔxΔy gives the cell surface, N_k(i, j) is the number of particles in the water layer k of cell (i, j) and d_k is the water layer thickness.

The radionuclide concentration in the bed sediment of cell (i, j) is:

where N_sed(i, j) is the number of particles in the bed sediment of cell (i, j), L is sediment thickness (set as 0.05 m) and ρ_s is sediment bulk density:

where ρ_m = 2,600 kg·m⁻³ is mineral particle density and por is sediment porosity. A number of parameters are defined within the model code, whose values are selected from standard ones or previous works. Thus, porosity is set as por = 0.6, the desorption kinetic rate as k₂ = 1.16 × 10⁻⁵ s⁻¹ and the sediment correction factor as ϕ = 0.1. A justification for parameter values may be seen in the papers describing the model applications. Of course, these values may be changed if desired.

The adsorption rate k₁ is deduced from the desorption rate k₂ and the radionuclide distribution (or partition) coefficient k_d for ocean waters, which can be obtained from the IAEA [37], following the procedure described in [38] for instance. The desorption rate k₂ given above was obtained from the experiments described in Nyffeler et al. [39], value which has been used in other modelling works (for instance [10,11,14,40] and the review in [1]).

Lagrangian transport models also allow the calculation of ages of the released radionuclides. The age is defined as the time elapsed since a given Lagrangian particle was released into the sea. It may provide useful information about circulation in a marine region. Age is a Lagrangian concept and is obtained simply attaching a clock to each particle, which is started in the moment when the particle is released. The age-averaging hypothesis, as introduced by Deleersnijder [41], was used: the mean age of particles which are within each grid cell is defined as the mass-weighted arithmetic average of the ages of the particles present there. The mass of particle i is defined as its radioactivity content R_i, thus the mean age, < age >, in a given grid cell is:

where N is the number of particles within the cell and τ_i is the reading of the clock attached to particle i. If the values of R_i were the same for all particles (as it occurs in the examples here presented), then the mean age is the arithmetic means of the clock readings:

The computer code was written in FORTRAN and runs on a desktop PC working over Ubuntu operating system. Running times depends on the simulation time horizon and spatial resolution of the domain, but typically ranges from two hours for the shortest simulations (90 days in the Indian Ocean) to four days for the 10 year long simulations in the Atlantic Ocean.

3. Results

The numerical transport model described in the previous section was applied to a number of radionuclides, as already indicated: ¹²⁹I [9], ¹³⁷Cs [10,17,34,38], ²³⁶U [11], ⁹⁰Sr [40] and ^239,240Pu [38]. Nevertheless, only ¹³⁷Cs results are summarized here. The ¹³⁷Cs distribution coefficient k_d was fixed as 4.0 m³·kg⁻¹, which is a recommended value by [37].

3.1 Atlantic Ocean

An exercise on modelling the radiological impact of radioactive wastes dumped in the Arctic was carried out under the auspices of the IAEA [42]. But an exercise on the dispersion of radionuclides leaking in deep areas of the Northeast Atlantic or in areas affected by strong currents (like the English Channel) was not carried out until [10].

A numerical experiment was carried out on the site where wastes were released in 1969, with coordinates 49.08°N, 17.08°W and 4,600 m water depth. Seven countries dumped 8.2 × 105 GBq of β − γ emitters. Charmasson [43] reports that ¹³⁷Cs accounts for 1.92% of such amount. It was supposed that ¹³⁷Cs was released in January 2012. Thus, correction by radioactive decay leads to a release of 5.86 × 10¹² Bq of ¹³⁷Cs. A worst case scenario, consisting of an instantaneous release of the whole inventory of this radionuclide, was simulated. Radionuclides were also supposed to be released in dissolved form to deal with a worst situation.

As an example, the calculated distributions of ¹³⁷Cs in the water column at the end of 2022 is presented in Fig. 5. Four water layers have been considered to calculate the radionuclide concentrations.

It may be seen that below 200 m, radionuclides are distributed over virtually all the north Atlantic Ocean. Radionuclides reach the surface layer (0–200 m) in the shelf seas (Celtic, Irish, North seas), but concentrations are several orders of magnitude below background levels. Some snapshots of the ¹³⁷Cs distribution in sediments are presented in Fig. 6. Since the release occurs close to the seabed, the sediment adsorb a significant fraction of ¹³⁷Cs. Thus, in some areas levels approach to those measured for instance in the English Channel due to the releases from Sellafield and La Hague. Initially, only sediments in the area of the dumping site are affected. As ¹³⁷Cs slowly moves upwards in the water column and is transported towards shelf waters, sediments of the shelf seas are also affected by the leakage as well.

A second experiment was carried out for the 1962 dumping site in the English channel: 49.83°N, 2.30°W, where water depth is 120 m. Activity dumped by two countries was 4.7 × 10³ GBq of β − γ emitters. After correction by radioactive decay to 2012 (the same date and conditions for the leakage as before were assumed, as well as the 1.92% ¹³⁷Cs fraction given in [43]), this number implies that 2.86 × 10¹⁰ Bq of ¹³⁷Cs were released.

Calculated concentrations in surface waters and bed sediments at the end of 2012 and 2017 are presented in Fig. 7 as examples. Radionuclides travel along the Channel towards the North Sea. Then part of the radionuclides enter the Danish Straits and the remaining follow the coast of Norway towards the Arctic Ocean. Bed sediments adsorb radionuclides as contaminated water travels above them. Calculated concentrations are two orders of magnitude lower than those of the first experiment (1969 dumping site) since the dumped inventory in this site in 1962 is two orders of magnitude lower. Concentrations in water are below today background levels.

3.2 Indian Ocean

A number of numerical experiments for the Indian Ocean were carried out in [14], in which ¹³⁷Cs releases of similar magnitude to those of Fukushima NPP after the accident were simulated for winter and summer monsoons. These hypothetical accidents were simulated for four NPPs located along the northern Indian Ocean coasts. Releases were supposed to last 90 days, as the acute phase of Fukushima accident. Some results are summarized here.

Calculated ¹³⁷Cs concentrations in the water surface layer, inventory in the deep layer (from the bottom of the surface layer to the seabed) and concentrations in the bed sediments can be seen in Fig. 8 for Tarapur NPP as an example. These maps are obtained after 90 days, i.e., just at the moment when releases stop.

As a general result it was found that contamination of the sediment is quite limited to the coastal area close to the release point. In the case of Tarapur NPP (Fig. 8) contamination extends to the border of the continental shelf, which is of the order of 100 km wide here, the widest in all the coast around Pakistan and India (only in the northern Bay of Bengal is similar). The narrow continental shelf prevents the offshore contamination of the sediment.

The radionuclide distributions in water resulting from releases starting in summer and winter monsoons are in general not so different as could be initially expected from the different circulation schemes, which may be seen in Fig. 2. The largest differences appear for the releases from Kudankulam NPP, which is located in the most southern point of India. Currents are strong in this region (Fig. 2) and radionuclides, once they reach the nearly flowing WMC or SMC, are efficiently transported to the west or east respectively. In the case of releases starting in summer, they are mainly transported to the central Bay of Bengal, while winter releases tend to flow northwards along the western India coast [14].

The radionuclide distributions in the deep water layer reproduce the obtained in the surface. Thus, the surface and deep circulations in the considered regions of the Indian Ocean are similar.

The age distributions resulting from releases occurring in Tarapur NPP are presented in Fig. 9. The colored parts of the map indicate the extension of contamination (age in a grid cell where there are not any particles is zero). These extensions are similar. However, the calculated age distributions help to visualize the circulation of pollutants. Thus, if releases start with summer monsoon pollutant circulation tends to be cyclonic: age increases progressively from Tarapur NPP towards the northwest, then along Africa coast to the southwest and finally towards the northeast, with particles traveling to the release point. This pattern is reversed to anticyclonic in the case of releases starting in winter monsoon, as may be seen in Fig. 9, where both schemes are represented by the black arrows. This seems to be in contradiction with circulation patterns in the Arabian Sea [19] consisting of anticyclonic in summer and cyclonic in winter. The key to understand this is that Fig. 9 is not a snapshot of winter/summer circulation. For instance, for the winter simulation the radionuclide release is supposed to start in December 15th. Then, the release lasts for three months and nine additional months (without releases) are simulated. Thus, these maps represent an integration of circulation during a whole year (with release only in the first 3 months). The final result is a apparent anticyclonic overall transport of the winter release, for instance, which seems to contradict the well–known circulation scheme, correctly reproduced by HYCOM ocean model (Fig. 2).

3.3 Pacific Ocean

The total amount of ¹³⁷Cs contained in the tanks which store the contaminated water used to cool reactors was fixed as 1.2 × 10⁹ Bq in 2022, as described in detail in [17]. In order of assess the ³H enhancement in the Pacific Ocean after releasing this stored water Zhao et al. [16] carried out 4 simulations (since the release procedure has not been announced yet) using different release durations. These release durations were from one month to 10 years. Maximum peak concentrations are obtained in the shortest release case. Thus, a ¹³⁷Cs release lasting one month was simulated to deal with the worst case situation.

Contamination of the bed sediment is limited to a local area along Fukushima coast, and does not significantly extend offshore [17]. Maximum concentrations calculated by the model are of the order of 10⁻² Bq·kg⁻¹, which is two orders of magnitude below concentrations measured in sediments after Fukushima accident and also lower than background ¹³⁷Cs concentrations in sediments due to global fallout (of the order of 10⁰ Bq·kg⁻¹). Sediments buffer radionuclides [44,45], which are slowly released back to the water column and thus the sediment acts as a long–term delayed radionuclide source. Actually, ¹³⁷Cs is still present after one year. This may be seen in Fig. 10: a band of sediments containing some ¹³⁷Cs can be observed along the coast. Note that this slow release of radionuclides from the sediment can be simulated since a dynamic model is used to describe water/sediment interactions. It is not possible with an equilibrium model based on the distribution coefficient, k_d, concept [45].

Calculated concentrations in surface waters are presented in Fig. 11 at 60 days, 120 days, 1 year and 2 years after starting the releases of Fukushima stored waters. It can be clearly seen that the patch of traced water travels offshore from Japan coast due to the Kuroshio current. The highest concentrations are located around the release area, as could be expected. However, this region remains with the highest concentrations even a couple of years after the releases have finished. This is due to the role of sediment as long–term delayed source of previously released radionuclides, as already commented. However, it should be noted that the maximum concentration 60 days after starting releases is about 0.4 Bq·m⁻³, value which is below the order of magnitude of pre-Fukushima background in the Pacific Ocean [46]. After one year the peak concentration is reduced in one order of magnitude, being ∼10⁻² Bq·m⁻³. The role of eddies present in Kuroshio current in producing turbulent mixing in the release can be clearly seen in Fig. 11.

The Lagrangian model, in its application to this particular study, was completed with a dynamic foodweb biota uptake model which was integrated [17] within it. Such basic dynamic foodweb model consists of four species [47]: phytoplankton, zooplankton, non-piscivorous and piscivorous fish. The ¹³⁷Cs concentrations in biota resulting from the releases of stored water was analyzed in [17].

It was found that due to the dynamic nature of the biota uptake model, radionuclide concentration increase in fish is slow and, consequently, peak concentrations are not found immediately after the release. Nevertheless, peak concentrations are very low. In non–piscivorous fish they are of the order of 10⁻³ Bq·kg⁻¹ (wet weight) 30 days after the release starting time and increase to 10⁻² Bq·kg⁻¹ (wet weight) after 120 days. After 1 year, peak concentration is of the order of 10⁻³ Bq·kg⁻¹ (wet weight). Finally, the peak magnitude decreases one more order of magnitude 10 years after the release start [17].

Peak concentrations in piscivorous fish are one order of magnitude lower that those of non– piscivorous fish in the initial phase of the releases, but after one year they remain of the same order of magnitude. In any case, calculated maximum concentrations are below pre–Fukushima background levels in biota [17].

4. Conclusions

A Lagrangian transport model is described. It may be used to assess the transport of radioactive releases in the sea at oceanic scales, due to waste disposal or to nuclear accidents in any nuclear facility. The model includes advection by three–dimensional currents, three–dimensional turbulent mixing, radioactive decay and, finally, interactions of radionuclides between water and sediments. A dynamic model based on kinetic transfer coefficients is used for such purpose. A stochastic method is used to numerically solve the processes of mixing, decay and water/sediment interactions. Water currents are obtained from operational global ocean circulation models. The numerical model was applied to simulate ¹³⁷Cs transport in the Atlantic, Indian and Pacific oceans. Nevertheless, it can be applied to any other radionuclide.

In the case of the Atlantic, the transport of radionuclides potentially leaked from nuclear dumped waste was analyzed. Some hypothetical accidents in NPPs located along the northern Indian Ocean coast were simulated as well. Finally, the release of stored water used to cool Fukushima reactors after the accident in 2011, as planned by the government of Japan, was studied in the Pacific Ocean.

The model is a useful tool which may be applied to support decision–making after an accident, and also to assess the long–term environmental implications of chronic or planned radionuclide releases from nuclear facilities.