Multi-compartment kinetic–allometric (MCKA) model of radionuclide bioaccumulation in marine fish

A model of the radionuclide accumulation in fish taking into account the contribution of different tissues and allometry is presented. The basic model assumptions are as follows. (i) A fish organism is represented by several compartments in which radionuclides are homogeneously distributed. (ii) The compartments correspond to three groups of organs or tissues: muscle, bones and organs (kidney, liver, gonads, etc.) differing in metabolic function. (iii) Two input compartments include gills absorbing contamination from water and digestive tract through which contaminated food is absorbed. (iv) The absorbed radionuclide is redistributed between organs or tissues according to their metabolic functions. (v) The elimination of assimilated elements from each group of organs or tissues differs, reflecting differences in specific tissues or organs in which elements were accumulated. (vi) The food and water uptake rates, elimination rate, and growth rate depend on the metabolic rate, which is scaled by fish mass to the $\mathrm{3}/\mathrm{4}$ power. The analytical solutions of the system of model equations describing dynamics of the assimilation and elimination of ¹³⁴Cs, ⁵⁷Co, ⁶⁰Co, ⁵⁴Mn and ⁶⁵Zn, which are preferably accumulated in different tissues, exhibited good agreement with the laboratory experiments. The developed multi-compartment kinetic–allometric model was embedded into the box model POSEIDON-R (Maderich et al., 2018b), which describes transport of radionuclides in water, accumulation in the sediment and transfer of radionuclides through the pelagic and benthic food webs. The POSEIDON-R model was applied for the simulation of the transport and fate of ⁶⁰Co and ⁵⁴Mn routinely released from Forsmark Nuclear Power Plant (NPP) located on the Baltic Sea coast of Sweden and for calculation of ⁹⁰Sr concentration in fish after the accident at Fukushima Dai-ichi NPP. Computed concentrations of radionuclides in fish agree with the measurements much better than calculated using standard whole-body model and target tissue model. The model with the defined generic parameters could be used in different marine environments without calibration based on a posteriori information, which is important for emergency decision support systems.

1 Introduction

Accumulation of radionuclides in marine organisms is a complicated process that is governed by uptake of radionuclides from water, sediment and food as well as by depuration. In turn, these processes depend on the chemical properties of elements, their roles in metabolic processes, the positions of organisms in the food web and marine environmental parameters. In the case of chronic exposure, the radiological assessment models often assumed an equilibrium approach (Carvalho, 2018), in which concentration in the organism relates to the concentration in water using a biological accumulation factor (BAF). However, to describe highly time-dependent transfer processes resulting from accidental releases, dynamic models for the uptake and retention of activity in marine organisms are necessary (Vives i Batlle et al., 2016). According to Takata et al. (2019), the effective half-lives of post-Fukushima Dai-ichi Nuclear Power Plant (FDNPP) accident disequilibrium of ¹³⁷Cs in biota ranged from 100 to 1100 d. The most commonly used bioaccumulation models are the whole-body models, where the organism is represented as a single box in which contamination is evenly distributed (e.g. Fowler and Fisher, 2004; Vives i Batlle et al., 2016). However, the distribution of radionuclides in organisms, and in particular in fish, is non-uniform. For example, the highest concentration of radiocaesium in fish is observed in the muscle, while the highest concentrations of the actinides, plutonium and americium, are measured in specific organs (Coughtrey and Thorne, 1983). Moreover, Vives i Batlle (2012) noted that elimination of activity from organisms occurred with different rates that can be interpreted as elimination from different tissues or organs with different metabolism. In a first approximation, this is used in the “target tissue” approach (Heling et al., 2002; Tateda et al., 2013; Maderich et al., 2014a, b; Bezhenar et al., 2016), where radionuclides are grouped into several classes depending on the type of tissues in which a specific radionuclide accumulates preferentially (target tissue). However, the contribution of other tissues with greater mass than the mass of the target tissue can be commensurate with the contribution of the target tissue to the amount of radioactivity in the body. Distribution of accumulated activities of isotopes Cs, Sr and Co in muscle, bone and liver estimated from previously reported data (Yankovich, 2003; Yankovich et al., 2010) are shown in Fig. 1. The accumulated activity in a given tissue was calculated as a ratio of tissue mass fraction (%) (Yankovich, 2003) to body-to-tissue concentration ratio (Yankovich et al., 2010). As seen in Fig. 1, the accumulated activity of Sr in muscle is not negligible in comparison with accumulated activity in bones, whereas the accumulated activity of Co is redistributed between muscle, bone and liver.

Figure 1 Distribution of accumulated activities of isotopes Cs, Sr and Co in muscle, bone and liver estimated from data (Yankovich, 2003; Yankovich et al., 2010).

A more general approach to the description of the radionuclide accumulation in the tissues of fish is using the physiologically based pharmacokinetic (PBPK) models (Barron et al., 1990; Thomann et al., 1997; Garnier-Laplace et al., 2000; Otero-Muras et al., 2010; Grech et al., 2019). In the PBPK models, the fish organism is represented as three groups of compartments: absorption compartments simulating uptake of contaminants, distribution compartments simulating tissues and organs, and excretion compartments. The exchange of contaminants between compartments is limited by blood flux perfusing compartments. However, these models require a significant number of parameters depending on elements, fish species and marine environments. They must be determined from laboratory experiments (Thomann et al., 1997) or by optimization procedures (Otero-Muras et al., 2010). Note that, with the exception of model Grech et al. (2019), PBPK fish models do not include scaling (allometric) relationships between metabolic rates and organism mass (West et al., 1997; Higley and Bytwerk, 2007; Vives i Batlle et al., 2007; Beresford et al., 2016). Therefore, there is a need to develop a generic model of intermediate complexity between the one-compartment model and the PBPK model taking into account (i) the heterogeneity of the distribution of contamination in fish tissues and (ii) the allometric relationships between metabolic rates and organism mass. Such a model can be used for accidental release simulations without local calibration, which is a complicated task in the circumstances of the accident.

In this paper, a new approach for predicting radionuclide accumulation in fish by taking into account the contributions of different tissues and allometry is presented. The developed multi-compartment kinetic–allometric (MCKA) model is embedded into the box model POSEIDON-R (Lepicard et al., 2004; Maderich et al., 2014a, b; 2018b; Bezhenar et al., 2016), which describes transport of radionuclides in water, accumulation in the sediment, and transfer of radionuclides through the pelagic and benthic food webs. The paper is organized as follows. The MCKA model is described in Sect. 2. The comparison with laboratory experiments is given in Sect. 3. The results of simulation of several radionuclides in the marine environment for regular and accidental releases are described in Sect. 4. The conclusions are presented in Sect. 5.

2 Model

2.1 Model equations

Here, a simple multi-compartmental model to simulate kinetics of radionuclides in the fish is described. The basic assumptions are as follows: (i) a fish organism is represented by several compartments in which radionuclides are homogeneously distributed; (ii) the compartments correspond to three groups of organs or tissues differing in metabolic function: muscle, bones and organs (kidney, liver, gonads, etc.); (iii) two input compartments include gills which absorb contamination from water and digestive tract through which contaminated food is absorbed; (iv) the absorbed radionuclide is redistributed between organs or tissues according to their metabolic functions; (v) the elimination of assimilated elements from each group of organs or tissues differs, reflecting differences in the specific tissues or organs in which elements were accumulated; and (vi) the food and water uptake rates, elimination rate, and growth rate depend on the metabolic rate, which is scaled by fish mass to the $\mathrm{3}/\mathrm{4}$ power following general theory (West et al., 1997). A schematic representation of the model is shown in Fig. 2.

