# Numerical modelling of 137Cs content in the pelagic species of the Japanese Pacific coast following the Fukushima Dai-ichi Nuclear Power Plant accident using a size-structured food-web model

Mokrane Belharet, Sabine Charmasson, Daisuke Tsumune, Mireille Arnaud, Claude Estournel

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Mokrane Belharet, Sabine Charmasson, Daisuke Tsumune, Mireille Arnaud, Claude Estournel. Numerical modelling of 137Cs content in the pelagic species of the Japanese Pacific coast following the Fukushima Dai-ichi Nuclear Power Plant accident using a size-structured food-web model. PLoS ONE, 2019, 14 (3), pp.1-31. 10.1371/journal.pone.0212616. [https://doi.org/10.1371/journal.pone.0212616](https://doi.org/10.1371/journal.pone.0212616). This is an open access article distributed under a Creative Commons Attribution 4.0 license.

HAL Id: hal-02346985

https://hal.science/hal-02346985v1

Submitted on 3 Jun 2020

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
RESEARCH ARTICLE

# Numerical modelling of ¹³⁷Cs content in the pelagic species of the Japanese Pacific coast following the Fukushima Dai-ichi Nuclear Power Plant accident using a size-structured food-web model

<figure><img src="image_3.png" /></figure>

Mokrane Belharet<sup>1,2,*</sup>, Sabine Charmasson<sup>1</sup>, Daisuke Tsumune<sup>1,3</sup>, Mireille Arnaud<sup>1</sup>,
Claude Estournel<sup>2</sup>

1 Institut de Radioprotection et de Süreté Nucléaire (IRSN), PSE-ENV/SRTE, Laboratoire de Recherche sur les Transferts de radionucléides dans les écosystèmes Aquatiques (LRTA), Saint-Paul Lez Durance, France,

2 Laboratoire d'Aérologie (LA), UMR 5560, CNRS-Université de Toulouse, UPS, Toulouse, France,

3 Environmental Science Research Laboratory, Central Research Institute of Electric Power Industry, Abiko, Japan

* Current address: Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di Milano, Italy

* mokrane.belharet@gmail.com

<figure><img src="image_4.png" /></figure>

OPEN ACCESS

**Citation:** Belharet M, Charmasson S, Tsumune D, Arnaud M, Estournel C (2019) Numerical modelling of ¹³⁷Cs content in the pelagic species of the Japanese Pacific coast following the Fukushima Dai-ichi Nuclear Power Plant accident using a size-structured food-web model. PLoS ONE 14(3): e0212616. https://doi.org/10.1371/journal.pone.0212616

**Editor:** Rachel A. Hovel, University of Maine at Farmington, UNITED STATES

**Received:** May 2, 2018

**Accepted:** February 6, 2019

**Published:** March 13, 2019

**Copyright:** © 2019 Belharet et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files. These data were collected from the Japanese Ministry of Agriculture, Forestry and Fisheries (http://www.jfa.maff.go.jp/e/inspection/index.html), and the Japanese Nuclear Regulation Authority (http://radioactivity.nsr.go.jp/en/list/205/list-1.html). All the global ocean physics reanalysis files are available for downloading from the Copernicus marine environment monitoring service (www.

## Abstract

As result of the great east Japan earthquake on March 2011 and the damages of the Fukushima Dai-ichi Nuclear Power Plant (FDNPP), huge amount of radionuclides, especially ¹³⁷Cs, were released to the Japanese Pacific coast. By consequence, several marine species have been contaminated by direct uptake of radionuclides from seawater or through feeding on contaminated preys. In the present study we propose a novel radioecological modelling approach aiming to simulate the radionuclides transfer to pelagic marine species by giving to the organism body-size a key role in the model. We applied the model to estimate the ¹³⁷Cs content in 14 commercially important species of the North-Western Pacific Ocean after the FDNPP accident. Firstly, we validated the model and evaluated its performance using various observed field data, and we demonstrated the importance of using such modelling approach in radioecological studies. Afterwards, we estimated some radioecological metrics, such as the maximum activity concentration, its corresponding time and the ecological half-life, which are important in assessment of the previous, current and future contamination levels of the studied species. Finally, we estimated the time duration required for each species to reach the pre-accident ¹³⁷Cs activity concentrations. The results showed that the contamination levels in the planktivorous species have generally reached the pre-accident levels since about 5 years after the accident (since 2016). While in the case of the higher trophic level species, although the activity concentrations are much lower than the regulatory limit for radiocesium in seafood in Japan (100 Bq kg⁻¹), these species still require another 6–14 years (2018–2026) to reach the pre-accident levels.
<figure><img src="image_5.png" /></figure>

marine.copernicus.eu). The bathymetry data were
obtained from the Japan Oceanographic Data
Centre (JODC).

**Funding:** This work was funded by the French Institute for Radioprotection and nuclear safety (IRSN), http://www.irsn.fr/FR/Pages/Home.aspx, and supported by the French National Research Agency (ANR) in the framework of AMORAD project (ANR-11-RSNR-0002). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

**Introduction**

As result of the great East Japan earthquake on March 11, 2011 and the subsequent tsunami, loss of electric power and hence the failure of the cooling systems of the nuclear reactors of the Fukushima Dai-ichi Nuclear Power Plant (FDNPP) caused the hydrogen explosion of three of them. By consequence, a large amount of radioactive materials entered the marine environment both as fallout and direct releases into the ocean, leading to the contamination of marine species. ¹³⁷Cs and ¹³⁴Cs were released in large quantities with estimates converging on a range of 15–20 PBq for each of them [1]. However, ¹³⁷Cs has a longer half-life time (2.07 yr for ¹³⁴Cs and 30.07yr for ¹³⁷Cs) and will last longer in the environment [1,2].

The released radiocesium reached marine organisms both through uptake from the surrounding water and food ingestion [2]. The results of the monitoring conducted by the Japanese Ministry of Agriculture, Forestry and Fisheries (MAFF) just after the accident, showed that several pelagic species presented ¹³⁷Cs activity concentrations higher than 100 Bq kg⁻¹ defined as the regulatory limit for radiocesium in seafood in Japan (limit established on 1st April 2012). However, the high contamination levels characterizing the organisms immediately after the accident were followed by a decreasing trend since summer 2011 [3]. To complete this monitoring process and to understand the dynamics of marine organisms contamination, various modelling studies have been performed [4–8], aiming to both estimate the contamination levels of different organisms following this FDNPP accident and provide more scientific elements helping to understand and interpret the response of marine organisms to the high contamination of their living environments.

The mechanism of radiocesium bioaccumulation by marine organisms is quite complex due to the diversity of contamination sources. For example, a predatory fish takes up radionuclide from water via its gills, and feeds on a variety of prey items belonging to different trophic levels, each of them acquired their contaminant burden via similar mechanisms as the predator [9]. These trophic interactions lead to the contamination of the entire wildlife food web by transfer of the contaminant between the different organisms (from prey to predator). By consequence, food-web bioaccumulation models are the most useful tools to handle this complexity to better characterize and understand the relationship between radionuclide activity concentration in abiotic environmental media and that in the marine organisms.

In the marine ecosystem, size is an important parameter governing biological processes [10]. The use of size spectra to characterize the marine ecosystem is an appealing method to reduce a complex food web to a simple representation [11]. Body size exerts a critical influence on various feeding-related behaviours of individuals, such as predation and predation avoidance, as well as constraining metabolic rate and affecting the rate at which interactions occur between predators and preys [12]. Given that these processes have a substantial effect on the radionuclides transfer between species, and the fact that previous studies have already highlighted the variability of radiocesium content with respect to the organism body size or its position in the trophic chain [13–16], it is important that modern radioecological models be able to involve this fundamental individual-specific trait (the body size) as a critical governing parameter.

In this study, we propose a food-web radioecological modelling approach to simulate the transfer of artificial radionuclides, such as ¹³⁷Cs, to marine pelagic species. We applied this approach to estimate the ¹³⁷Cs activity concentrations in 14 commercially important species of the pelagic food chain of the Japanese Pacific coast, following the FDNPP accident. The purpose of this study was to show the importance of considering the organism body size as a key parameter in the radioecological models, and to better represent the inter- and intra-species variability in terms of ¹³⁷Cs activity concentration due to the size difference between the
individuals. A secondary aim was to estimate some radioecological metrics to well understand the contaminant transfer along the trophic chain and to assess the contamination levels of the studied species at different spatial and temporal scales. Finally, we estimated the time duration required for each species to reach the pre-accident activity concentrations.

## Model description and parameterization

In contaminated marine environment, the living organisms accumulate and eliminate simultaneously the radionuclides. For heterotrophs, the radionuclides can be accumulated both from the dissolved phase directly and from ingested food [17]. The transferred radionuclides are, afterwards, eliminated from the organism's body through biological (excretion and egestion) and physical (radioactive decay) processes.

The dynamics of radionuclides accumulation and elimination by a marine heterotroph organism belonging to species *s* and cohort *a* (a group of individuals from the same age-class) can be expressed by the following equation:

$$ \frac{d[R]_{s,a}}{dt} = \mu [R]_w + AE \frac{C_{s,a}}{W_{s,a}} \sum_{i=1}^{N} \sum_{j=1}^{n(i)} \alpha_{i,j} [R]_{i,j} - (\lambda_{s,a} + \lambda_R + G_{s,a}) [R]_{s,a} \quad (1) $$

where $[R]_{s,a}$ ($Bq \ g^{-1}$) represents the radionuclide activity concentration in the organism, $[R]_{i,j}$ ($Bq \ g^{-1}$) is the radionuclide activity concentration in the prey belonging to species *i* and cohort *j*, *N* is the total number of species assumed as preys, $n(i)$ represents the total number of cohorts related to the prey-species *i*. $[R]_w$ is the radionuclide activity concentration in seawater ($Bq \ l^{-1}$), while $\mu$ ($lg^{-1} \ d^{-1}$) represents the aqueous uptake rate of radionuclides by the organism directly from seawater.

$W_{s,a}$ ($g$) is the weight of the organism, and $C_{s,a}$ ($g \ d^{-1}$) being its daily ration. $AE$ (dimensionless) is the assimilation efficiency of the radionuclide by the organism, while $\alpha_{i,j} \in [0,1]$ (dimensionless) represents the proportion of prey *j* in the diet of the predator *i*.

$\lambda_{s,a}$ ($d^{-1}$) represents the biological elimination rate of radionuclide by the organism, $\lambda_R$ ($d^{-1}$) is the rate of the physical radioactive decay of the radionuclide, and $G_{s,a}$ ($d^{-1}$) being the organism growth rate.

## Accumulation from water

Uptake of ¹³⁷Cs from water mainly occurs through the gills before being distributed to the whole organism's body. The uptake rate of radiocesium directly from aqueous phase is reported to be constant and independent of the ambient Cs activity concentration [18]. Due to the lack of precise information about the influence of various environmental or biological factors such as temperature and body size on the ¹³⁷Cs uptake from seawater, we assumed, in this study, a constant uptake rate estimated from the literature, and we speculate that any possible bias generated by this parameter will not grossly impact the findings of this study, since the amount of radiocesium accumulated by the fish directly from water is negligible compared to that accumulated from diet [19].

The uptake rate estimated by [18] for the piscivorous fishes was about $14.5 \times 10^{-4} \lg^{-1} \ d^{-1}$, while [20] reported values ranging from $4 \times 10^{-4}$ to $10 \times 10^{-4} \lg^{-1} \ d^{-1}$ for the juvenile stage of the sea bream (*Sparus auratus*) exposed to ¹³⁴Cs (which has the same chemical characteristics as ¹³⁷Cs). The value used by [21] and [22] in their models for planktivorous an piscivorous fishes was $0.4 \times 10^{-4} \lg^{-1} \ d^{-1}$. Therefore, the uptake rate value used in this study was $7 \times 10^{-4} \lg^{-1} \ d^{-1}$ corresponding to the average of these reported values.
Accumulation from food

Fish must consume enough food to satisfy their energetic requirements for maintenance, growth, activity and reproduction [23]. In contaminated marine environment, the food is gen- erally more or less contaminated leading to the transfer of contaminants from preys to their predators, and the quantity of absorbed radionuclides is mainly depending on the quantity of consumed food and its composition.

The individual-based modelling strategy used in this study to quantify the absorbed radio- nuclides from ingested food is based on the estimation of: (i) the daily ration of the food con- sumed by the organism, and (ii) the diet composition (proportion of each prey species in the predator diet).

**Estimation of the daily food ration.** To estimate the daily ration of the food consumed by the organism, we used a simple energy budget equation in which energy consumed by a fish is balanced by total metabolism, waste losses and growth:

$$ \textit{Consumption} = \textit{Metabolism} + \textit{Wastes} + \textit{Growth} $$

The equation can therefore be written as:

$$ C = (R + SDA) + (U + F) + G \qquad (2) $$

Where, the metabolic energies R and SDA are the Respiration and Specific Dynamic Action rates respectively. F and U represent the waste losses due to egestion (feces) and excretion (nitrogenous wastes) processes, respectively. G represents the energy allowed to support organismal growth (growth rate). All the variables represented in this equation are expressed in grams per day (g d⁻¹).

The Respiration rate R is a function of the Standard Metabolic Rate (SMR), an activity fac- tor (R<sub>A</sub>) and a water temperature function:

$$ R = (SMR \times R_A) \times conv \times f(T^o) $$

where, conv represents the coefficient that converts the respiration rate unit from mg(O₂) h⁻¹ to g d⁻¹:

$$ conv = \frac{434J}{32 \ mg(O_2)} \times \frac{1}{CAL_f \ g \ fish^{-1}} \times \frac{24 \ h}{1 \ d} $$

and CAL<sub>f</sub> is the fish energy density.

The standard (or basal) metabolic rate SMR represents the minimum rate of energy expen- diture needed to keep a fish alive. It is represented by a simple allometric equation:

$$ SMR = a_R W^{b_R} $$

where the intercept of the allometric equation a<sub>R</sub> represents the respiration rate at 0°C of 1 g resting fish. b<sub>R</sub> is the slope of the allometric mass function for standard metabolism, and W represents the weight of the organism.

The activity factor is represented by an exponential equation, which is function of the fish swimming speed:

$$ R_A = \exp(d_r \times V) $$

with d<sub>r</sub> being a coefficient relating swimming speed to metabolism (s m⁻¹), and V (m s⁻¹) the average swimming speed, which in turn is depending on the fish size L (m):

$$ V = \varphi \times L, \text{ with } \varphi \text{ (s}^{-1}\text{)} \text{ being a constant.} $$
The temperature dependence of respiration is represented for the fish species as a simple exponential relationship:

$$f(T^\circ) = \exp(c_R \times T)$$

where T is the water temperature (°C) and $c_R = \frac{lnQ_{10}}{10}$ ($Q_{10}$ the rate at which the function increases over relatively cool water temperature).

The specific dynamic action (SDA) is the daily energetic cost of digestion, assimilation and protein turnover. We formulated SDA as a constant proportion of assimilated energy:

$$SDA = \gamma(C - F)$$

where $\gamma$ represents the specific dynamic action coefficient.

Egestion (F) is modelled as a constant proportion of consumption, while excretion (U) is formulated as a fraction of assimilated energy:

$$F = \delta \times C$$

$$U = \beta(C - F)$$

where $\delta$ and $\beta$ are constant probabilities.

With a simple algebraic rearrangement, the equation of consumption ration of an organism belonging to cohort a and species s can be written as follows:

$$C_{s,a} = \frac{R_{s,a} + G_{s,a}}{1 - \eta_s} \qquad (3)$$

Where,

$$\eta_s = (\gamma + \beta)(1 - \delta) + \delta$$

The growth rate is the change over time of the organism body mass. It is expressed as:

$$G_{s,a} = \frac{dW}{da} = \frac{W_{s,a} - W_{s,(a-da)}}{da} \qquad (4)$$

The fish body mass is related to its length, the length-weight relationship of the fish is worked out as per cubic law given by [24]:

$$W_{s,a} = q \times L_{s,a}^b \qquad (5)$$

where $L_{s,a}$ is the total fish length (cm) at age a, q (g cm⁻⁶) a constant and b an exponent usually lying between 2.7 and 3.4 [25].

To estimate the length of studied organisms, we applied the Von-Bertalanffy growth model (VBGM). This model was originally used to describe the relationship between the Age and the Length of the fish assuming a maximum length that might eventually be attained ($L_∞$). The model is expressed as:

$$L_{s,a} = L_\infty(1 - \exp(-k(a - a_0))) \qquad (6)$$

where k is a constant that defines the rate at which the growth curve approaches the asymptote, a is the fish age, and $a_0$ a constant that determines the hypothetical age at which the size of the fish is zero.

The set of ecological parameters used in this model, as well as, their collection sources are summarized in S1 and S2 Tables.
**Estimation of the organism diet composition.** Under nuclear accident conditions, radio-nuclide activity concentrations in different marine organisms are generally very heterogeneous reflecting the environmental conditions (radiocesium activity concentrations in food and water, temperature, food availability) that could be encountered in their living ecosystem or along their migration trajectories. A good knowledge of the organism diet composition is therefore crucial to properly estimate the quantity of radionuclides ingested by the predator. Usually, radioecological models used constant proportions to represent the diet composition of the organisms [21,22,26], while, in reality, the organism diet is often changing over time with respect to its ontogenetic evolution. Describing the transfer of contaminants between organisms due to feeding interactions is, therefore, challenging as food webs are complex and vary as a function of time and space.

The trophic relationships represent any alimentary relationship involving at least two individuals (predator and prey) belonging to the same ecosystem. The set of these relationships represent the trophic web, which is composed of different trophic levels, and whose base is logically occupied by the primary producers. The trophic interactions between a predator and its different preys can be directed either by a “preferential selection” where the predator selects its preys with respect to their qualities in terms of energy benefits, or by an “opportunistic selection” where the preys are selected only in function of their size and their relative abundance.

In this study, predation is supposed to be totally opportunistic and only derived by the ratio of sizes between organisms. Therefore, all organisms can be potentially predators and preys at the same time, depending on their relative size. To calculate the probability that an organism *i* of size *L*<sub>*i*</sub> can prey on another organism *j* of size *L*<sub>*j*</sub>, we used the equation reported by [27]

$$S_{i,j} = (1 + \exp\left(\beta_1\left(\rho_1 - \frac{L_i}{L_j}\right)\right))^{-1} * \left(1 - \left(1 + \exp\left(\beta_2\left(\rho_2 - \frac{L_i}{L_j}\right)\right)\right)^{-1}\right) \quad (7)$$

With $\rho_1, \rho_2, \beta_1$ and $\beta_2$, being constant positive parameters characterizing both the half saturation and the flatness of the sigmoid functions (Fig 1).

The diet composition of a predator *i* was estimated by normalizing the probabilities as follows:

$$\alpha_{i,j} = \frac{S_{i,j}}{\sum_{j=1}^{N} S_{i,j}} \quad (8)$$

The Trophic Levels (TL) of a given individual *i* could be derived from the diet composition as follows:

$$TL_i = 1 + \sum_{j}^{n} TL_j * \alpha_{i,j} \quad (9)$$

where *n* being the number of preys *j* assigned to the predator *i*.

## Assimilation efficiency of the radionuclide

Radionuclide Assimilation Efficiency (AE) is defined, here, as the assimilated fraction of the radionuclide contained in the ingested food. Reported ¹³⁷Cs Assimilation efficiency by fish varies between 0.4 and 0.96 [18,20,28,29]. For the purpose of this study, due to the lack of precise information about the variability of this parameter with respect to the fish taxa and/or its size, we assumed a constant Assimilation Efficiency rate (AE) equal to 0.75, as used in previous similar modelling studies [21,22].
<figure><img src="image_6.png" /><figcaption>Fig 1. Probability of selectivity versus prey length and predator length, calculated from Eq (7).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g001

## Biological elimination

The ¹³⁷Cs elimination rate from the organism is represented in the model by the biological and the physical elimination rates, as well as, the growth rate that plays a role in the dilution of radionuclide activity concentration in the organism body. For ¹³⁷Cs, the physical half-life is fairly high (30 years), leading to a negligible physical elimination rate compared to the biological elimination.

To estimate the ¹³⁷Cs biological elimination rate in this study, we used the allometric relationship reported in [30]. This relation estimates the elimination rate as a function of the organism weight and the water temperature (Fig 2), and can be written, after rearrangements, as follows:

$$ \lambda_{s,a} = W_{s,a}^{-0.28} * 10^{(9.03 - 0.47/k^2)} \quad (10) $$

with k being the Boltzmann constant (~8.62 * 10⁻⁵ eV K⁻¹), and T the water temperature.

## Model application

### Study area

Our study area covers the eastern coast of Japan located between 36 and 38.5°N, and between the shore and the isobath 500 m (Fig 3A). This area was subdivided into 30 evenly spaced sites (Fig 3B) to take into account the spacial variability of ¹³⁷Cs activity concentration in seawater and of water temperature. The distance between two successive is 0.25° (about 28 km). However, a large part of the data available for validation is located in the area between 36.8° – 38° N
<figure><img src="image_7.png" /><figcaption>**Fig 2.** Variability of the biological elimination rate (left) calculated using Eq (10) and the biological half-life (right), with the water temperature and the organism weight.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g002

and between 0–200 m depth. Therefore, for the purpose of model validation and for more precision, the spatial resolution in this coastal area was increased (about 0.11°) as shown in Fig 3C.

**Studied species**

The studied ecosystem consists of 14 pelagic marine species ranging from the small planktivorous fishes such as sardine and anchovy to the higher predatory fishes such as the tuna species. These species are of great importance in this region in terms of abundance and/or commercial interests. The considered species are shown in Fig 4 and listed in Table 1 by their scientific and common names.

For each species, one cohort by year was considered taking into account the spawning period. The spawning dates corresponding to the starting date for each cohort were randomly selected from a normal distribution in which the mean value represents the date at which the spawning is generally at the highest.

**Numerical simulations**

Under the assumption of a limited movement of the studied species (the species are assumed to keep their position throughout the simulation period), we carried out 11 years of simulation from January 2002 to February 2013, assuming a time step of 1 day chosen to take into account the high temporal variability of ¹³⁷Cs activity concentration in seawater after the accident.
<figure><img src="image_8.png" /><figcaption>**Fig 3.** (A) Map of Japan showing the location of the study area in the Pacific coast. (B) Positions of the 30 simulated sites (numbered from 1 to 30) used for the validation of the model in the case of the non-georeferenced data. (C) Zoom of the area with a higher spatial resolution, showing the positions of the sites used for the validation of the model in the case of the georeferenced data.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g003
<figure><img src="image_9.png" /><figcaption>Fig 4. Estimated consumption rates for the 14 studied species, confronted with the variability ranges of the consumption rates reported in the literature (blue bars). The colours of the dots correspond to the different cohorts (2002–2013). The observed data are collected from the references [39–50]. The studied species are shown in an increasing order of their asymptotic maximum Length.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g004

During the last decade before the accident (11 March 2011), the ¹³⁷Cs activity concentrations in the North Western Pacific (NWP) Ocean ranged from 1 to 3 x 10⁻³ Bq l⁻³ [31,32]. Therefore, the ¹³⁷Cs activity concentrations in seawater used in the simulations during this period were randomly selected from a uniform distribution within the interval [1, 3] x 10⁻³ Bq l⁻¹.

After the accident we used the results of the numerical simulation of ¹³⁷Cs distribution in the NWP Ocean carried out by [33], extended to cover our simulation period [34]. To simplify, and for the purpose of the modelling, we considered in this study only one group of phytoplankton with size ranging from 0.02 µm to 2 mm and only one group of zooplankton with size ranging from 0.2 mm to 3 cm. The ¹³⁷Cs activity concentrations in these two groups (Bq g⁻¹) are calculated using their activity concentration ratios derived from the simulation results reported by
<table><caption>Table 1. List of the 14 pelagic fish species considered in this study.</caption><thead><tr><td>Scientific name</td><td>Common name</td></tr></thead><tbody><tr><td>Engraulis japonicus</td><td>Japanese anchovy</td></tr><tr><td>Sardinops melanostictus</td><td>Japanese sardine</td></tr><tr><td>Clupea pallasii</td><td>Pacific herring</td></tr><tr><td>Psenopsis anomala</td><td>Rudderfish</td></tr><tr><td>Scomber japonicus</td><td>Chub mackerel</td></tr><tr><td>Cololabis saira</td><td>Pacific saury</td></tr><tr><td>Scomber australisicus</td><td>Spotted mackerel</td></tr><tr><td>Trachurus japonicus</td><td>Jack mackerel</td></tr><tr><td>Decapterus macarellus</td><td>Mackerel scad</td></tr><tr><td>Katsuwonus pelamis</td><td>Skipjack tuna</td></tr><tr><td>Scomberomorus niphonius</td><td>Japanese Spanish mackerel</td></tr><tr><td>Thunnus alalunga</td><td>Albacore</td></tr><tr><td>Thunnus albacares</td><td>Yellowfin tuna</td></tr><tr><td>Thunnus orientalis</td><td>Bluefin tuna</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0212616.t001

[35] and the ¹³⁷Cs activity concentration in seawater as follow:

$$[^{137}\text{Cs}]_{pk} = CR_{pk} \times [^{137}\text{Cs}]_{w}$$

$$[^{137}\text{Cs}]_{zk} = CR_{zk} \times [^{137}\text{Cs}]_{w}$$

$CR_{pk}$ represents the weighted mean of the activity concentration ratios obtained for the 2 groups of phytoplankton (Small phytoplankton PS and Large phytoplankton PL), and $CR_{zk}$ is the weighted mean of the activity concentration ratios estimated for the 3 groups of zooplankton (Small zooplankton ZS, Large zooplankton ZL and Predatory zooplankton ZP). For more details see [35].

For the purpose of this study, we used daily average seawater temperature data at ¼ degree horizontal resolution obtained from the global ocean physics reanalysis generated by the NEMO ocean model [36]. These data were downloaded directly from the Copernicus marine environment monitoring service (www.marine.copernicus.eu) over the period from January 2002 to February 2011. From March 2011, we used the same temperature data used by [33] extracted from JCOPE2 reanalysis data (Japan Coastal Ocean Prediction Experiment 2, [37]).

## Field data

To validate the model and to test its performance, we used two sources of field data: (1) the georeferenced data for which the geographic coordinates of the place where the fish was caught are indicated, and (2) the non-georeferenced data for which only the prefecture corresponding to the place where the fish was caught are indicated without any information about the precise geographic coordinates. For both types of data, no information on the size of sampled organisms was provided.

The non-georeferenced and a part of georeferenced data were collected from the website of the Japanese Fisheries Agency of Ministry of Agriculture, Forestry and Fisheries (http://www.jfa.maff.go.jp/e/inspection/index.html), which regularly published the results of the monitoring on radioactive materials in fishery products since March 23rd, 2011. The other part of the georeferenced data was extracted from the data reported by [38].
In this study, all the field data and the simulation results are expressed in Bq kg⁻¹ wet-weight.

## Results

### Consumption rates and trophic levels

Estimated consumption rates, using Eq (3) for the species-related cohorts are presented in Fig 4. The few available observation data collected from the literature for some species are also presented as ranges of variability.

The estimated values vary from one species to another and within the same species (between the different cohorts). These estimations range from 1 to 6% BW (Body weight) per day. They are, in general, close to the observed data reported in the literature [39–50], especially in the case of the Japanese sardine (S. melanostictus) and the Japanese saury (C. saira).

The trophic levels estimated for the studied species at the end of 2012 using the Eq (9) are shown in Fig 5. The estimated values are in general consistent with the observed data reported in the literature (Fishbase website). The small differences may be ascribed to the simplicity of our approach that assumed an ecosystem of only 14 species, while in reality the pelagic ecosystem is much more complex and may contain many other species able to be part of the organism diet. Added to that, the fact that in our approach we considered that the predation is fully opportunistic and is only determined by the predator/prey size ratio, while in the nature, some species may be selective by presenting certain preferences towards some specific preys.

In general, the results obtained in this study for both consumption rates and diet compositions were satisfactory, and we speculate that the small differences observed between the field data and simulation results would not significantly impact the findings of this study.

### Model validation

**Pre-accident period.** The field data related to the ¹³⁷Cs activity concentration in fish before the accident are very scarce and are not available for all the studied species. The validation process involving all these species is, therefore, not possible.

Values of simulated ¹³⁷Cs activity concentrations and Concentration Factor (ratio between ¹³⁷Cs activity concentration in fish and its activity concentration in seawater) for the pre-accident period are shown in Fig 6 for the 14 studied species. The mean and standard deviation represented for each species are calculated on the 1st of March 2011 involving all the species-related cohorts throughout our study area. The standard deviation corresponds to the intra-species variability (the variability between the different cohorts of each species).

The simulated ¹³⁷Cs activity concentrations for the different species range from about 0.10 ±0.05 Bq kg⁻¹ in the case of small pelagic species such as the Japanese sardine (S. melanostictus) to about 0.49±0.25 Bq kg⁻¹ for large predators such as the Yellowfin tuna (T. albacares). The high standard deviation calculated for the longest fish, reflecting the higher intra-species variability, compared to the smaller species, is due to the fact that, for the longest species, the size difference between the individuals belonging to the different cohorts is quite high, leading to a clear difference between them in terms of diet composition and biological elimination rates. In the case of the small pelagic species, however, the size difference between the cohorts is small and the diet composition is globally not different (plankton), leading to very small intra-species variability.

[51] reported a time series of ¹³⁷Cs activity concentrations in skipjack tuna (K. pelamis) between 2000 and 2010. At the end of this time series (one year before the accident), the activity concentrations were about 0.15–0.20 Bq kg⁻¹ and are, therefore, in good agreement with our finding regarding this species (0.27 ±0.12).
<figure><img src="image_10.png" /><figcaption>**Fig 5.** Comparison between the estimated trophic levels and the observed values for the 14 studied species. Black error bars represent the intra-species variability (averages and standard deviation) calculated over all the species-related cohorts. The blue bars are the reported observations in Fishbase database. The studied species are shown in an increasing order of their asymptotic maximum Length.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g005

For some species the predicted values are slightly higher than observed ones. The observed ¹³⁷Cs activity concentrations in the rudderfish *P. anomala* during the period 1984–1995 are between 0.10 and 0.26 Bq kg⁻¹ with an observed ecological half-life of about 8.6 years [14]. In 2011 (just before the accident), these observed activity concentrations should be ranging between about 0.05 and 0.12 Bq kg⁻¹ (when applying the observed ecological half-life of 8.6 years for a period of 10 years), and are therefore two times lower than the predicted values (0.21 ±0.06 Bq kg⁻¹). The predicted value of ¹³⁷Cs activity concentration in Japanese sardine *S. melanostictus* is 0.10 ± 0.02 Bq kg⁻¹, which is about two times higher than the observed activity concentrations (i.e. 0.05 ±0.02 Bq kg⁻¹) in one species of the same gender (*Sardinops sp.*) in the Japanese coasts [52].
<figure><img src="image_11.png" /><figcaption>Fig 6. Simulated ¹³⁷Cs activity concentrations (Bq kg⁻¹ wet weight) and concentration factors (l kg⁻¹) for the 14 studied species in the steady state conditions before the accident. The studied species are shown in an increasing order of their asymptotic maximum Length.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g006
The simulated Concentration Factors (CF) are ranging between 60 ± 14 l kg⁻¹ and 287 ± 150 l kg⁻¹ depending upon the species (Fig 6). The minimum value is estimated for the Japanese sardine (S. melanostictus), while the maximum value is found for the Yellowfin tuna (T. albacares). The simulated CF for the small pelagic fishes range between about 55 and 124 l kg⁻¹, and are therefore in good agreement with the CF reported by [53] for some Japanese coastal fishes over the period 1984–1990, i.e. from 14 to 133 l kg⁻¹. The average value of the simulated CF for the small pelagic fish is about 93 l kg⁻¹, which therefore well match with the International Atomic Energy Agency recommended value for the marine pelagic fish, which is 100 l kg⁻¹ [54].

**Post-accident period.** Temporal dynamics of ¹³⁷Cs activity concentrations in the studied species for which the observation field data are available are shown in Figs 7 and 8 for the period from Jan 2011 to Feb 2013. The simulation results are represented for all the species-related cohorts and are confronted with observed field data. In the Fig 7 each cohort curve illustrates the spatial mean of the highest spatial resolution area, while in Fig 8, each cohort curve corresponds to the spatial average in the prefecture to which the field data are assigned (see Fig 3). Due to the lack of information about the size of the observed fish species, the observations are confronted with the simulation results of all the species-related cohorts.

Comparisons between observed and simulated ¹³⁷Cs activity concentrations are shown in Fig 9. The results show a good agreement between observed and simulated activity concentrations. This finding means that our radioecological model successfully predicts the contamination levels observed in the marine pelagic species in the 2 years following the FDNPP accident.

## Radioecological metrics

After an accidental situation with a high radioactive contamination of the seawater, the shape of the curve representing the dynamics of the level of contamination in marine organisms generally resembles that shown in Fig 10. Although the accumulation and elimination processes occur simultaneously, the predominance of one of the two determines the shape of the curve. This curve is characterized by two successive phases separated by $T_{max}$ corresponding to the time, from the beginning of the contamination process, at which the ¹³⁷Cs activity concentration in the marine organism reaches its maximum value. These two phases are: (1) the accumulation phase that corresponds to the period over which the accumulation process predominates, leading to a rapid increase of the contaminant activity concentration in the organism; and (2) the elimination phase that represents the period over which the elimination process predominates, leading to a progressive decrease of the ¹³⁷Cs activity concentration in the organism.

**Maximal activity concentrations.** The spatial distribution of mean and standard deviation values of the maximum ¹³⁷Cs activity concentrations in the different species-related cohorts is represented in Fig 11A. The values change with respect to the location, and range between 4.7 ± 1.5 Bq kg⁻¹ (found at site 5) and 238.0 ± 36.5 Bq kg⁻¹ (in front of FDNPP at site 12). This spatial distribution is similar to that of the maximum activity concentration of ¹³⁷Cs in seawater (Fig 11C). The high activity concentrations estimated in front of FDNPP (sites 11, 12 and 18) are mainly due to the strong confinement of the contaminated waters in the vicinity of FDNPP during the high input period [33,55]. Furthermore, the maximum activity concentrations found in the northern part (Miyagi prefecture) are, in general, higher than those found in the southern part (Ibaraki prefecture), although the maximum ¹³⁷Cs activity concentrations in seawater were almost similar in the two parts (Fig 11C). This can be explained by the residence time of the contaminated waters, which is higher in the northern sites compared to the southern ones (see S1 Fig).
<figure><img src="image_12.png" /><figcaption>Fig 7. Time series of the predicted ¹³⁷Cs activity concentrations (Bq kg⁻¹ wet weight) in the different cohorts (individuals born in the same year), confronted with the observed field data (coloured symbols) of 5 studied species (for which the data are available). For each species, 11 cohorts are represented (Age1 includes individuals aged from 0 to 1 years, Age2 includes individuals aged from 1 to 2 years, etc). Each curve represents the dynamics of ¹³⁷Cs activity concentration in only one cohort. For each cohort, the shown ¹³⁷Cs activity concentration corresponds to the median of the 35 sites of the high spatial resolution area (see Fig 3C).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g007
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>
**Fig 8.** Time series of the predicted ¹³⁷Cs activity concentrations (Bq kg⁻¹ wet weight) in the different cohorts (coloured curves), confronted with the observed field data (coloured symbols). For each species, 11 cohorts are represented (Age1 includes individuals aged from 0 to 1 years, Age2 includes individuals aged from 1 to 2 years, etc). Each curve represents the dynamics of ¹³⁷Cs activity concentration in only one cohort. For each cohort, the shown ¹³⁷Cs activity concentration corresponds to the median of the sites belonging to the prefecture in which the observed individuals are caught (see Fig 3B).

https://doi.org/10.1371/journal.pone.0212616.g008

In Fig 12A and 12B the mean values of the maximum ¹³⁷Cs activity concentrations calculated for the species-related cohorts, are plotted as function of their corresponding Weights and estimated Trophic Levels respectively. These figures show moderate positive linear correlations between the maximum ¹³⁷Cs activity concentrations in the organism and the corresponding Weight and TL. This correlation is less good especially when the organism weights about 10 to 100 g (or 3 < TL < 3.7), showing a steady trend rather than an upward trend.

**Duration of the accumulation phase.** Fig 11B and 11D show, respectively, the spatial distribution of $T_{max}$ averaged over the 14 studied species, and the spatial distribution of $T_{max}$ corresponding to the maximum ¹³⁷Cs activity concentration in seawater. The species-related $T_{max}$ range between 109 and 237 days, while the seawater-related $T_{max}$ vary between 12 and 95 days depending upon the site location. The spatial distribution of species-related $T_{max}$ shows highest values in Miyagi prefecture sites (especially sites 2, 4, 5 and 6) although the seawater-related $T_{max}$ in the same sites is the lowest. This can be explained by the presence, in these sites, of a second lower peak that arrived some days after the first one (see S1 Fig) extending, thereby, the accumulation process and delaying the peak of the ¹³⁷Cs activity concentration in the organisms. In the rest of the study area, the spatial distribution of the species-related $T_{max}$ is almost the same as that of seawater-related $T_{max}$.

For each species, the cohorts-related $T_{max}$ are averaged over the whole study area and are plotted as functions of their corresponding weight (Fig 12C) and estimated Trophic Levels (Fig 12D). These figures show significant one-term exponential correlations with high coefficients of determination (R²), indicating a very good correlation between $T_{max}$ and Weight and/or TL. The growing exponential trend shows that $T_{max}$ estimated for the small organisms is significantly lower compared to the larger ones. This can be explained by the difference, in terms of diet composition, between these organisms in the sense that the smaller ones feed mainly on plankton, which requires only short time to reach equilibrium with seawater [35], leading to shorter $T_{max}$ (i.e. shorter accumulation phase). While in the case of larger organisms, the diet composition is more diverse and contains large individuals that require longer time to reach their maximum activity concentrations, leading to an extension of the accumulation period and thus to longer $T_{max}$.

**The elimination phase.** Fig 13 exhibits, for the simulated species and sites, the ratio between ¹³⁷Cs activity concentrations at the end of the simulation period (Feb 2013) and those just before the accident (March 2011). For each species, the mean value of its related cohorts is shown in the figure. The predicted ratios vary from 2.38±0.03 for the Japanese Anchovy *E. japonicus* at site 4 to about 570±142 at site 12 for the Yellowfin tuna *T. albacares*, meaning that, in the extreme situations (if the species movements are very limited throughout the period of the highest contamination of seawater), ¹³⁷Cs activity concentrations in the pelagic species, two years after the accident, might be more than 2 orders of magnitude higher than its activity concentrations before the accident, depending upon species and sites. The highest values for the studied species are mainly found at the sites located close to the damaged FDNPP (especially 12, 18, 11), in which the inter-species variability is also the highest. The results show also that, predicted ratios are generally increasing with the asymptotic maximum Length of the studied species (the studied species are shown in an increasing order of their asymptotic
<figure><img src="image_20.png" /><figcaption>Fig 9. Predicted versus observed ¹³⁷Cs activity concentrations (Bq kg⁻¹) in 6 species shown in Figs 7 and 8.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g009
<figure><img src="image_21.png" /><figcaption>Fig 10. Typical curve of the ¹³⁷Cs activity concentration dynamics in the marine organism, illustrating the accumulation and elimination phases as well as the different radioecological metrics.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g010

Lengths). This trend is most noticeable in the sites with the highest contamination levels (the closest to the FDNPP).

To measure the decreasing speed of ¹³⁷Cs concentration in the organism, we calculated the ecological half-lives that represents, here, the time (from $T_{max}$) required for the ¹³⁷Cs activity concentration in the organism to reach the half of its maximum value. The elimination phase of the contamination curve (see Fig 10) is fitted following a one-term exponential model, leading to estimation of the elimination rate. The ecological half-life ($E$) is, afterwards, calculated using the simple relationship $E = \frac{ln2}{\lambda}$. Results of the estimated ecological half-lives are shown in Fig 14 as function of the simulated sites and species. Generally, the species with the highest asymptotic lengths (Tuna species for example) show values higher than those related to the small species. In the case of the small species, the values range between 110 and 300 days with a low spatial variability (coefficient of variation: 4.4% ≤ CV ≤ 7.2%). These results are in
<figure><img src="image_22.png" /><figcaption>**Fig 11.** Spatial distribution of: (A) the maximum ¹³⁷Cs activity concentration in the organism (Bq kg⁻¹), (B) its corresponding Tₘₐₓ (in days), (C) the maximum ¹³⁷Cs activity concentration in seawater (in Bq l⁻¹), (D) the Tₘₐₓ corresponding to the maximum concentration in seawater (in days). The values shown in this Figure are calculated over the cohorts related to all the studied species.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g011

general consistent with values reported in the previous studies [56,57]. Spatial average of the simulated ecological half-life for the Japanese jack mackerel *Trachurus japonicus* is 215.7±13.0 days, this value matches perfectly the value observed off Fukushima (229±52 days) [56]. However, in the case of the large species the values vary from 360 to about 1400 days with a high spatial variability (6% ≤ CV ≤ 28%).

From the results of ¹³⁷Cs activity concentration ratios (Fig 13) and those of the ecological half-lives (Fig 14), we estimated the time, in years, required for the studied species to reach the ¹³⁷Cs activity concentrations of the pre-accident steady-state period. Fig 15 shows higher values for the large species compared to the small ones. The estimated values for the small species range between 0.33 and 5.19 years depending upon species and sites. The spatial mean values calculated for these species are ranging between 0.89±0.29 and 3.39±0.77 years. In the case of the large species, the estimated values vary from 3.45 to 16 years, and the spatial mean values lie between 6.36±0.97 and 13.53±1.1 years.

These results mean that the ¹³⁷Cs activity concentrations in the small species have already reached the pre-accident values by 2016 at the latest. While in the case of the large predatory species, even though the activity concentrations are much lower than the regulatory limit for radiocesium in seafood in Japan (100 Bq kg⁻¹), they require, on average, a duration of 6–14 years to reach the activity concentrations characterizing the pelagic organisms in the pre-accident period.
<figure><img src="image_23.png" /><figcaption>Fig 12. Predicted maximum activity concentrations (in Bq kg⁻¹) and Tmax (in days) of the species-related cohorts, plotted as a function of the their weights and trophic levels. Significance level of the correlation: p>0.05 (), p<0.05 (*), p<0.01 (**).</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g012

Discussion/Conclusion

Our model in the context of previous studies

One of the most important challenges in marine radioecology is to estimate the contamination
levels of the edible marine species in order to assess the risk related to their contamination
and, therefore, to think about the efficient strategies to protect consumers in case of uncon-
trolled releases of radionuclides in the environment. Modelling is among the important tools
being increasingly used by radioecologists to reach these objectives and to document our
knowledge on the interactions between the chemistry of these contaminants and the biology
(including physiology and ecology) of the organism.

Here we propose a novel modelling approach to simulate the dynamics of radiocesium
activity concentration in pelagic marine organisms. The model is afterwards applied to study
the contamination levels, by ¹³⁷Cs, of 14 commercially important pelagic species in the Pacific
coast of Japan following of the FDNPP accident, which is characterized by a huge
<figure><img src="image_24.png" /><figcaption>Fig 13. Ratio between the ¹³⁷Cs activity concentration in the organism at the end of the simulation (28 Feb 2013) and the activity concentration in the same organism two years before (28 Feb 2011), shown as function of studied species (x-axis) and sites (y-axis). The studied species are shown in an increasing order of their asymptotic lengths.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g013

contamination of the marine environment. This study has been strongly motivated by two main findings: (1) the lack of more sophisticated radioecological models being able to take into account some biological and environmental factors widely reported in the literature as quite important in the determination of marine organism contamination levels (e.g. the size). (2) The limited nature of the information that can help to understand the mechanism of pelagic species contamination following the FDNPP accident, due to the lack of available data and their low coverage at spatial and temporal scales.

Since it started back to the end of 1950s and the beginning of 1960s, radioecological modeling evolved and became flexible enough to incorporate environmental variability in radionuclide sources, food availability, and organism growth rates in their predictions of organism radionuclide levels [17]. In spite of these progresses, existing radioecological models are still far from proposing a realistic configuration of main processes affecting radionuclide contents in aquatic organisms. The common features of all these models lie in the following points:
<figure><img src="image_25.png" /><figcaption>Fig 14. Estimated ecological half-lives (in days) shown as function of studied species (x-axis) and sites (y-axis). The studied species are shown in an increasing order of their asymptotic lengths.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g014

* Each species in the model is considered as only one individual with constant size, while in reality, one species is generally composed of many cohorts of different sizes;

* The organism growth rate, which is explicitly taken into account in some models [21,58,59], does not generate any change in their diet. Given that ontogenetic growth is often accompanied with dietary shifts in life history parameters and creates variability in resource use and the strength of trophic interactions [60], this assumption could significantly bias the model results especially in the case of organisms with high growth rate such as the large fish species.

* The diet composition as well as the ingestion rate are considered as constant parameters that are not changing over time, while in reality, these parameters vary in parallel with the organism growth. Indeed, during their life cycle, individuals not only change in size but often also in food requirements, behaviour, trophic interaction, ...etc. [61].

It is noteworthy that many previous studies have already highlighted the existing relationship between the ¹³⁷Cs content in the organism and its body size and trophic position. For
<figure><img src="image_26.png" /><figcaption>Fig 15. Estimated time duration (in years) required for the species to reach the pre-accident activity concentrations shown as function of studied species (x-axis) and sites (y-axis). The studied species are shown in an increasing order of their asymptotic lengths.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0212616.g015

example, [62] found that ¹³⁷Cs content in brown trout from Wastewater lake after the Chernobyl accident increased with the individual weight. [13] reported that the ¹³⁷Cs activity concentrations in European hake *Merluccius merluccius* from the Mediterranean exhibit an increasing trend with the fish size in relation to change in diet composition. [63] observed the same trend for the silver carp *Hypophthalmichthys molitrix* from the cooling pond of the Chernobyl NPP, while [14] highlighted the presence of a clear positive correlation between mean weight and activity concentration of ¹³⁷Cs in 276 fish samples from the Japanese coasts. [15] and [64] highlighted the year-class-related differences in radiocesium activity concentrations in Pacific cod *Gadus macrocephalus* and Japanese flounder *Paralichthys olivaceus*, respectively, following the FDNPP accident, as result of ontogenetic and spatiotemporal changes in their diet. Other studies highlighted the good correlation between the fish size and the ¹³⁷Cs biological or ecological half-lives (so the ¹³⁷Cs retention), and they proposed, therefore, some allometric relationships connecting this parameter to the organism size [30,65,66].
Our modelling approach is thus proposed in order to complete and improve the existing radioecological models [7,21,22,67]. In this context, our model was developed based on these aforesaid existing models while giving to the organism size a critical role. Our model characteristics can be summarized in the following points:

* Individual-Based model that describes individual organisms as autonomous, unique entities.

* Age-structured model: the species cohorts are explicitly represented in the model. The main parameters affecting the content of ¹³⁷Cs in the marine organism such as the food consumption rate, the diet composition as well as the biological elimination rate are all related to the body size with respect to the findings reported in previous studies.

* Multi-species: an important number of species can be considered in the model, the interactions between the different species-related individuals are systematically determined by the model according to the predator/prey size ratio under the assumption of opportunistic predation. In our study, we considered 14 species belonging to different levels of the pelagic trophic chain of the Japanese Pacific coast, but it is quite possible, if required, to consider more species if the availability of parameters and computer resources, allow.

* Generic: although this model is until now applied only to study the transfer of one radionuclide, namely ¹³⁷Cs, it is noteworthy that this general modelling approach can be applied to other radionuclides.

In general terms, in addition to the different existing radioecological models, which have already demonstrated their efficiencies in the modelling of radionuclides transfer to the marine biota, our modelling approach is an important achievement for the radioecological modelling by including the organism size as an explicit parameter, allowing one to explicitly represent not only the variability between the different species (inter-species variability), but also among the different cohorts belonging to the same species (intra-species variability). Furthermore, this approach has demonstrated its effectiveness in both accidental and steady state situations, making it a good tool to use for forecasting and risk assessment.

## Current and future situations of the contamination levels in the pelagic species

Application of the model to simulate the ¹³⁷Cs content in 14 commercially important species in the Pacific coast of Japan highlighted the robustness of the model in both steady-state and accidental situations (before and after the accident). The estimated ¹³⁷Cs concentrations were generally in good agreement with the observations in the case of the small species, while in the case of the large species whose movements are less limited, further field observations are required to validate the results obtained in this study. It should also be noted that in spite of the spatial explicitness of our results, it is wiser that interpretations related to the ¹³⁷Cs content in a given species (especially the large species), be based on the spatial mean at a regional scale rather than on predicted values at a local scale (at a given site) to minimize the biases due to the limited nature of the fish movements assumed in this study.

The results of the simulations clearly point up the return to “normal” conditions (activity concentrations similar to those before the accident) of the majority of small pelagic species about 5 years after the accident (~ from March 2016), and that another 6–14 years would be necessary before seeing the larger species reach the pre-accident levels. This difference is mainly due to the ecological half-life between the species, which is longer in old and larger individuals. In general terms, the dynamics of ¹³⁷Cs transfer along the studied food chain following the FDNPP accident is characterized by a cascading effect due to the gradual trophic
interactions. In the first time, the contamination has mainly affected the lower levels of the trophic chain (e.g. sardine, anchovy, herring, ...) that are, in general, plankton- and larva-feeding, and which reached their maximum values, in average, 3-4 months after the accident. Afterwards, the higher trophic levels of the trophic chain (feeding on the organisms of the lower levels of the trophic chain) are affected by the contamination, reaching their maximum values 6.5-9 months after the accident. The maximum values reached by all the species were practically similar or slightly higher for the largest species (highest Trophic Levels). However, the rapidity in the contamination process characterizing the species of the lowest trophic levels (small species) compared to the species of the highest trophic levels (large species), is also accompanied by fast elimination process, leading to a more rapid return to the normal conditions in the case of the small species compared to the large ones.

It is important to note that although these results were validated with the observed field data, the forecasts related to this study are carried out under certain assumptions such as the ¹³⁷Cs activity concentrations in seawater that were obtained from dispersion simulations and the limited movement of the individuals, hence the interest to be careful when interpreting these results.

Although the results of this study are, in general, satisfactory, the uncertainty remains high, especially regarding the large species known for their long-distance migratory movements. Consequently, taking into account the fish movement in the radioecological modelling approaches, may greatly improve the predictions especially in the case of an accidental situation with a high radioactive contamination of seawater. Indeed, unlike in a steady-state situation, where the radionuclide activity concentrations are generally spatially homogeneous in both seawater and biota, in an accidental situation, activity concentrations can be very heterogeneous. Consequently, the fish species can, according to their movement trajectories, cross more or less contaminated areas, leading to larger or smaller amounts of accumulated radio-nuclides. Therefore, despite of its complexity, it is a challenge for the further modelling approaches to explicitly take these movements into account, by coupling an individual based radioecological model with a fish movement model to allow estimation of the radionuclides accumulation in parallel with the fish displacement.

**Supporting information**

S1 Fig. Two years time series of ¹³⁷Cs concentrations in the 30 sites considered in this study.

(PDF)

S1 Table. Numerical values of the parameters related to the Von Bertalanffy growth equation and the length-weight relationship used in this study for the 14 studied species.

(PDF)

S2 Table. Numerical values of the parameters used to estimate the species metabolic rates.

(PDF)

**Acknowledgments**

This work has been supported by the French National Research Agency (ANR) in the framework of AMORAD project (ANR-11-RSNR-0002). We are most grateful to Dr Yutaka Tateda of CRIEPI (Japan) for his help and support. We are also grateful for discussions and input from Hervé Thébault, Céline Duffa, Patrick Marsaleix, Thomas Duhaut, Cyril Nguyen, Caroline Ulses and Scott Fowler, that improved this paper.
Author Contributions

**Conceptualization:** Mokrane Belharet, Sabine Charmasson, Claude Estournel.

**Data curation:** Mokrane Belharet, Daisuke Tsumune, Mireille Arnaud.

**Formal analysis:** Mokrane Belharet.

**Funding acquisition:** Sabine Charmasson, Claude Estournel.

**Investigation:** Mokrane Belharet.

**Project administration:** Sabine Charmasson, Claude Estournel.

**Supervision:** Sabine Charmasson, Claude Estournel.

**Validation:** Sabine Charmasson, Daisuke Tsumune, Claude Estournel.

**Writing – original draft:** Mokrane Belharet.

**Writing – review & editing:** Mokrane Belharet.

References

1. Buesseler KO. Fukushima and Ocean Radioactivity. Oceanography. 2014; 27: 92–105.

2. Buesseler K, Dai M, Aoyama M, Benitez-Nelson C, Charmasson S, Higley K, et al. Fukushima Daiichi-Derived Radionuclides in the Ocean: Transport, Fate, and Impacts. Annu Rev Mar Sci. 2017; 9: 173–203. https://doi.org/10.1146/annurev-marine-010816-060733 PMID: 27359052

3. MEXT. Report on the monitoring of radionuclides in fishery products (March 2011–March 2014). Fisheries Agency of Japan; 2014.

4. Fiévet B, Bailly-du-Bois P, Laguionie P, Morillon M, Arnaud M, Cunin P. A dual pathways transfer model to account for changes in the radioactive caesium level in demersal and pelagic fish after the Fukushima Dai-ichi nuclear power plant accident. PLoS ONE. 2017; 12. https://doi.org/10.1371/journal.pone.0172442 PMID: 28248982

5. Tateda Y, Tsumune D, Tsubono T. Simulation of radioactive cesium transfer in the southern Fukushima coastal biota using a dynamic food chain transfer model. J Environ Radioact. 2013; 124: 1–12. https://doi.org/10.1016/j.jenvrad.2013.03.007 PMID: 23639689

6. Tateda Y, Tsumune D, Tsubono T, Misumi K, Yamada M, Kanda J, et al. Status of ¹³⁷Cs contamination in marine biota along the Pacific coast of eastern Japan derived from a dynamic biological model two years simulation following the Fukushima accident. J Environ Radioact. 2016; 151 Pt 2: 495–501. https://doi.org/10.1016/j.jenvrad.2015.05.013 PMID: 26070950

7. Vives i Batlle J. Dynamic modelling of radionuclide uptake by marine biota: application to the Fukushima nuclear power plant accident. J Environ Radioact. 2016; 151, Part 2: 502–511. https://doi.org/10.1016/j.jenvrad.2015.02.023 PMID: 25773012

8. Walters WJ, Christensen V. Ecotracer: analyzing concentration of contaminants and radioisotopes in an aquatic spatial-dynamic food web model. J Environ Radioact. 2018; 181: 118–127. https://doi.org/10.1016/j.jenvrad.2017.11.008 PMID: 29145014

9. Gobas FAPC. Food-Web Bioaccumulation Models. Encyclopedia of Ecology. 2008. pp. 1643–1652. https://doi.org/10.1016/B978-008045405-4.00398-0

10. Hildrew A, Raffaeli D, Edmonds-Brown R. Body Size: The Structure and Function of Aquatic Ecosystems. FRESHWATER BIOL. 2007. https://doi.org/10.1017/CBO9780511611223

11. Woodward G, Blanchard J, Lauridsen R, Edwards F, Jones J, Figueroa D, et al. Individual-Based Food Webs: Species Identity, Body Size and Sampling Effects. Adv Ecol Res. 2010; 43: 211–266.

12. Cohen JE, Pimm SL, Yodzis P, Saldaña J. Body Sizes of Animal Predators and Animal Prey in Food Webs. J Anim Ecol. 1993; 62: 67–78. https://doi.org/10.2307/5483

13. Harmelin-Vivien M, Bodiguel X, Charmasson S, Loizeau V, Mellon-Duval C, Tronczyński J, et al. Differential biomagnification of PCB, PBDE, Hg and Radiocesium in the food web of the European hake from the NW Mediterranean. Mar Pollut Bull. 2012; 64: 974–983. https://doi.org/10.1016/j.marpolbul.2012.02.014 PMID: 22404970

14. Kasamatsu F, Ishikawa Y. Natural variation of radionuclide ¹³⁷Cs concentration in marine organisms with special reference to the effect of food habits and trophic level. Mar Ecol Prog Ser. 1997; 160: 109–120. https://doi.org/10.3354/meps160109
15. Narimatsu Y, Sohtome T, Yamada M, Shigenobu Y, Kurita Y, Hattori T, et al. Why do the Radionuclide Concentrations of Pacific Cod depend on the Body Size? SpringerLink. 2015; 123–137. https://doi.org/10.1007/978-4-431-55537-7_10

16. Heldal HE, Føyn L, Varskog P. Bioaccumulation of ¹³⁷Cs in pelagic food webs in the Norwegian and Barents Seas. J Environ Radioact. 2003; 65: 177–185. https://doi.org/10.1016/S0265-931X(02)00095-4 PMID: 12527234

17. Fowler SW, Fisher NS. Radionuclides in the biosphere. Marine Radioactivity. Internation Atomic Energy Agency, Marine Environment Laboratory. Principality of Monaco: Hugh D. Livinston; 2004. pp. 167–203.

18. Zhao XG, Wang WX, Yu KN, Lam PKS. Biomagnification of radiocesium in a marine piscivorous fish. Mar Ecol Prog Ser. 2001; 222: 227–237. https://doi.org/10.3354/meps222227

19. Rowan DJ, Rasmussen JB. Bioaccumulation of Radiocesium by Fish: the Influence of Physicochemical Factors and Trophic Structure. Can J Fish Aquat Sci. 1994; 51: 2388–2410. https://doi.org/10.1139/f94-240

20. Mathews T, Fisher N, Jeffree R, Teyssié J. Assimilation and retention of metals in teleost and elasmo-branch fishes following dietary exposure. Mar Ecol Prog Ser. 2008; 360: 1–12. https://doi.org/10.3354/meps07462

21. Thomann R. Equilibrium-Model of Fate of Microcontaminants in Diverse Aquatic Food-Chains. Can J Fish Aquat Sci. 1981; 38: 280–296.

22. Brown J, Dowdall M, Gwynn JP, Børretzen P, Selnæs ØG, Kovacs KM, et al. Probabilistic biokinetic modelling of radiocaesium uptake in Arctic seal species: verification of modelled data with empirical observations. J Environ Radioact. 2006; 88: 289–305. https://doi.org/10.1016/j.jenvrad.2006.02.008 PMID: 16650920

23. Killen S. Energetics of Foraging Decisions and Prey Handling. Encyclopedia of fish physiology: From Genome to Environment. 2011. pp. 1588–1595.

24. Le Cren ED. The Length-Weight Relationship and Seasonal Cycle in Gonad Weight and Condition in the Perch (Perca fluviatilis). J Anim Ecol. 1951; 20: 201–219. https://doi.org/10.2307/1540

25. Froese R. Cube law, condition factor and weight–length relationships: history, meta-analysis and recommendations. J Appl Ichthyol. 2006; 22: 241–253. https://doi.org/10.1111/j.1439-0426.2006.00805.x

26. Wang WX, Ke CH, Yu KN, Lam PKS. Modeling radiocesium bioaccumulation in a marine food chain. Mar Ecol Prog Ser. 2000; 208: 41–50. https://doi.org/10.3354/meps208041

27. Maury O, Faugeras B, Shin Y-J, Poggiale J-C, Ari TB, Marsac F. Modeling environmental effects on the size-structured energy flow through marine ecosystems. Part 1: The model. Prog Oceanogr. 2007; 74: 479–499. https://doi.org/10.1016/j-pocean.2007.05.002

28. Hewett CJ, Jefferies DF. The accumulation of radioactive caesium from food by the plaice (Pleuro-nectes platessa) and the brown trout (Salmo trutta). J Fish Biol. 1978; 13: 143–153. https://doi.org/10.1111/j.1095-8649.1978.tb03422.x

29. Rowan DJ, Rasmussen JB. Measuring the bioenergetic cost of fish activity in situ using a globally dispersed radiotracer (¹³⁷Cs). Can J Fish Aquat Sci. 1996; 53: 734–745. https://doi.org/10.1139/f95-046

30. Doi H, Takahara T, Tanaka K. Trophic position and metabolic rate predict the long-term decay process of radioactive cesium in fish: a meta-analysis. PloS One. 2012; 7: e29295. https://doi.org/10.1371/journal.pone.0029295 PMID: 22279534

31. Aoyama M, Hirose K, Nemoto K, Takatsuki Y, Tsumune D. Water masses labeled with global fallout ¹³⁷Cs formed by subduction in the North Pacific. Geophys Res Lett. 2008; 35. https://doi.org/10.1029/2008gl035090

32. Hirose K, Igarashi Y, Aoyama M. Analysis of the 50-year records of the atmospheric deposition of long-lived radionuclides in Japan. Appl Radiat Isot. 2008; 66: 1675–1678. https://doi.org/10.1016/j.apradiso.2007.09.019 PMID: 18502651

33. Tsumune D, Tsubono T, Aoyama M, Uematsu M, Misumi K, Maeda Y, et al. One-year, regional-scale simulation of ¹³⁷Cs radioactivity in the ocean following the Fukushima Dai-ichi Nuclear Power Plant accident. Biogeosciences. 2013; 10: 5601–5617. https://doi.org/10.5194/bg-10-5601-2013

34. Tsumune D, Tsubono T, Aoyama M, Hirose K. Distribution of oceanic ¹³⁷Cs from the Fukushima Dai-ichi Nuclear Power Plant simulated numerically by a regional ocean model. J Environ Radioact. 2012; 111: 100–108. https://doi.org/10.1016/j.jenvrad.2011.10.007 PMID: 22071362

35. Belharet M, Estournel C, Charmasson S. Ecosystem model-based approach for modeling the dynamics of ¹³⁷Cs transfer to marine plankton populations: application to the western North Pacific Ocean after the Fukushima nuclear power plant accident. Biogeosciences. 2016; 13: 499–516. https://doi.org/10.5194/bg-13-499-2016
36. Madec G. NEMO ocean engine. Paris, France: Institut Pierre-Simon Laplace (IPSL); 2008. Report No.: 27.

37. Miyazawa Y, Zhang R, Guo X, Tamura H, Ambe D, Lee J-S, et al. Water mass variability in the western North Pacific detected in a 15-year eddy resolving ocean reanalysis. J Oceanogr. 2009; 65: 737. https://doi.org/10.1007/s10872-009-0063-3

38. Wada T, Fujita T, Nemoto Y, Shimamura S, Mizuno T, Sohtome T, et al. Effects of the nuclear disaster on marine products in Fukushima: An update after five years. J Environ Radioact. 2016; 164: 312–324. https://doi.org/10.1016/j.jenvrad.2016.06.028 PMID: 27552655

39. Sharp G D., Francis R C. An energetics model for the exploited yellowfin tuna, Thunnus albacares, population in the eastern Pacific Ocean. Fish Bull. 1976; 74: 36–51.

40. Ferriss BE, Essington TE. Can fish consumption rate estimates be improved by linking bioenergetics and mercury mass balance models? Application to tunas. Ecol Model. 2014; 272: 232–241. https://doi.org/10.1016/j.ecolmodel.2013.10.010

41. Kurita Y, Sugisaki H. Changes in the daily ration size of Pacific saury, *Cololabis saira*, with seasons and body sizes. Bull Jpn Soc Fish Oceanogr Jpn. 2004;

42. Ménard F, Stéquert B, Rubin A, Herrera M, Marchal É. Food consumption of tuna in the Equatorial Atlantic ocean: FAD-associated versus unassociated schools. Aquat Living Resour. 2000; 13: 233–240. https://doi.org/10.1016/S0990-7440(00)01066-4

43. Nikolioudakis N, Palomera I, Machias A, Somarakis S. Diel feeding intensity and daily ration of the sardine Sardina pilchardus. Mar Ecol Prog Ser. 2011; 437: 215–228. https://doi.org/10.3354/meps09275

44. Noguchi M, Takahashi Y, Hirota Y. The relationship between food consumption and growth of sardine Sardinops melanostictus by the feeding experiment. Bull Jpn Sea Natl Fish Res Inst Jpn. 1990;

45. O'Connell CP, Zweifel JR. A laboratory study of particulate and filter feeding of the Pacific mackerel, Scomber japonicus Fishery. Fish Bull. 1972; 70: 973–981.

46. Okunishi T, Yamanaka Y, Ito S. A simulation model for Japanese sardine (*Sardinops melanostictus*) migrations in the western North Pacific. Ecol Model. 2009; 220: 462–479. https://doi.org/10.1016/j.ecolmodel.2008.10.020

47. Olson RJ, Boggs CH. Apex Predation by Yellowfin Tuna (*Thunnus albacares*): Independent Estimates from Gastric Evacuation and Stomach Contents, Bioenergetics, and Cesium Concentrations. Can J Fish Aquat Sci. 1986; 43: 1760–1775. https://doi.org/10.1139/f86-220

48. Pillar SC, Barange M. Feeding habits, daily ration and vertical migration of the cape horse mackerel off South Africa. South Afr J Mar Sci. 1998; 19: 263–274. https://doi.org/10.2989/025776198784126683

49. Sakurai Y, Ikeda Y, Shimizu M, Shimazaki K. Feeding and growth of captive adult Japanese common squid, Todarodes pacificus, measuring initial body size by cold Anesthesia. Recent advances in cephalopod fisheries biology. Tokai University press. Tokyo; 1993. pp. 467–476.

50. Young JW, Lansdell MJ, Campbell RA, Cooper SP, Juanes F, Guest MA. Feeding ecology and niche segregation in oceanic top predators off eastern Australia. Mar Biol. 2010; 157: 2347–2368. https://doi.org/10.1007/s00227-010-1500-y

51. Povinec PP, Hirose K, Aoyama M. Pre-Fukushima Radioactivity of the Environment. In: Aoyama PPPH, editor. Fukushima Accident. Boston: Elsevier; 2013. pp. 277–323.

52. JCAC. Radioactivity survey data in Japan: Environmental and dietary materials. Chiba, Japan: Japan Chemical Analysis Center; 2004. Report No.: 139.

53. Tateda Y, Koyanagi T. Concentration factors for 137Cs in Japanese coastal fish (1984–1990). J Radiat Res (Tokyo). 1996; 37: 71–79.

54. IAEA. Sediment distribution coefficients and concentration factors for biota in the marine environment. Vienna: International Atomic Energy Agency; 2004. Report No.: 422.

55. Estournel C, Bosc E, Bocquet M, Ulses C, Marsaleix P, Winiarek V, et al. Assessment of the amount of cesium-137 released into the Pacific Ocean after the Fukushima accident and analysis of its dispersion in Japanese coastal waters. J Geophys Res-Oceans. 2012; 117: C11014. https://doi.org/10.1029/2012JC007933

56. Iwata K, Tagami K, Uchida S. Ecological half-lives of radi cesium in 16 species in marine biota after the TEPCO's Fukushima Daiichi nuclear power plant accident. Environ Sci Technol. 2013; 47: 7696–7703. https://doi.org/10.1021/es400491b PMID: 23768046

57. Wada T, Nemoto Y, Shimamura S, Fujita T, Mizuno T, Sohtome T, et al. Effects of the nuclear disaster on marine products in Fukushima. J Environ Radioact. 2013; 124: 246–254. https://doi.org/10.1016/j.jenvrad.2013.05.008 PMID: 23831549

58. Mathews T, Fisher NS. Trophic transfer of seven trace metals in a four-step marine food chain. Mar Ecol Prog Ser. 2008; 367: 23–33. https://doi.org/10.3354/meps07536
59. Reinfelder JR, Fisher NS, Luoma SN, Nichols JW, Wang WX. Trace element trophic transfer in aquatic organisms: a critique of the kinetic model approach. Sci Total Environ. 1998; 219: 117–135. PMID: 9802246

60. Nakazawa T, Ushio M, Kondoh M. Scale dependence of predator–prey mass ratio: determinants and applications. Adv Ecol Res Role Body Size Multispecies Syst. 2011; 45: 269.

61. Grimm V, Railsback S. Individual-based Modeling and Ecology: (Princeton Series in Theoretical and Computational Biology). Princeton University Press; 2005.

62. Elliott JM, Hilton J, Rigg E, Tullett PA, Swift DJ, Leonard DRP. Sources of Variation in Post-Chernobyl Radiocaesium in Fish from Two Cumbrian Lakes (North-West England). J Appl Ecol. 1992; 29: 108–119. https://doi.org/10.2307/2404354

63. Koulikov AO, Ryabov IN. Specific cesium activity in freshwater fish and the size effect. Sci Total Environ. 1992; 112: 125–142. PMID: 1574700

64. Kurita Y, Shigenobu Y, Sakuma T, Ito S. Radiocesium Contamination Histories of Japanese Flounder (*Paralichthys olivaceus*) After the 2011 Fukushima Nuclear Power Plant Accident. SpringerLink. 2015; 139–151. https://doi.org/10.1007/978-4-431-55537-7_11

65. Rowan DJ, Rasmussen JB. The Elimination of Radiocaesium from Fish. J Appl Ecol. 1995; 32: 739–744. https://doi.org/10.2307/2404813

66. Vives i Batlle J, Wilson RC, McDonald P. Allometric methodology for the calculation of biokinetic parameters for marine biota. Sci Total Environ. 2007; 388: 256–69. https://doi.org/10.1016/j.scitotenv.2007.07.048 PMID: 17878058

67. Erichsen AC, Konovalenko L, Møhlenberg F, Closter RM, Bradshaw C, Aquilonius K, et al. Radionuclide transport and uptake in coastal aquatic ecosystems: a comparison of a 3D dynamic model and a compartment model. Ambio. 2013; 42: 464–475. https://doi.org/10.1007/s13280-013-0398-2 PMID: 23619804