The POSEIDON-R compartment model for the prediction of transport and fate of radionuclides in the marine environment

Graphical abstract


Method details
Despite the rapid development of numerical hydrodynamic models, the compartment (box) models continue to play an important role in the modelling of transport and fate of radioactivity in marine environments because of its relative simplicity and flexibility [1]. Originally, POSEIDON model [2] was developed to assess the radiological consequences of instantaneous or continuous releases of a mixture of radionuclides in European sea waters. For this purpose a box modelling approach as described in [3] was adopted. The model calculated the radioactivity transfer in the water, bottom sediments and biota to estimate doses for humans resulting from the consumption of contaminated seafood. The modified model became a part of European Decision Support System for emergency response to nuclear accidents RODOS [4,5] under the name POSEIDON-R. In this model a transfer of radionuclides to marine organisms was described by means of a dynamical uptake model BURN [6], taking into account the trophic level of the organisms. Simultaneously, the list of sources of radionuclides due the accidental releases was also essentially extended.
In the paper we describe an advanced version of POSEIDON-R that is supplemented with a new biological uptake model describing migration of activity through the pelagic and benthic food webs [7] and modified dose modules with extended pathways for human exposure from marine releases of radioactivity. The model was applied to the northeast Atlantic shelf seas and the modelling results presented in this work were compared with measurements.