Figure 2 Schematic of the multi-compartment kinetic–allometric model.

The equation for concentration of radionuclide C₁ [Bq kg⁻¹] wet weight (WW) in the gill compartment (i=1) is written as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\text{d}{m}_{\mathrm{1}}{C}_{\mathrm{1}}}{\mathrm{d}t}}=K_{\text{w}}m{C}_{\text{w}}-k_{\mathrm{1}}{m}_{\mathrm{1}}{C}_{\mathrm{1}}-{\mathit{\lambda }}_{\mathrm{1}}{m}_{\mathrm{1}}{C}_{\mathrm{1}}.\end{array}$$

The equation for concentration of unabsorbed radionuclide C₂ [Bq kg⁻¹] (WW) in the digestive tract compartment (i=2) is

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\text{d}{m}_{\mathrm{2}}{C}_{\mathrm{2}}}{\mathrm{d}t}}=K_{\text{f}}m{C}_{\text{f}}-k_{\mathrm{2}}{m}_{\mathrm{2}}{C}_{\mathrm{2}}-{\mathit{\lambda }}_{\mathrm{2}}{m}_{\mathrm{2}}{C}_{\mathrm{2}}.\end{array}$$

The equations for concentrations of radionuclide C_i [Bq kg⁻¹] (WW) in the muscle (i=3), bones (i=4) and organs (i=5) are

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\text{d}{m}_i{C}_i}{\mathrm{d}t}}=k_{\mathrm{2}i}{m}_{\mathrm{2}}{C}_{\mathrm{2}}+k_{\mathrm{1}i}{m}_{\mathrm{1}}{C}_{\mathrm{1}}-{\mathit{\lambda }}_i{m}_i{C}_i.\end{array}$$

Here, t is time; m_i is the mass of the ith tissue; m [kg] is the total mass of fish; C_w [Bq m⁻³] is the concentration of radionuclide in water; C_f [Bq kg⁻¹] (WW) is the concentration of radionuclide in food; K_w [m³(kg d)⁻¹] is a water uptake rate per unit fresh weight of fish and K_f [kg(kg d)⁻¹] is a food uptake rate per unit fresh weight of fish; λ₁ is a loss rate from gills to water [d⁻¹]; k₁, k_1i, k₂, and k_2i [d⁻¹] are transfer rates between tissues; and λ₂ [d⁻¹] is the egestion rate from the digestive tract; λ_i [d⁻¹] is the absorbed radionuclide elimination rate from different tissues or organs ($i=\stackrel{\mathrm{‾}}{\mathrm{3},\mathrm{5}}$).

The activity concentration in the food C_f is expressed by the following equation, summing for a total of n prey types

$$\begin{array}{}\text{(4)}& C_{\text{f}}=\sum _{j=\mathrm{0}}^n{C}_{\text{prey},j}{P}_j{\displaystyle \frac{{\text{drw}}_{\text{pred}}}{{\text{drw}}_{\text{prey},j}}},\end{array}$$

where C_prey,j is the activity concentration in prey of type $\mathrm{0}\le j\le n$, P_j is preference for prey of type j, drw_pred is the dry weight fraction of fish and drw_prey,j is the dry weight fraction of prey of type j. The mean whole-body concentration of activity in the organism C_wb and whole-body activity A_wb [Bq] are calculated as

$$\begin{array}{}\text{(5)}& A_{\text{wb}}=m{C}_{\text{wb}}=\sum _{i=\mathrm{1}}^{\mathrm{5}}{m}_i{C}_i=m\sum _{i=\mathrm{1}}^{\mathrm{5}}{\mathit{\mu }}_i{C}_i,\end{array}$$

where μ_i represents weighting factors $\left({\sum }_{i=\mathrm{1}}^{\mathrm{5}}{\mathit{\mu }}_i=\mathrm{1}\right)$.

Transfer rates k₁ and k₂ are related with tissue transfer rates k_1i and k_2i as

$$\begin{array}{}\text{(6)}& k_{\mathrm{1}}=\sum _{i=\mathrm{3}}^{\mathrm{5}}{k}_{\mathrm{1}i},k_{\mathrm{2}}=\sum _{i=\mathrm{3}}^{\mathrm{5}}{k}_{\mathrm{2}i}.\end{array}$$

Summing Eqs. (1) to (3) yields the equation for total concentration of activity in fish C_wb as

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{C}_{\text{wb}}}{\mathrm{d}t}}=K_{\text{w}}{C}_{\text{w}}+K_{\text{f}}{C}_{\text{f}}-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{1}}{C}_{\mathrm{1}}-{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\mu }}_{\mathrm{2}}{C}_{\mathrm{2}}\\ & -\sum _{i=\mathrm{3}}^{\mathrm{5}}{\mathit{\lambda }}_i{\mathit{\mu }}_i{C}_i-{\mathit{\lambda }}_{\text{g}}{C}_{\text{wb}},\end{array}\end{array}$$

where λ_g is the organism growth rate, defined as

$$\begin{array}{}\text{(8)}& {\mathit{\lambda }}_{\text{g}}={\displaystyle \frac{\mathrm{1}}{m}}{\displaystyle \frac{\mathrm{d}m}{\mathrm{d}t}}.\end{array}$$

The growth dilution can be ignored in the model calculations when λ_g≪λ_i. For short-lived radionuclides, λ_i should be corrected by taking into account the physical decay. The assimilation efficiencies of elements from water AE_w and food AE_f (Pouil et al., 2018) can be introduced, assuming that uptake from water and food is equilibrated by loss to the water from gills and through egestion. The corresponding relations are

$$\begin{array}{}\text{(9)}& {\displaystyle }{\text{AE}}_{\text{w}}={\displaystyle \frac{K_{\text{w}}{C}_{\text{w}}-{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\mu }}_{\mathrm{1}}{C}_{\mathrm{1}}}{K_{\text{w}}{C}_{\text{w}}}},\text{(10)}& {\displaystyle }{\text{AE}}_{\text{f}}={\displaystyle \frac{K_{\text{f}}{C}_{\text{f}}-{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\mu }}_{\mathrm{2}}{C}_{\mathrm{2}}}{K_{\text{f}}{C}_{\text{f}}}}.\end{array}$$

Taking into account the relations Eqs. (9)–(10) for constant AE_w and AE_f, Eq. (7) will be similar to the standard whole-body single-compartment equation

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}{C}_{\text{wb}}}{\mathrm{d}t}}={\text{AE}}_{\text{w}}{K}_{\text{w}}{C}_{\text{w}}+{\text{AE}}_{\text{f}}{K}_{\text{f}}{C}_{\text{f}}-({\mathit{\lambda }}_{\text{wb}}+{\mathit{\lambda }}_{\text{g}})C_{\text{wb}}\end{array}$$

if elimination terms in Eq. (7) are replaced by a single term λ_wbC_wb, assuming that λ_wb is a single whole-body elimination rate (Fowler and Fisher, 2004).

Table 1 Parameters in allometric relations, standard deviation (SD) of parameters and number of measurements N.

The food and water uptake rates, elimination rate, and growth rate depend on the metabolic rate, which in turn is known to scale by the organism mass. Here, we employed quarter-power scaling for uptake, elimination and growth rates derived from general theory (West et al., 1997). This theory predicts for all organisms the $\mathrm{3}/\mathrm{4}$ power law for metabolic rates. It describes transport of essential materials through space-filling fractal networks of branching tubes assuming that the energy dissipation is minimized and that the terminal branch of the network is size invariant. The scaling relations are

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& K_{\text{w}}\left(m\right)={\mathit{\alpha }}_{\text{w}}{m}^{-\mathrm{1}/\mathrm{4}},K_{\text{f}}\left(m\right)={\mathit{\alpha }}_{\text{f}}{m}^{-\mathrm{1}/\mathrm{4}},{\mathit{\lambda }}_{\text{g}}\left(m\right)\\ & ={\mathit{\alpha }}_{\text{g}}{m}^{-\mathrm{1}/\mathrm{4}},{\mathit{\lambda }}_i\left(m\right)={\mathit{\alpha }}_i{m}^{-\mathrm{1}/\mathrm{4}},\end{array}\end{array}$$

where ${\mathit{\alpha }}_{\text{w}},{\mathit{\alpha }}_{\text{f}},{\mathit{\alpha }}_{\text{g}},{\mathit{\alpha }}_i$ ($i=\stackrel{\mathrm{‾}}{\mathrm{1},\mathrm{5}}$) are constants. These parameters can also depend on temperature, salinity and fish age (e.g. Belharet et al., 2019; Heling and Bezhenar, 2009). Notice that a number of laboratory experiments (Thomas and Fisher, 2010) showed that temperature exerts no major influence on uptake and elimination, whereas the effect of salinity varies for elements (Heling and Bezhenar, 2009; Jeffree et al., 2017). Here, we did not analyse these factors requiring separate consideration. The values of constants ${\mathit{\alpha }}_{\text{w}},{\mathit{\alpha }}_{\text{f}},{\mathit{\alpha }}_{\text{g}},{\mathit{\alpha }}_i$ ($i=\stackrel{\mathrm{‾}}{\mathrm{1},\mathrm{5}}$) estimated from laboratory experiments and marine data are provided in Table 1, whereas values of weighting factors μ₁=0.01, μ₂=0.01, μ₃=0.78, μ₄=0.12, and μ₅=0.08 were chosen using estimates from Yankovich (2003).

2.2 Kinetics in equilibrium state

The model parameters can be estimated using measurement data and applying the kinetic equations under equilibrium conditions. Equations (1)–(3) rewritten for radionuclide concentrations in the equilibrium state are

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{1}}{C}_{\mathrm{1}}={\displaystyle \frac{K_{\text{w}}{C}_{\text{w}}}{k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}+{\mathit{\lambda }}_{\text{g}}}},\text{(14)}& {\displaystyle }{\mathit{\mu }}_{\mathrm{2}}{C}_{\mathrm{2}}={\displaystyle \frac{K_{\text{f}}{C}_{\text{f}}}{k_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\text{g}}}},\text{(15)}& {\displaystyle }{k}_{\mathrm{2}i}{\displaystyle \frac{K_{\text{f}}{C}_{\text{f}}}{k_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\text{g}}}}+k_{\mathrm{1}i}{\displaystyle \frac{K_{\text{w}}{C}_{\text{w}}}{k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}+{\mathit{\lambda }}_{\text{g}}}}={\mathit{\lambda }}_i{\mathit{\mu }}_i{C}_i+{\mathit{\lambda }}_{\text{g}}{\mathit{\mu }}_i{C}_i.\end{array}$$

Then, using Eqs. (13)–(14), the assimilation efficiencies of elements from water AE_w and food AE_f are rewritten as

$$\begin{array}{}\text{(16)}& {\text{AE}}_{\text{w}}={\displaystyle \frac{k_{\mathrm{1}}+{\mathit{\lambda }}_{\text{g}}}{k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}+{\mathit{\lambda }}_{\text{g}}}},{\text{AE}}_{\text{f}}={\displaystyle \frac{k_{\mathrm{2}}+{\mathit{\lambda }}_{\text{g}}}{k_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{2}}+{\mathit{\lambda }}_{\text{g}}}}.\end{array}$$

When λ₁≫λ_g and λ₂≫λ_g, we determine from Eqs. (16) and (12), approximately, that

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& k_{\mathrm{1}}={\displaystyle \frac{{\text{AE}}_{\text{w}}{\mathit{\lambda }}_{\mathrm{1}}}{\mathrm{1}-{\text{AE}}_{\text{w}}}}={\displaystyle \frac{{\text{AE}}_{\text{w}}{\mathit{\alpha }}_{\mathrm{1}}}{(\mathrm{1}-{\text{AE}}_{\text{w}})m^{\mathrm{1}/\mathrm{4}}}},k_{\mathrm{2}}={\displaystyle \frac{{\text{AE}}_{\text{f}}{\mathit{\lambda }}_{\mathrm{2}}}{\mathrm{1}-{\text{AE}}_{\text{f}}}}\\ & ={\displaystyle \frac{{\text{AE}}_{\text{f}}{\mathit{\alpha }}_{\mathrm{2}}}{(\mathrm{1}-{\text{AE}}_{\text{f}})m^{\mathrm{1}/\mathrm{4}}}}.\end{array}\end{array}$$

This assumption also imposes a requirement on the modelling of radionuclides with decay constant λ≪λ_i. Eq. (17) relations are used to relate kinetic coefficients of the model with experimentally determined parameters AE_w, AE_f, λ₁ and λ₂. Eq. (15) can be rewritten as

$$\begin{array}{}\text{(18)}& {\text{AE}}_{\text{f}i}{\mathit{\alpha }}_{\text{f}}{\text{BAF}}_{\text{food}}+{\text{AE}}_{\text{w}i}{\mathit{\alpha }}_{\text{w}}={\mathit{\mu }}_i({\mathit{\alpha }}_i+{\mathit{\alpha }}_{\text{g}}){\displaystyle \frac{{\text{BAF}}_{\text{wb}}}{{\text{CR}}_i}},\end{array}$$

where AE_wi and AE_fi are assimilation efficiencies for tissue i,

$$\begin{array}{}\text{(19)}& \sum _{i=\mathrm{3}}^{\mathrm{5}}{\mathit{\mu }}_i{\text{AE}}_{\text{w}i}\approx {\text{AE}}_{\text{w}},\sum _{i=\mathrm{3}}^{\mathrm{5}}{\mathit{\mu }}_i{\text{AE}}_{\text{f}i}\approx {\text{AE}}_{\text{f}}.\end{array}$$

The assimilation efficiencies are expressed through kinetic coefficients as

$$\begin{array}{}\text{(20)}& {\mathit{\mu }}_i{\text{AE}}_{\text{w}i}={\displaystyle \frac{k_{\mathrm{1}i}}{k_{\mathrm{1}}}}{\text{AE}}_{\text{w}},{\mathit{\mu }}_i{\text{AE}}_{\text{f}i}={\displaystyle \frac{k_{\mathrm{2}i}}{k_{\mathrm{2}}}}{\text{AE}}_{\text{f}}.\end{array}$$

We define the bioconcentration factor (BCF) as a ratio of whole body of fish to water concentrations with no dietary intake, bioaccumulation factor (BAF) as ratio of whole body of fish to water concentrations with dietary intake, and body-to-tissue concentration ratio (CR_i) as the ratio of whole body to ith tissue concentrations, whereas the ratio of food to water concentrations is indicated as BAF_food. The parameters BAF, CR_i and BAF_food are described as