Dispersion of activity in water and sediments
The POSEIDON-R 3D compartment model simulates transfer of radioactivity in the water column and bottom sediment (Fig. 1). The water column compartment is vertically subdivided into layers. The suspended matter is settling in the water column. The radionuclide concentration in the water compartment is governed by a set of differential equations that describe i. the temporal variation in the concentration, ii. the exchange of radionuclides between adjacent compartments and between radionuclides in suspension and in the bottom sediment, iii. and radioactivity sources and decay. Exchanges among the water column layers are described by radionuclide fluxes due to advection, sediment settling, and turbulent diffusion. The bottom sediments are divided into three layers, and the transfer of radioactivity between the upper sediment layer and the water column resulting from resuspension, diffusion and bioturbation, and between the upper and middle sediment layers, resulting from diffusion only, are described as shown in Fig. 1. Downward burial processes operate in all three sediment layers.
The POSEIDON-R equations are obtained by averaging the three dimensional transport equations (see [8]) for the dissolved radionuclide concentration in the water column and the concentration in three layers of the bottom sediment assuming that partitioning between dissolved and particulate fractions is described by a distribution coefficient K d (m 3 t À1 ). The equations for the water column layers, upper and middle sediment layer read as follows: Here C 0ik is the spatially averaged concentration of radionuclide (Bq m À3 ) in the water column layer i of box k; i = 1 corresponds to the near bottom water column layer; C 1k and C 2k are the averaged concentration of radionuclide (Bq m À3 ) in the upper and middle sediment layers of box k, respectively; l (y -1 ) is the radionuclide decay constant; F ijkm is the water flux (t y -1 ) from layer i of box k to layer j of box m; V 0ik is the volume (m 3 ) of layer i of box k; h ik is the thickness (m) of the water column layer i of box k ; L t,k and L m,k are the thicknesses (m) of top and middle bottom sediment layers of box k, respectively; Q sik is the source of the activity (Bq y -1 ) in layer i of box k; g 0ik . . . g 5k are the transfer rates: g 0ik is a transfer rate of radioactivity due to sinking of the suspended sediments from the upper layer in the water column in box k,g 1ik is a transfer rate of radioactivity due to sinking of the suspended sediments to the lower layer in the water column or to the top sediment layer in box k, g 2k is a transfer rate of radioactivity from the top sediment layer to the near bottom layer in the water column in box k, g 3k is a transfer rate of radioactivity from the top to middle sediment layer in box k, g 4k is a transfer rate of radioactivity from the middle to top sediment layer in box k, g 5k is a transfer rate of radioactivity from the middle to deep sediment layer in box k. For the surface water layer, these coefficients are as follows: For the layers in the water column below the surface water layer, the coefficients are defined as: In the near bottom layer located in the water column just above the bottom sediment, the coefficients are defined as: Here L b,k is the thickness (m) of the bottom boundary layer in the water column of box k, SS k is the concentration of suspended sediments (t m À3 ) in box k, obtained from observations or model simulation, W gk is the settling velocity (m y À1 ) of suspended sediments in box k calculated as function of particles size; SSW k =SS k W gk is the sedimentation flux (t m À2 y À1 ) in box k; D is the coefficient of vertical diffusion (m 2 y À1 ) in the bottom layers; B k is the coefficient of bioturbation (m 2 y À1 ) in the top sediment layer of box k; e k is the porosity of the bottom sediment in box k; r k is the density of sediment particles (t m À3 ) in box k; and the coefficient R k is defined as The mass concentration of radioactivity in dry sediment for top and middle sediment layersC 1k and C 2k (Bq kg À1 ), respectively, is calculated as It was found [7] that the standard parameterization of exchange between water column and bottom sediments, which includes sorption/desorption reactions, molecular diffusion and bioturbation, downplays fluxes of activity between water column and bottom sediments in shallow areas with energetic currents. The same difficulties occurred for the Irish and North Seas with strong tidal currents. Therefore, the model was complemented by a parameterization of resuspension mechanism allowing more accurate description of the upper bottom sediments self-cleaning. The resuspension was parameterized as an additional flux directed from the upper sediment layer to the water (last term in the Eq. (12)), which is proportional to the near-bottom velocity in shallow boxes. Parameter β 0 was estimated as 2 Á 10 À8 and 4 Á 10 À8 for seas without and with strong tidal currents, respectively, |V k | is the velocity module (m yr À1 ) in near bottom water layer of box k defined as a sum of all fluxes of water between the near bottom water layer of box k and neighbouring water layers F kn divided by areas of faces between these layers S kn Equations for the daughter products are similar to Eqs. The recent biota model intercomparison [9] shows that dynamic biota models, which handle situations out from equilibrium, perform better than equilibrium models for the radioecological dose assessment after nuclear accidents. The biota model in POSEIDON-R is a dynamic food web model based on the approach developed by Heling et al. [6]. In the food web model, marine organisms are grouped into classes according to trophic level and species type (Fig. 2). Radionuclides are also grouped into classes according to the fish tissue type in which they are preferentially accumulated (e.g., 137 Cs tends to accumulate in muscle). These simplifications allow for a limited number of standard input parameters. The scheme of transfer of radionuclides through the marine food web is shown in Fig. 2. The different food chains exist in the pelagic zone and in the benthic zone. Pelagic organisms are grouped into a primary producer (phytoplankton) and consumers: zooplankton, forage (non-piscivorous) fish and piscivorous fish. The pelagic food web was implemented in the compartmental POSEIDON-R model [4,10,11]. The benthic food web includes three primary pathways for radionuclides: (i) transfer from water to macroalgae, then to grazing invertebrates; (ii) transfer through the vertical flux of detritus and zooplankton faces to detritus-feeding invertebrates; and (iii) transfer through contaminated bottom sediments to deposit-feeding invertebrates. External boxes in Fig. 2 show the concentrations of radionuclides in the water and in the organic deposit, which is in instantaneous equilibrium with the upper layer of the bottom sediment, calculated by the above described POSEIDON-R model. In the benthic food chain the radioactivity is transferred from the deposit feeding invertebrates to the demersal fish, and to the bottom predators. The components of this system are crustaceans (e.g detritus-feeders), molluscs (filter-feeders) and coastal predators feeding in the whole water column of shallow coastal waters [7]. Along with the food web, all organisms take radionuclides directly from water.
Due to the rapid uptake from water and the short retention time of radioactivity, the concentration of radionuclides in phytoplankton is calculated using the Biological Concentration Factor (BCF) approach [12]. For the macroalgae, a dynamic model is used to describe radionuclide concentrations due to the longer retention times where C w and C ma are the radionuclide concentration in the water and macroalgae, respectively, CF ma is the corresponding BCF, T 0.5,ma is the biological half-life of the radionuclide in the macroalgae, t is the time. The concentration of a given radionuclide in other considered marine organisms is described by the following differential equation: where C i and C f,i are the radionuclide concentration in the i-th marine organisms and their food, respectively, a i is the food extraction coefficient (assimilation rate), b i is the water extraction coefficient, K f,i is the food uptake rate, K w,i is the water uptake rate and T 0.5,i is the biological half-life of the radionuclide in the organism. The activity concentration in the food of a predator C f is expressed by the following equation, summing for a total of n prey types, where C prey,i is the activity concentration in prey of type i, P prey,i is preference for prey of type i, drw pred is the dry weight fraction of predator, and drw prey is the dry weight fraction of prey of type i. The index "0 00 corresponds to the organic deposit in bottom sediments. Values of the model parameters are discussed in [7] and are given in Tables 1-3.
It is well known that the uptake of caesium and strontium decreases with increasing salinity due to the increase in concentration of competing ions of potassium and calcium, respectively. For caesium it  [7]. The radionuclide transfers among marine food web compartments are given for 11 types of marine organisms. was taken into account when introducing the salinity-dependent correction factor F K for phytoplankton and macroalgae because caesium enters the food web primarily through the lowest trophic level whereas the contribution of direct uptake from water is minor [13]. Instead of using fixed concentration factors, the BCF for 137 Cs is related to potassium concentration via the electrochemical competition for which the parameters are based on laboratory experiments with marine plants. Then the BCFs for phytoplankton and macroalgae can be expressed by:  Table 3 a See text for definitions of parameters.  where BCF Ã ph = 20 L kg À1 and BCF Ã ma = 50 L kg À1 are standard BCFs for 137 Cs in marine environment [12], K + is the potassium concentration (mgÁL À1 ) and T is temperature (K).The concentration of potassium varies from the typical value for rivers of 0.5 mg L À1 till the typical value in the marine environment of 400 mg L À1 . For water with a K + concentration of above 1.5 mg L À1 , the potassium concentration could be linked to the salinity using the following linear relationship: where S is the salinity (g L À1 ). For strontium, the direct gill uptake is more important due to the lack of bioaccumulation through the food web. The gill extraction coefficient b i is based on empirical correlations derived from measured equilibrium levels in the seas and can be expressed by: According to the review of radiological data [14,15], every radionuclide is mainly accumulated in a specific tissue (target tissue). It can be assumed that the target tissue controls the overall elimination rate of the nuclide (T 0.5 ) in the organism. The radioactivity in the food for the predator is then the activity concentration in the target tissue diluted by the remaining body mass of the prey fish, calculated by multiplying the predicted level in the target tissue by its weight fraction. To calculate the concentration in the edible part of fish (flesh) from the calculated levels in the target tissues, a target tissue modifier (TTM) is introduced. This is also based on tissue distribution information as reported by [14,15]. Values of described parameters for the dynamic food-chain model are listed in Table 3.