$$\begin{array}{}\text{(21)}& \text{BAF}={\displaystyle \frac{C_{\text{wb}}}{C_{\text{w}}}},{\text{CR}}_i={\displaystyle \frac{C_{\text{wb}}}{C_i}},{\text{BAF}}_{\text{food}}={\displaystyle \frac{C_{\text{food}}}{C_{\text{w}}}}.\end{array}$$

Values of BAF for different radionuclides are named as CF (Concentration Factor) in IAEA (2004). Yankovich et al. (2010) provide CR_i based on aggregate experimental data for marine fish. We assume that the kinetics of assimilation in fish tissues are similar for radionuclides absorbed from water and food; i.e.

$$\begin{array}{}\text{(22)}& {\displaystyle \frac{{\text{AE}}_{\text{w}i}}{{\text{AE}}_{\text{w}}}}={\displaystyle \frac{{\text{AE}}_{\text{f}i}}{{\text{AE}}_{\text{f}}}}.\end{array}$$

Notice that assimilation for some elements can be considered route dependent (Reinfelder et al., 1999), and so Eq. (22) is only a first approximation. Inserting Eq. (22) into Eq. (18) yields

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& {\displaystyle \frac{{\text{AE}}_{\text{f}i}}{{\text{AE}}_{\text{f}}}}\left({\text{AE}}_{\text{f}}{\mathit{\alpha }}_f{\text{BAF}}_{\text{food}}+{\text{AE}}_{\text{w}}{\mathit{\alpha }}_{\text{w}}\right)\\ & =\text{BAF}{\displaystyle \frac{{\mathit{\mu }}_i({\mathit{\alpha }}_i+{\mathit{\alpha }}_{\text{g}})}{{\text{CR}}_i}},\end{array}\end{array}$$

Summing Eq. (23) for $i=\mathrm{3},\mathrm{4},\mathrm{5}$ yields

$$\begin{array}{}\text{(24)}& {\text{AE}}_{\text{f}}{\mathit{\alpha }}_{\text{f}}{\text{BAF}}_{\text{food}}+{\text{AE}}_{\text{w}}{\mathit{\alpha }}_{\text{w}}=\text{BAF}\sum _{i=\mathrm{3}}^{\mathrm{5}}{\displaystyle \frac{{\mathit{\mu }}_i({\mathit{\alpha }}_i+{\mathit{\alpha }}_{\text{g}})}{{\text{CR}}_i}}.\end{array}$$

Using Eqs. (23) and (24), the ratio AE_fi $/$ AE_f can be written as

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{{\text{AE}}_{\text{f}i}}{{\text{AE}}_{\text{f}}}}={\mathit{\mu }}_i({\mathit{\alpha }}_i+{\mathit{\alpha }}_{\text{g}})/\left({\text{CR}}_i\sum _{i=\mathrm{3}}^{\mathrm{5}}{\displaystyle \frac{{\mathit{\mu }}_i({\mathit{\alpha }}_i+{\mathit{\alpha }}_{\text{g}})}{{\text{CR}}_i}}\right).\end{array}$$

The values of kinetic coefficients k₁ and k₂ and assimilation efficiencies for tissues AE_fi were calculated from Eqs. (17) and (25) using assimilation efficiencies from experimental data (Pouil et al., 2018). The values of AE_fi for several radionuclides are given in Table 2.

Table 2 The food assimilation efficiency AE_f (Pouil et al., 2018), assimilation efficiency of elements from water AE_w, tissue assimilation efficiencies AE_fi and trophic transfer factor (TTF) for several elements.

^a The values were estimated from experiments by Mathews and Fisher (2008). ^b The value was estimated using BAF from IAEA (2004). ^c Heling and Bezhenar (2009) for seawater environment.

Equation (24) can be rearranged to express the ratio of BAF_wb to BAF_food by taking into account dominance of dietary intake over water intake (Mathews and Fisher, 2009). This ratio is the trophic transfer factor (TTF), written as

$$\begin{array}{}\text{(26)}& \text{TTF}={\text{AE}}_{\text{f}}{\mathit{\alpha }}_{\text{f}}/\sum _{i=\mathrm{3}}^{\mathrm{5}}{\displaystyle \frac{{\mathit{\mu }}_i({\mathit{\alpha }}_i+{\mathit{\alpha }}_{\text{g}})}{{\text{CR}}_i}}.\end{array}$$

TTF > 1 indicates possible biomagnification, and TTF < 1 indicates that biodiminution is likely. As follows from Eq. (26), the TTF value does not depend on the mass of fish. The TTF values calculated by using Eq. (26) are given in Table 2. Among the considered elements, only caesium (TTF > 1) may be biomagnified in the food chain, in agreement with Kasamatsu and Ishikawa (1997), where it was found that the BAF of ¹³⁷Cs increased with increasing trophic level. The concentrations of other radionuclides in Table 2 decreased with the increase of trophic level (TTF < 1), which is consistent with the findings presented by Cardwell et al. (2013), where an inverse relationship was obtained between trophic levels and the concentration of inorganic metals in water chains.

Notice that values of AE_f, AE_w and BAF do not depend on the fish mass. The literature data reveal diverse relationships between fish mass and both radionuclide bioconcentration (BCF) and bioaccumulation (BAF) factors. In particular, data of laboratory experiments (Mathews et al., 2008) showed that there is no significant relationship between bioconcentration factor BCF and fish size for most studied aqueous metals. The BAF in larger and older fish of the same species can differ from smaller and younger fish due to the change of habitat and diet with age (e.g. Kasamatsu and Ishikawa, 1997; Ishikawa et al., 1995; Kim et al., 2019); however, in this study we did not consider the change of prey preference along the fish growth.

3 Comparison with laboratory experiments

3.1 Depuration of radionuclides after pulse-like feeding

Retention of absorbed elements in fish after single feeding was often used to estimate AE_f and depuration rate (e.g. Jeffree et al., 2006; Mathews and Fisher, 2008; Mathews et al., 2008; Pouil et al., 2017). According to Goldstein and Elwood (1971), single feeding can be approximated by a delta function δ(t) at t=0 as

$$\begin{array}{}\text{(27)}& K_{\text{f}}m{C}_{\text{f}}=A_{\text{f}}\mathit{\delta }\left(t\right),\end{array}$$

where A_f is the total amount of ingested activity. The solutions of Eqs (1)–(3) for activities A₂=m₂C₂ and A_i=m_iC_i and for initial conditions $A_{\mathrm{2}}=A_i=\mathrm{0}$ are

$$\begin{array}{}\text{(28)}& {\displaystyle \frac{A_{\mathrm{2}}}{A_{\mathrm{f}}}}={\displaystyle \frac{m_{\mathrm{2}}{C}_{\mathrm{2}}}{A_{\mathrm{f}}}}=\exp (-(k_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{2}}\left)t\right),\text{(29)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{A_i}{A_{\mathrm{f}}}}={\displaystyle \frac{m_i{C}_i}{A_{\mathrm{f}}}}={\displaystyle \frac{k_{\mathrm{2}i}}{k_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{2}}-{\mathit{\lambda }}_i}}\left[\exp \right(-{\mathit{\lambda }}_{i}t)\\ & -\exp (-(k_{\mathrm{2}}+{\mathit{\lambda }}_{\mathrm{2}}\left)t\right)].\end{array}\end{array}$$

As follows from these solutions, the decay of activity in the fish organisms includes a fast component with decay constant representing transfer of activity to the fish body and unabsorbed element egestion from the digestive tract, along with a slow component which is governed by elimination constants for i tissues. These solutions are generalized solutions of the equations of the sequentially linked two-compartment model by Goldstein and Elwood (1971), whereas in the D-DAT model (Vives i Batlle et al., 2008) an organism was represented by two boxes with parallel kinetics, also describing “slow” and “fast” exchange processes.

The analytical solutions (Eqs. 28–29) can be compared with laboratory experiments in which depuration of metals from the fish after single feeding was studied. In the experiment by Mathews and Fisher (2008), the retention of several radioisotopes in juvenile sea bream (Sparus auratus) was considered. The average wet weight of the fish was 0.0001 kg. These fish were fed radiolabelled Artemia salina nauplii. The fish were allowed to feed for 1 h, after which metal retention was observed in clean water over a 15 d period. The solutions with parameters corresponding to the fish mass and metal AE (Pouil et al., 2018) were compared with experimental data for ¹³⁴Cs, ⁶⁰Co, ⁵⁴Mn, and ⁶⁵Zn in Fig. 3. As seen in Fig. 3, both model and experiments showed two phases (fast and slow) of radionuclide elimination. Most of the activity is contained in muscle; however, the first 5 d after feeding the concentrations of ⁶⁰Co and ⁵⁴Mn in organs are much greater than in the muscle.

Figure 3 Retention of radionuclides in whole body and tissues of juvenile sea bream (Sparus auratus). The simulations are compared with whole-body measurements by Mathews and Fisher (2008).

The solutions were also compared with laboratory experiments for predator fish (Mathews et al., 2008). In these experiments, the retention of several radioisotopes in sea bream (Sparus auratus), turbot (Psetta maxima) and spotted dogfish (Scyliorhinus canicula) was studied. Immature S. auratus (wet weight 0.012 kg), P. maxima (wet weight 0.027 kg) and S. canicula (wet weight 0.008 kg) were fed radiolabelled prey fish (juvenile S. auratus). After this feeding, fish were fed unlabelled prey fish for 3 weeks. In Fig. 4, the solutions Eqs. (28)–(29) were compared with the laboratory experiments in which prey fish were labelled by ¹³⁴Cs, ⁵⁷Co, ⁵⁴Mn and ⁶⁵Zn. The analytical solutions (Eqs. 28–29) and experiment agreed, demonstrating general dependence of the depuration process on fish mass. Differences between experimental data for different species may be due to differences in anatomy and physiology, as discussed by Jeffree et al. (2006) for P. maxima and S. canicula. The model, unlike the situation for prey fish (Fig. 3), underestimates the total concentrations of ⁵⁷Co and ⁵⁴Mn in comparison with experiments, which is probably due to the neglect of other factors, except body weight, for the bioaccumulation kinetics. Parameters of the MCKA model for fish from experiments (Mathews and Fisher, 2008; Mathews et al., 2008) are given in Table S1, whereas Tables S2–S4 show dependence on radionuclides of the transfer rates k_2,i in different fishes.

Figure 4 Retention of radionuclides in whole bodies of predator fish: sea bream (Sparus auratus), turbot (Psetta maxima) and spotted dogfish (Scyliorhinus canicula). The simulations are compared with whole-body measurements by Mathews et al. (2008).

3.2 Bioconcentration of dissolved radionuclides from sea water

Uptake and absorption in fish of elements from water were studied in several laboratory experiments (e.g. Jeffree et al., 2006; Mathews and Fisher, 2008; Mathews et al., 2008). The modelling of the absorption of elements can be used to estimate an assimilation efficiency AE_w. An analytical solution of Eqs. (1) and (3) with initial conditions C_i=0 at t=0 is written as

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{{\mathit{\mu }}_{\mathrm{1}}{C}_{\mathrm{1}}}{C_{\text{w}}}}={\displaystyle \frac{K_{\text{w}}}{k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}}}\left[\mathrm{1}-\exp (-(k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}\left)t\right)\right],\text{(31)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{{\mathit{\mu }}_i{C}_i}{C_{\text{w}}}}={\displaystyle \frac{k_{\mathrm{1}i}{K}_{\text{w}}}{(k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}){\mathit{\lambda }}_i}}[\mathrm{1}-{\displaystyle \frac{k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}}{k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}-{\mathit{\lambda }}_i}}\exp (-{\mathit{\lambda }}_{i}t)\\ & +{\displaystyle \frac{{\mathit{\lambda }}_i}{k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}-{\mathit{\lambda }}_i}}\exp (-(k_{\mathrm{1}}+{\mathit{\lambda }}_{\mathrm{1}}\left)t\right)],(i=\stackrel{\mathrm{‾}}{\mathrm{3},\mathrm{5}}).\end{array}\end{array}$$

Figure 5 Simulated BCF in marine fish during the period of exposition in water. The simulations are compared with isotope measurements by Mathews et al. (2008) in juvenile sea bream (Sparus auratus) and measurements by Jeffree et al. (2006) in turbot (Psetta maxima) and spotted dogfish (Scyliorhinus canicula).

These solutions were compared with laboratory experiments for prey fish (Mathews et al., 2008) and for predator fish (Jeffree et al., 2006). In the experiment by Mathews et al. (2008), the uptake of several radioisotopes by juvenile S. auratus (wet weight 0.0002 kg) was studied during 25 d of exposure, whereas in experiments by Jeffree et al. (2006) immature P. maxima (wet weight 0.0061 kg) and S. canicula (wet weight 0.0067 kg) were used for study during a 15 d period. Parameters of the MCKA model for these fishes are given in Table S5, whereas Tables S6–S8 show dependence on radionuclides of the transfer rates k_1,i in different fishes. The comparison of the analytical solutions (Eqs. 30–31) with experimental data with respect to bioconcentration factor (BCF = C${}_{\text{wb}}/C_{\text{w}}$ [L kg⁻¹]) is presented in Fig. 5. The values of AE_w were selected to approximate the experiment for small prey fish (Mathews et al., 2008). They differ for different metals. For ¹³⁴Cs, the value of AE_w was 0.001, whereas for ⁵⁷Co and ⁶⁰Co, ⁵⁴Mn, and ⁶⁵Zn, these values were 0.0025, 0.0045, 0.0065, respectively. This contrasted with AE_f, which is larger for ¹³⁴Cs than for ⁵⁷Co, ⁶⁰Co, ⁵⁴Mn, and ⁶⁵Zn (Pouil et al., 2017). The estimated values of AE_w are presented in Table 2. The above-mentioned values of AE_w were used to calculate BCF for larger predator fish in experiments by Jeffree et al. (2006). As seen in Fig. 5, values of AE_w for considered elements are of the order 10⁻³, which is in agreement with most models. However, comparison with larger fish highlighted some differences between species of fish, as discussed by Jeffree et al. (2006), and differences between model and experiment for a constant value of AE_w. At the same time, it is known that dietary intake of metals dominates over water intake (Mathews and Fisher, 2009). Therefore, deviations in values of AE_w would not be significantly affected by the full uptake of elements from the marine environment.

The parameter α₁ can be estimated from the relations Eqs. (13) and (16) according to the experimental data for equilibrium conditions. An average value of gill BCF₁ is approximately 10 L kg⁻¹ for ⁵⁸Co, ⁵⁴Mn, ¹³⁴Cs and ⁶⁵Zn (Jefferies and Heweit,1971; Pentreath, 1973). Then for AE_w=0.001, we obtained α₁=800. With the selected value of α₁, the process of adaptation of the gill tissue to changes in the concentration of radioactivity in water is much faster than for other tissues of the fish. In addition, as follows from the analytical solutions (Eqs. 30–31) at λ₁≫λ_i, the contribution of gill contamination to the whole-body contamination is small.