Sources of activity
The POSEIDON-R model can deal with four types of routine and accidental radioactive releases: (i) atmospheric deposition directly on the sea surface; (ii) runoff of land deposited radionuclide; (iii) point sources associated with routine releases of nuclear facilities, located either directly at the coast or inland at river systems; (iv) point sources associated with accidental releases located in any box of the model domain.
For coastal discharges, it is useful to provide a more detailed description in the area close to the release point. For this purpose, the additional "coastal" boxes are nested into the large ("regional") boxes in the box system of the considered region. There are some assumptions and restrictions to the approach. These are as follows: (i) a coastal box has one vertical layer for the water column; (ii) a coastal box interacts with the surface layer of the surrounding regional box only, the depth of a coastal box is therefore less or equal to that of the surface layer of regional box; (iii) the exchange fluxes with the adjacent regional box are equal in both directions, i.e. only lateral diffusion is taken in account; (iv) only one coastal box can be added per regional box; (v) a coastal box contains at most one source of radioactivity. When calculating the radionuclide concentration in fish in small coastal boxes, the random fish migration should be taken into account. Following [7], the right hand side of Eq. (20) for radionuclide concentration in fish, both in the inner and outer compartments, is extended by the term ðC out f ish À C in f ish Þ=T migr for the coastal compartment and by the term ÀðC out f ish À C in f ish Þ=ðdT migr Þ, for the outer compartment. Here T migr is the characteristic time of fish migration from a coastal compartment, depending on compartment scale and fish species, and d is the ratio between the volumes of the outer and the coastal compartments.
POSEIDON-R has also the possibility to deal with off-shore point releases (e.g. for evaluation of the impact of sunken vessels, nuclear submarines, and off-shore waste dumping). In that case, it is possible to use a so-called "local" box. The off-site local boxes have the following features: (i) a local box can be placed at any point in the surrounding regional box at any depth; (ii) the volume and thickness of the local box are calculated as proportional parts of the outer regional box; (iii) as in the case of the coastal box, the exchange flows between the local box and the surrounding regional box are assumed to be equal in both directions.