3.3 Simplification of the model based on the results of analytical solutions and laboratory experiments

Comparison of the model against laboratory experiments on the retention of absorbed elements in fish after single feeding demonstrated the need to include the kinetic characteristics of the digestive tract in the model when highly non-equilibrium transfer dynamics are expected. However, for modelling of food uptake in marine environment with multiple feedings the simple equilibrium assumption Eq. (14) can be used. At the same time, the analytical solution describing the bioconcentration due to the uptake and absorption of elements from water in fish, as well as the results of the laboratory experiment, shows that the contribution to the gills (Eq. 13) is negligible. Therefore, for modelling of uptake from water, the equilibrium assumption can also be used as it is shown above. The corresponding simplified equations for muscle, bone and organs can be rewritten as

$$\begin{array}{}\text{(32)}& \begin{array}{rl}& {\mathit{\mu }}_i{\displaystyle \frac{\mathrm{d}{C}_i}{\mathrm{d}t}}={\text{AE}}_{\text{f}i}{K}_{\text{f}}{C}_{\text{f}}+{\text{AE}}_{\text{w}i}{K}_{\text{w}}{C}_{\text{w}}-{\mathit{\lambda }}_i{\mathit{\mu }}_i{C}_i\\ & -{\mathit{\lambda }}_{\text{g}}{\mathit{\mu }}_i{C}_i.\end{array}\end{array}$$

The uncertainty in calculations using Eq. (32) arises due to (i) limited experimental data to define allometry constants α_f, α_w, α₃, α₄, α₅; (ii) large intervals in the values of known AE_f coefficients; (iii) unknown AE_f values for many radionuclides; (iv) a lack of experimental data about AE_w; and (v) limited experimental data for whole-body to tissue concentration ratios, CR_i, in the marine fish. Therefore, a sensitivity analysis is necessary to estimate uncertainty of the simulation results. We estimated the effects of variations in above parameters in the Eq. (24) on the value of BAF in equilibrium state. The simple local sensitivity analysis and the one-at-a-time method (Pianosi et al., 2016) were used. The sensitivity of model output was estimated using a sensitivity index (SI) calculated following Hamby (1994) as

$$\begin{array}{}\text{(33)}& \text{SI}={\displaystyle \frac{D_{\text{max}}-D_{\text{min}}}{D_{\text{max}}}},\end{array}$$

where D_max and D_min are the outputs corresponded to maximal and minimal input parameter values, respectively. Similarly to Bezhenar et al. (2016), the range for every parameter was defined as follows: minimum value was set to half the reference value and maximum value was set to twice the reference value. Calculated SI for three radionuclides, which are preferably accumulated in different tissues (¹³⁷Cs, ⁹⁰Sr and ⁶⁰Co), are given in Fig. S1.

The results of the sensitivity study suggest that model results are most sensitive to variations of AE_f and α_f for ¹³⁷Cs and ⁶⁰Co, whereas for ⁶⁰Co they are almost not sensitive to variations of AE_w and α_w. Note that model results are also sensitive to the variations of parameters related to tissues where the radionuclide is mainly accumulated: α₃ and CR₃ for ¹³⁷Cs and α₅ and CR₅ for ⁶⁰Co. For ⁹⁰Sr, the model results are moderately sensitive to the variations of most above-considered parameters.

4 Model applications

4.1 Modified POSEIDON-R box model

In order to predict the accumulation of radionuclides in fish in the marine environment using the MCKA model described above, it is necessary to calculate changes in concentration in water and in bottom sediments and to calculate the transport of radionuclides through food chains. The POSEIDON-R box model (Lepicard et al., 2004; Maderich et al., 2014a, b; 2018b; Bezhenar et al., 2016) can be used to simulate the marine environment as a system of 3D boxes for the water column, bottom sediment and food web. The water column box is vertically subdivided into layers. The suspended matter settles in the water column. The bottom sediment box is divided into three layers (Fig. S2). The downward burial processes operate in all three sediment layers. Maderich et al. (2018b) described the POSEIDON-R model in detail.

A food web model that includes pelagic and benthic food chains is implemented within the POSEIDON-R box model (Bezhenar et al., 2016). In the food web model, marine organisms are grouped into classes according to trophic level and species type (Fig. S3). The food chains differ between the pelagic zone and the benthic zone. Pelagic organisms comprise primary producers (phytoplankton) and consumers (zooplankton, non-piscivorous (forage) fish and piscivorous fish). In the benthic food chain, radionuclides are transferred from algae and contaminated bottom sediments to deposit-feeding invertebrates, demersal fish and benthic predators. Bottom sediments include both organic and inorganic components. Radioactivity is assumed to be assimilated by benthic organisms from the organic components of the bottom deposits. Other food web components are crustaceans (detritus feeders), molluscs (filter feeders) and coastal predators, which feed throughout the water column in shallow coastal waters. All organisms take in radionuclides both via the food web and directly from the water. Table S9 in the Supplement contains food preferences for organisms in the food web which are used in the model. Details of the transfer of radiocaesium through the marine food web are presented by Bezhenar et al. (2016) and Maderich et al. (2018b).

The POSEIDON-R model can handle different types of radioactive releases: including atmospheric fallout and point sources associated with routine releases from nuclear facilities located directly on the coast or point sources associated with accidental releases (Lepicard et al., 2004). For coastal discharges occurring in the large (“regional”) boxes, “coastal” release boxes are nested into the regional box system. The intermediary boxes between coastal and regional boxes are called “inner” boxes.

Figure 6 Box system around the Forsmark NPP. The numbers denote regional boxes in the box system of the Baltic Sea (Bezhenar et al., 2016). An additional “Inner box” is separated from regional box 68 by the blue line. The coastal box (red rectangle) surrounds the area where cooling water from NPP is discharged.

4.2 Release of radionuclides during normal operation of the Forsmark Nuclear Power Plant

This section presents the simulation results of ⁶⁰Co and ⁵⁴Mn routine release into the marine environment from Forsmark NPP, located on the Baltic Sea coast of Sweden. The POSEIDON-R model with embedded food web MCKA model, one-compartment fish model and bioaccumulation factor (equilibrium) model was customized for the Baltic Sea, as described by Bezhenar et al. (2016) and Maderich et al. (2018a). The nested boxes (inner and coastal boxes) inside the regional box no. 68 in the Baltic Sea box system were added to resolve the radionuclide concentration in the near field (Fig. 6). Parameters of the inner and coastal boxes are based on data from Aquilonius (2010). The main parameters of boxes when zooming in to the NPP are presented in Table S10 in the Supplement. The simulation results for ⁶⁰Co and ⁵⁴Mn were compared with measurements for the smallest coastal box, where measurement data for bottom sediments and fish were available (Forsmark, 2014). The measurement data were compared with simulations for two species of fish: herring (Clupea harengus membras) as a non-piscivorous fish and pike (Esox lucius) as a coastal predator. There is no information on the mass of fish caught in the vicinity of Forsmark NPP. Therefore, we used estimates of the masses typical for prey and predatory fish, which are given in Table S11. In the one-compartment model, two parameters must be prescribed: assimilation efficiency and biological half-life $T_{\mathrm{0.5}}=\ln \left(\mathrm{2}\right){\mathit{\lambda }}_{\text{wb}}^{-\mathrm{1}}$. Assimilation efficiency (see Table 2) was obtained from Pouil et al. (2018). Baudin et al. (1997) used T_0.5=21 d for ⁶⁰Co in the one-compartment model. The average value for T_0.5 in predatory marine fish is 40 d (Beresford et al., 2015). Therefore, for ⁶⁰Co, we used T_0.5=20 d for prey fish and T_0.5=40 d for predatory fish. There are very limited data for T_0.5 values in marine fish for ⁵⁴Mn. According to Beresford et al. (2015), T_0.5 is in the range between 20 and 40 d. Therefore, we used the same values of T_0.5 for ⁵⁴Mn, as for ⁶⁰Co. According to Jeffree et al. (2017), the uptake and depuration kinetics of ⁶⁰Co and ⁵⁴Mn for fish species in marine, brackish and freshwater environments are very similar. Therefore, we can apply the model parameters defined for marine environment (see Table S11) to reconstruct the herring and pike contamination by the above-mentioned radionuclides in the area near the Forsmark NPP with low salinity (3–5 PSU, denoting practical salinity unit).

The release rates of ⁶⁰Co from the Forsmark NPP (Forsmark, 2014) are plotted in Fig. 7.

Figure 7 Release rates of ⁶⁰Co and ⁵⁴Mn according to measurements (Forsmark, 2014).

As seen in Fig. 8a, the results of simulation for the concentration of ⁶⁰Co in the bottom sediments are in good agreement with the measurements (Forsmark, 2014) in the coastal box for the wide range of employed values of sediment distribution coefficient K_d: from $K_{\text{d}}=\mathrm{3}\times {\mathrm{10}}^{\mathrm{5}}÷\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ L kg⁻¹ for margin seas to $K_{\text{d}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{7}}$ L kg⁻¹ for open ocean (IAEA, 2004). The benthic food web (Bezhenar et al., 2016a), which describes the transfer of radioactivity from bottom sediments to deposit-feeding invertebrates and finally to fish, is quite important in this range of K_d values. The results from the POSEIDON-R calculations obtained with the MCKA model are compared (Fig. 8b, c) with measurements and results of calculations obtained with a one-compartment model and with an equilibrium approach using a standard BAF value (IAEA, 2004). Whereas one-compartment and equilibrium models underestimated the concentration of ⁶⁰Co in fish, the MCKA model using generic parameters yields better agreement with measurements for both non-piscivorous fish (Fig. 8b) and coastal predator feeding by pelagic and benthic organisms in the coastal area (Fig. 8c).

Similarly, the behaviour of ⁵⁴Mn in the marine environment near the Forsmark NPP is modelled. The release rate of ⁵⁴Mn from the NPP is also plotted in Fig. 7 using data from Forsmark (2014). Comparison of the simulated concentration of ⁵⁴Mn in bottom sediments with measurements (Forsmark, 2014) for the Forsmark coastal box is given in Fig. 8d. Good agreement was obtained when a standard value of $K_{\text{d}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{6}}$ L kg⁻¹ for ⁵⁴Mn in margin seas (IAEA, 2004) was used. This means that, as in the case of ⁶⁰Co, a significant fraction of radionuclide is deposited at the bottom, and the benthic food web should be considered. Similarly to the ⁶⁰Co case, obtained results of simulation are also compared with results obtained using the one-compartment model and equilibrium approach. Again, the MCKA model yields the best agreement with measurements for both non-piscivorous fish (Fig. 8e) and coastal predators (Fig. 8f); however, this agreement is slightly worse than in the ⁶⁰Co case. Notice that the BAF in the equilibrium approach can be locally estimated using a posteriori data. However, the MCKA model provided good agreement with measurements using only a priori information, which is important in the case of accidents, as considered in the next section.

Figure 8 Comparison between calculated and measured ⁶⁰Co (Forsmark, 2014) concentrations in bottom sediment (a), non-piscivorous fish (herring) (b) and coastal predator fish (pike) (c) and ⁵⁴Mn concentrations in bottom sediment (d), non-piscivorous fish (herring) (e) and coastal predator fish (pike) (f) for the coastal box.

4.3 Accumulation of ⁹⁰Sr in the fish after the Fukushima Dai-ichi accident

Following caesium, ⁹⁰Sr is the second most important radiologically long-lived radionuclide released as a result of the FDNPP accident. It is highly soluble in water and exhibits a relatively high ability for assimilation by marine organisms due to similar chemical properties with calcium. The atmospheric deposition of ⁹⁰Sr is usually not taken into account due to its low volatility. Most of the ⁹⁰Sr released from the FDNPP was directly released to the ocean, with estimates of total inventory in the range from 0.04 to 1.0 PBq (Buesseler et al., 2017). Here, we extended the simulation by Maderich et al. (2014b) of transfer and fate of ⁹⁰Sr resulting from the FDNPP accident using the POSEIDON-R model complemented by the food web model (Bezhenar et al., 2016) and MCKA fish model.

The POSEIDON-R model was customized for the north-western Pacific and adjacent seas (the East China Sea, the Yellow Sea, and the East Sea or Sea of Japan) as in (Maderich et al., 2018a). A total of 188 boxes covered this region. The boxes around the FDNPP are shown in Fig. 9 with an additional 4 km×4 km coastal box in the vicinity of FDNPP. The coastal box was included for the model validation with ⁹⁰Sr measurements conducted in this specific area. Details of customization are given in Maderich et al. (2014a, b; 2018a).