Numerical solution
The problem is described by a set of ordinary differential equations, which may be written in a vector-matrix notation as: where C is the concentration vector; A is the coefficient matrix that includes water fluxes between boxes, parameters of the food-chain model, etc; Q re is the vector for the release term.
Step-like variations of the release in time are assumed, and the implicit Matrix Exponential Method [16] is used to solve a set of Eq. (26). A brief description of the method reads as follows: An unforced (homogeneous) matrix equation, has the solution: Over a computational time interval ðt n ; t nþ1 Þ this solution can be expressed as or, letting t be a variable representing the step-size, The matrix exponential e At is defined operationally by a truncated Taylor series where I is the identity matrix (diagonal elements = 1, off-diagonal elements = 0) with the same number of rows and columns as A. For a forced system (26) the general incremental solution is: whose exact solution in the case where Q is constant over the step-length is: This is the basic equation of the method. The symbol A À1 is the inverse of matrix A (i.e., Along the instant values of the radionuclide concentration there is a need to compute the time integrated concentration for effective committed dose estimation: For t 2 t n ; t n þ t Introducing the notation of the matrices: E2ðtÞ E3ðtÞ :¼ e At À I Eqs. (33) and (37) were rewritten as: Cðt n þ tÞ ¼ E1ðtÞCðt n Þ þ E2ðtÞQ ðt n Þ; ICðt n þ tÞ ¼ ICðt n Þ þ E2ðtÞCðt n Þ þ E3ðtÞQ ðt n Þ: Matrices E1, E2, E3 are linked with the recurrence relation: E2ðtÞ ¼ E3ðtÞA þ tI; The series in (42) can be substituted by the sum with requisite precision.

Dose module
The POSEIDON-R model includes dose module to assess individual and collective doses to the population due to the regular and accidental releases of radionuclides. The exposure pathways that are considered in the model include: internal exposure through ingestion of seafood and inhalation via sea spray and external exposure through swimming, boating and beach occupancy (Fig. 3).
The annual dose from consumption of marine products E marine,k (SvÁyr À1 ) from the ingestion of 8 categories (f) of marine products (piscivorous and non-piscivorous fish, demersal, bottom predator, coastal predator, crustaceans, molluscs and macro-algae) for a given box k is described as follows: where C f,k (Bq kg À1 ) is the activity concentration of the radionuclide in the marine product of type f, CR f,k is the marine food intake rate (kg y À1 ), DC ing (Sv Bq À1 ) is the dose coefficient due to ingestion from marine products given by [18]. The annual dose due to inhalation from sea-spray E inh,k (Sv y À1 ) for a given box k is described as follows: where C air,k (Bq m À3 ) is the nuclide specific concentration in air at a given distance from the coast line, calculated from the water concentration using an empirical relationship, R inh (m 3 h -1 ) is the inhalation rate for an individual, which is assumed 7300 m 3 y -1 . The DC inh (Sv Bq -1 ) is the dose coefficient for inhalation [18] and T inh is the occupancy time (h y -1 ). The external dose from beach material E beach,k (Sv y À1 ) for a given box k is described as follows: where C 1;k (Bq m À2 ) is the nuclide specific surface activity concentration in the shore and beach sediment, DC soil (Sv h -1 per Bq kg -1 ) is the dose coefficient for external exposure [18], and 60 is the areal density of the sediment layer (kg m À2 ), T beach is the occupancy time on the beach (h y -1 ). The external annual dose from swimming E swim,k (Sv y À1 ) for a given box k is described as follows: where DC subm (Sv h À1 per Bq m -3 ) is the dose coefficient for full submerging in water [18], C w,k is the concentration of radionuclide in the surface water layer (Bq m -3 ) and T swim is the occupancy time for swimming (h y À1 ). The external annual dose for boating E boat,k (Sv y À1 ) for a given box k is analogue to the dose for swimming E swim and is described as follows: where T boat is the occupancy time for boating (h y À1 ). As a conservative assumption the dose coefficient is taken as half the dose coefficient for full submerging in water, to consider only exposure from the water surface below the boat. The total annual dose E t (Sv y À1 ) is calculated by summing all individual pathways together. The effective committed dose can be calculated by integration in time for a prescribed time period. The doses for non-human biota are calculated using the ERICA Tool (v1.0) methodology [19].

Case study
In the case study we considered transport and fate of 137 Cs in the northeast Atlantic shelf seas (Fig. 4) released in the period 1950-2020 due to the global deposition after nuclear weapon testing and Chernobyl accident (Fig. 5a) and due to the routine discharges from the Sellafield and La Hague reprocessing plants (Fig. 5b). The box system for the North-Eastern Atlantic includes 108 boxes and covers North Sea, Irish Sea, English Channel, Biscay Bay and adjacent ocean areas. Volume and average depth for each box was calculated using the bathymetry data from [20]. Deep boxes were subdivided on three vertical layers to describe the vertical structure of the radioactivity transport in the upper layer (0-100 m), intermediate layer (100-500 m), and deeper layer (> 500 m). These boxes are marked by blue in the Fig. 4. Water fluxes between boxes were calculated by averaging over 10 years of threedimensional currents calculated from reanalysis [20]. Notice that the water balance should be checked for each box to satisfy the mass conservation. The particular attention should also be paid to an accurate description of advective and diffusive water fluxes through narrows, such as the Northern Passage in the Irish Sea, which control the transfer of contamination to the open ocean.
The simulation of transport and fate of 137 Cs in the North-Eastern Atlantic was carried out for the period 1950-2020. The main sources of 137 Cs as included in this simulation are: (i) global deposition from the weapon testing and from Chernobyl accident [21], (ii) release from the Sellafield and La Hague reprocessing plants [22], (ii) river runoff, and (iv) flow of activity through boundaries. Temporal variations of the sources (a) and (b) are shown in Fig. 5. The flux of 137 Cs from five main rivers (Elbe, Rhine, Seine, Loire, and Garonne) was estimated by using a generic river runoff model [23]. The flows of 137 Cs  activity through open boundaries of the computational domain were estimated by using observations from the MARiS (Marine Information System) database [24].
Results of simulation for period 1950-2020 are given for different boxes in Fig. 6. As seen in Fig. 6a and b, accounting for resuspension processes has significantly improved the agreement between measurement data and model predictions in the shallow Irish Sea. Measurements in the surface water in the modelling domain were used from the MARiS and OSPAR databases [24,25], whereas measurements of 137 Cs concentration in the bottom sediments of the Irish Sea were collected from the publications [26][27][28][29][30][31]. The modelling agrees well with measurements [25] of 137 Cs concentrations both in the demersal fish (flounder, Pleuronectes platessa) and in the coastal predator (cod, Gadus morhua) (Fig. 6c,d). Observations also confirm the importance of the activity pathway from organic sediments to the demersal fish, bottom and coastal predators (Fig. 1b). The geometric mean (GM) of the simulated-to-observed ratios for concentration in the water and sediment are 0.99 and 1.4, respectively. The geometric standard deviation (GSD) for concentration in water is 1.85 for a total number of observations N = 790 in the whole modelling domain, whereas the corresponding value for concentration in the sediment is of 1.5 for a total number of observations N = 29 in the Irish Sea. The GM for the different species of fishes is in the range 1.19-1.27, whereas the corresponding GSD was in the range 1.25-1.34 for the whole modelling domain and N in the range 33-67.
The detailed comparison of simulation results and measurements in the Baltic and Black seas and off the Pacific coast of Japan during 1945-2020 due to the weapon testing and accidents at the Chernobyl and Fukushima Dai-ichi nuclear power plants was given in [32]. It was shown that results of simulations conducted with generic parameters agreed well with measurements of 137 Cs concentrations in the water, bottom sediments, and in fish for the basins with very different marine environment.

Summary
A detailed description of the three-dimensional compartment model POSEIDON-R for the prediction of transport and fate of radionuclides and their daughter products in the marine environment is given. This includes the equations for transfer of radionuclides in the water and bottom sediment compartments along with a dynamical food chain model for radioactivity transfer in the pelagic and benthic food web. Furthermore, details on the numerical method to solve the model equations are also presented. Novel features of the POSEIDON-R model include: a flexible box system that can handle several types of routine and accidental radioactive releases, "local" boxes to describe off-shore releases and a dynamic food web model which describes both pelagic and benthic pathways of radioactivity transfer in biota. These features are essential for use in emergency response software to nuclear accidents, while models based on a standard box modelling methodology [3] lack such properties. Future development on the dynamic food web model should include effects of fish migration, which till now has been described in the model as a random diffusion process only. Comparison of the model predictions with measurements in the northeast Atlantic seas demonstrates good agreement and confirms detailed model validation [32], which was done for three different marine environments (the Baltic and Black seas, and Japan shelf) using the same generic model parameters. We conclude that the POSEIDON-R is robust software suitable for the decision support systems for emergency response to nuclear accidents (e.g. RODOS [4,5]).

Conflicts of interest
The authors declare that they have no conflicts of interest.
side. This work was also supported by IAEA CRP K41017"Behaviour and Effects of Natural and Anthropogenic Radionuclides in the Marine Environment and their use as Tracers for Oceanography Studies" and KIOST Major ProjectPE99615. We are grateful to Mee Kyung Kim for drawing a graphical abstract. This article benefited from the comments and suggestions of two anonymous reviewers.