The historical contamination due to global atmospheric deposition in the period from 1945–2010 was simulated according to Maderich et al. (2014b) with data from the Marine Radioactivity Information System database (MARIS, 2020). The value of the accidental release was estimated as 160 TBq (16 TBq d⁻¹ during 10 d) that was consistent with the range reported by Buesseler et al. (2017). In the post-accidental period, the continuous leakage of ⁹⁰Sr due to groundwater transport of radioactivity from the NPP site was monitored (Castrillejo et al., 2016). Therefore, in the simulation the conservative scenario was used for release of ⁹⁰Sr in the post-accidental period (Fig. 9); the release of ⁹⁰Sr was assumed equal to ¹³⁷Cs release (Maderich et al. (2018a).

Figure 9 (a) Boxes along the eastern coast of Japan with fine resolution in the area around the FDNPP (Maderich et al., 2018a). Deep boxes with three vertical layers in the water column are coloured by blue, boxes with two vertical layers are marked by vertical stripes, and shallow one-layer boxes are white. (b) Release rates of ⁹⁰Sr in the accidental (bold line) and post-accidental (dashed line) periods.

Comparison between calculated and measured ⁹⁰Sr concentrations in water, bottom sediment and piscivorous fish for the coastal box and box no. 173 are shown in Fig. 10. Measured concentrations of ⁹⁰Sr in the water and bottom sediments before the accident were obtained from the MEXT database (MEXT, 2020). Concentrations of ⁹⁰Sr after the accident at TEPCO (Tokyo Electric Power Company) sampling points near the discharging canals (T-1 and T-2) and at different distances offshore (TD-5 inside the coastal box area and T-7, TD-1 and TD-9 for outer box no. 173) are available in the NRA (Japan Nuclear Regulation Authority) database (NRA, 2020). The results of simulation show that the concentration of ⁹⁰Sr in the seawater reaches the maximum just after the accidental release (Fig. 10a, d). Notice that the box model gives the average concentrations for each box, which means that local concentrations may be larger or smaller than the average (Fig. 10). Especially large differences can occur during the accidental release under strongly non-equilibrium conditions. Further large dispersion of measured concentrations shows that non-equilibrium conditions remained for a long time in the area close to NPP (Fig. 10a). The agreement between calculated and measured concentrations (Fig. 10a, d) could be a confirmation of the correctness of estimation of the source term. This is also confirmed by agreement of measured and simulated concentrations of ⁹⁰Sr in the bottom sediment (Fig. 10b, e). Notice that there is no measurement data for bottom sediment in box no. 173 (Fig. 10e).

Figure 10 Comparison between calculated and measured ⁹⁰Sr concentrations in water (a), bottom sediment (b) and piscivorous fish (c) for the coastal box and in water (d), bottom sediment (e) and piscivorous fish (f) for box no. 173.

The calculated concentration of ⁹⁰Sr in piscivorous fish was compared with measurement data for fat greenling (Hexagrammos otakii) before the accident (JCAC, 2020) and with data from NRA (2020) and Miki at al. (2017) after the accident. The NRA data for ⁹⁰Sr are very limited. Therefore, different species of piscivorous fish were considered, such as sharks (Triakis scyllium, Squatina japonica), rockfish (Sebastes cheni) and sea bass (Lateolabrax japonicus). Parameters of the MCKA model for these fishes are given in Table S11. The simulation results with the MCKA model agree well with the experimental observations (Fig. 10c, f), while the target tissue (TT) approach underestimates the concentration of ⁹⁰Sr in the fish.

Comparison of the Eqs. (1)–(3) and (7) of the MCKA model and Eq. (11) of the standard whole-body model demonstrates that the main difference is found in the description of the whole-body elimination rate λ_wb. Whereas in the whole-body model the value of λ_wb is constant, in the MCKA model it is the ratio of activity-weighted tissue elimination rates to whole-body activity:

$$\begin{array}{}\text{(34)}& {\mathit{\lambda }}_{\text{wb}}={\displaystyle \frac{{\sum }_{i=\mathrm{3}}^{\mathrm{5}}{\mathit{\mu }}_i{\mathit{\lambda }}_i{C}_i}{C_{\text{wb}}}}.\end{array}$$

Therefore, in the MCKA model, the value of λ_wb can vary over time, depending on the uptake of radionuclide and the tissue elimination rates. The time variation of λ_wb computed from Eq. (34) is shown in Fig. 11 for three different cases: (a) pulse-like feeding experiment (Mathews and Fisher, 2008), (b) release of ⁶⁰Co during normal operation of the Forsmark NPP and (c) accumulation of ⁹⁰Sr in the fish due to the FDNPP accident. As seen in the plots, λ_wb varies considerably when there is non-equilibrium, such as in the case of a pulse-like feeding or an accident. Even in case of a routine release, λ_wb follows any changes in the release rate (Fig. 11b).

Figure 11 The calculated λ_wb for three scenarios: (a) the pulse-like feeding experiment (Mathews, Fisher (2008), (b) the release of ⁶⁰Co during normal operation of the Forsmark NPP and (c) the accumulation of ⁹⁰Sr in the fish due to the FDNPP accident.

In case of the FDNPP accident, the calculated λ_wb shows some tendency towards an equilibrium value, but after a pulse-like release of ⁹⁰Sr in 2011, λ_wb doubled following the release of activity and then slowly converged to the quasi-equilibrium state governed by the global deposition. Notice that in this case we extended the simulation period to 2040 by extrapolating deposition data and FDNPP release data in Fig. 9b. Therefore, the whole-body model with a constant λ_wb, which is calibrated using observational data, cannot correctly describe such transient processes in the organism. This is confirmed in Fig. 10 by comparing the results from the MCKA model and the one-compartment model. Here the AE value in both models are the same, whereas the equilibrium value λ_wb was calculated in the MCKA model using the value before 2011 for piscivorous fish (λ_wb=0.0027 d⁻¹). The one-compartment model simulation results using parameters from the MCKA model are close to the MCKA model results at the initial stage of accidental release. However, over time, the concentration of ⁹⁰Sr in fish tends to equilibrium faster than the MCKA model predicts, which is explained by the time-dependent behaviour of λ_wb in the MCKA model. These results suggest also that the MCKA model could be effectively used to find equivalent parameters for the one-compartment model as was done here.

5 Conclusions

A new approach to predicting the accumulation of radionuclides in fish by taking into account heterogeneity of distribution of contamination in the organism and dependence of metabolic process rates on the fish mass was developed. The fish organism was represented by compartments for three groups of tissues or organs (muscle, bone, organs) and two input compartments representing gills and the digestive tract. The absorbed elements are redistributed between organs or tissues and then eliminated according to their metabolic function. The food and water uptake rates, elimination rate, and growth rate depend on the metabolic rate, which is scaled by the fish mass to the $\mathrm{3}/\mathrm{4}$ power, but they do not depend on the radionuclide. At the same time, the activity is distributed between the different tissues and organs according to the tissue assimilation efficiencies, which differ per radionuclide (Table 2), but do not depend on fish mass. Therefore, the transfer rates can be associated with a specific radionuclide and fish mass as shown, for example, in Tables S2–S4 and S6–S8. The position of the fish species in the trophic level also affects the concentration of activity in the organism.

This model is of intermediate complexity and provides an alternative for the basic/simplistic whole-body models and the highly advanced PBPK models. The main difference between the MCKA and whole-body models was found in the description of the whole-body elimination rate λ_wb. Whereas in the whole-body model the value of λ_wb is constant, in the MCKA model it is the ratio of activity-weighted tissue elimination rates to whole-body activity as described by Eq. (34). The elimination rate λ_wb varies considerably in non-equilibrium state of fish, such as in the case of a pulse-like feeding or an accident.

The trophic transfer factors (TTF) were calculated for five elements using assimilation efficiencies AE_f obtained from laboratory data. Among considered elements, only the caesium (Cs) level may elevate in the predator fish of the food chain (TTF > 1). This is in agreement with measurements. The concentrations of other elements (Sr, Co, Mn, Zn) decrease with the increase of trophic level (TTF <1). The kinetics of the assimilation and elimination of ¹³⁴Cs, ⁵⁷Co, ⁶⁰Co, ⁵⁴Mn and ⁶⁵Zn, which are preferably accumulated in different tissues, were analysed using the analytical solutions of a system of model equations. These solutions exhibited good agreement with the laboratory experiments for the depuration process after single feeding of fish with radiolabelled prey and with respect to uptake of activity from water. Notice that for relatively slow processes in the marine environment, transfer processes in the gills and digestive tract can be close to equilibrium, which allows for consideration of only a three-compartment (muscle, bone, organs) version of the model.

The developed MCKA model was embedded into the box model POSEIDON-R, which describes the transfer of radionuclides through the pelagic and benthic food webs. The POSEIDON-R model was applied for the simulation of the transport and fate of ⁶⁰Co and ⁵⁴Mn routinely released from Forsmark NPP located on the Baltic Sea coast of Sweden and for calculation of ⁹⁰Sr concentration in fish after the accident at Fukushima Dai-ichi NPP. Predicted concentrations of radionuclides in fish agreed well with the measurements in both case studies. It is shown that the MCKA model with the defined generic parameters could be used in different marine environments without calibration based on a posteriori information, which is important for emergency decision support systems (Periáñez et al., 2019).