Modelling the Effects of Radioactive Effluent on Thunnus orientalis and Oncorhynchus gorbuscha

Abstract

The contamination of the Pacific Ocean by the radioactive pollutants released from the Fukushima Daiichi Nuclear Power Plant has raised legitimate concerns over the viability of marine wildlife. We develop a modified Crank-Nicholson method to approximate a solution to the diffusion-advection-decay equation in time and three spatial dimensions to explore the extent of the effects of the radioactive effluent on two marine species: the Pacific Bluefin Tuna (Thunnus orientalis) and the Pacific Pink Salmon (Oncorhynchus gorbuscha).

1. Introduction

Following the 9.0 magnitude earthquake and consequent tsunami on March 11, 2011, in the Pacific Ocean, a nuclear crisis struck Japan. Tokyo Electric Power Company (Tepco), operators of the Fukushima Daiichi Nuclear Power Plant, struggled to contain the fallout from three melted nuclear fuel rods, nuclear fuel exposure, and compromised core reactor integrity [1–8]. To cool the plant, seawater was pumped into the reactor pool. This water was later found to have been leaking into the Pacific Ocean. Because of the brevity of the period of time during which radioactive material was discharged (a matter of days) in comparison to the length of time during which the effect is simulated (a matter of years), the effluence was assumed to be instantaneous.

The discharge of contaminated water into the sea resulted in nuclear pollution of the marine environment surrounding Fukushima Daiichi. However, because of the strength of the Kuroshio current, this pollution was spread through the entire Pacific Ocean. The average velocity of this current was used as the advection constant. The radioactive material was assumed to diffuse along a one-dimensional path consistent with the known habitats of the two populations. In response, the World Health Organization and the Food and Agriculture Organization of the United Nations put into place cautionary measures to prevent the distribution of potentially contaminated seafood [9]. Data released by Tepco and the Government of Japan regarding the nature and extent of the nuclear pollution of the Pacific waters around the Fukushima Daiichi Nuclear Power Plant allows for an analysis of the effects of the radioactive materials. It was assumed that those effects would be measured by the radioactive exposure, both physical and biological, on individuals of the two populations in mass proportion to the effect the same exposure would have had on a human being. It was also assumed that the rate of growth of the two populations was directly proportional to the difference between the respective birth rates and respective death and fishing effort rates.

A numerical solution to the diffusion-advection-decay equation using a modification to the Crank-Nicholson method will simulate the effects of the Fukushima Daiichi nuclear disaster on the Pacific Bluefin Tuna population, a species threatened by overfishing, and the Pink Salmon population, a species now in decline. This simulation has far-reaching implications for decisions related to the location of nuclear power plants as well as to fishing policy.

2. Nuclear Contaminants

Although the major nuclear contaminants dispersed from the Fukushima Daiichi facilities vary widely in their effects on the marine environment, the cumulative effects of iodine-131, caesium-134, and caesium-137 are not negligible [10].

Iodine-131 decays into a stable atom of xenon with primary beta radiation of 606 keV and primary gamma radiation of 364 keV. The physical half-life of iodine-131 is 8.04 days, whereas its biological half-life, the rate of decay while present in organic systems, ranges from 120 to 138 days [11].

The rate of decay of iodine-131 at time 𝑡 is directly proportional to the quantity of iodine-131, 𝐼(𝑡), present at time 𝑡: 𝑑𝐼(𝑡)/𝑑𝑡=−𝑘⋅𝐼(𝑡), where 𝑘>0.

With a half-life of 8.04 days, 𝑘=ln2/8.04. Thus, 𝐼(𝑡)=𝐼₀⋅2^−𝑡/8.04, where 𝑡 is measured in days, or 𝐼(𝑡)=𝐼₀⋅2^{−𝑡/.022012}, (1) where 𝑡 is measured in years.

The biological decay of iodine-131 is 𝐼_𝑏(𝑡)=𝐼_𝑏⋅2^−𝑡/𝑘, (2) where 𝐼_𝑏 is the initial amount absorbed by an organic system, 𝑡 is measured in years, and 𝑘 is the biological half-life parameter of iodine-131 such that 0.328542≤𝑘≤0.377823.

Similarly, caesium-134 decays into a stable atom of barium with primary beta radiation of 160 keV and primary gamma radiation of 1600 keV. The physical half-life of caesium-134 is 2.0652 years, whereas, is biological half-life ranges from 20 to 60 days [12, 13].

Since the rate of decay of caesium-134 is directly proportional to the quantity of caesium-134, 𝐶(𝑡), present at time 𝑡 and given a half-life of 2.0652 years, 𝐶(𝑡)=𝐶₀⋅2^{−𝑡/2.0652}, (3) where 𝐶₀=𝐶(0) and 𝑡 is measured in years. Similarly, the biological decay of caesium-134 is given by 𝐶_𝑏(𝑡)=𝐶_𝑏⋅2^−𝑡/𝑙, (4) where 𝐶_𝑏 is the initial amount absorbed by an organic system, 𝑡 is measured in years, and 𝑙 is the biological half-life parameter such that 0.054757≤𝑙≤0.164271.

Again, caesium-137 decays via beta emission into a radioactive isotope barium-137m with primary beta radiation of 190 keV. barium-137m then decays via isomeric transition into a stable atom of barium with primary beta radiation of 65 keV and primary gamma radiation of 662 keV. The physical half-life of caesium-137 is 30.22 years, whereas its biological half-life is 70 days. The physical rate of decay of barium-137m is 153 seconds with a similar biological rate of decay [13]. Only 94.4% of caesium-137 decays into barium-137m, while the remaining 5.6% of caesium-137 decays directly into a stable, nonradioactive atom of barium via beta emission [13].

Since the rate of decay of caesium-137 is directly proportional to the quantity of caesium-137, 𝑆(𝑡), present at time 𝑡, and since the rate of decay of barium-137m is directly proportional to the difference between the quantity of caesium-137 present at time 𝑡 that decays into barium-137m and the quantity of barium-137m, 𝐵(𝑡), at that time that decays into a stable, nonradioactive atom of barium-137, the system of ordinary differential equations 𝑑𝑆(𝑡)(𝑡)(𝑡)𝑑𝑡=−𝑘⋅𝑆𝑑𝐵(𝑡)(𝑡)𝑑𝑡=0.944⋅𝑘⋅𝑆−𝑗⋅𝐵 (5)(6) must be solved simultaneously.

The solution to (5) is 𝑆(𝑡)=𝑆₀⋅2^{−𝑡/30.22},where𝑆₀=𝑆(0), (7) where 𝑡 is measured in years.

Similarly, the biological decay is given by 𝑆_𝑏(𝑡)=𝑆_𝑏⋅2^{−𝑡/0.19165}. (8)

Substitution of (7) into (6) and the addition of 𝑗𝐵(𝑡) to both sides result in the linear first-order nonhomogeneous ordinary differential equation (𝑑𝐵(𝑡)/𝑑𝑡)+𝑗⋅𝐵(𝑡)=0.944⋅𝑘⋅(𝑆₀⋅2^{−𝑡/30.22}) which can be solved using the integrating factor 𝑒^{∫𝑗𝑑𝑡}=𝑒^𝑗𝑡, producing 𝐵(𝑡)=((𝑘⋅𝑆₀⋅𝑒^−𝑘𝑡)/(𝑗−𝑘))+𝐵₀. Given 𝑘=ln2/30.22, 𝑗=ln2/0.000005, and the assumption that 𝐵(0)=0, 𝐵₀=−.944⋅(ln2/30.22)⋅𝑆₀()−(),𝐵(𝑡)ln2/0.000005ln2/30.25=𝑆₀⋅𝜑⋅2^{−𝑡/30.22}−1, (9)(10) where ()𝜑=.944⋅ln2/30.22()−().ln2/0.000005ln2/30.22 (11)

Similarly, using (8), the biological decay is 𝐵_𝑏(𝑡)=𝑆_𝑏(𝑡)⋅𝜑⋅2^{−𝑡/0.19165}−1, (12) where ()𝜑=.944⋅ln2/0.19165()−().ln2/0.000005ln2/0.19165 (13)

Here, we use Table 1 to present the previous information.

Table 1 

3. Radioactivity

Radioactivity is a measure of the number of beta or gamma emissions per particle per second. Although this is a probabilistic, quantum measure, the radioactivity will be assumed constant, as per the law of large numbers. Conversion of the mass of the substance into the number of particles and multiplication by the number of decays per particle result in the formula 𝑅𝑎𝑑𝑖𝑜𝑎𝑐𝑡𝑖𝑣𝑖𝑡𝑦=(𝑚/𝑚_𝑎)⋅𝑁_𝑎⋅(ln2/𝑡_1/2), where 𝑚 is the mass of the radioisotope in grams, 𝑚_𝑎 is the atomic mass of the radioisotope, 𝑁_𝑎=6.022141793×10²³ is Avogadro’s number, 𝑡_1/2 is the half-life of the radioisotope, and radioactivity is measured in Beequerels [14].

3.1. Physical Radioactivity

Let 𝑅_𝑝(𝑡) represent the total physical radioactivity, measured in Becquerels, of the marine system at time 𝑡, where 𝑡=0 corresponds to the instance of contamination, that is, 10 April, 2011. Physical radioactivity due to iodine-131 will be denoted by 𝑅_𝐼(𝑡), physical radioactivity due to caesium-134 will be denoted by 𝑅_𝐶(𝑡), physical radioactivity due to caesium-137 will be denoted by 𝑅_𝑆(𝑡), and physical radioactivity due to barium-137m will be denoted by 𝑅_𝐵(𝑡).

The total physical radioactivity, 𝑅_𝑝(𝑡), is the summation of the functions of physical radioactivity of the constituent components, or 𝑅_𝑝(𝑡)=𝑅_𝐼(𝑡)+𝑅_𝐶(𝑡)+𝑅_𝑆(𝑡)+𝑅_𝐵(𝑡), (14) where 𝑅_𝐼(𝑡)=𝐼₀⋅2^{−𝑡/0.022012}131⋅6.022141793×10²³⋅ln2,𝑅0.022012_𝐶(𝑡)=𝐶₀⋅2^{−𝑡/2.0652}134⋅6.022141793×10²³⋅ln2,𝑅2.0652_𝑆(𝑡)=𝑆₀⋅2^{−𝑡/30.22}137⋅6.022141793×10²³⋅ln2,𝑅30.22_𝐵(𝑡)=𝑆₀𝜑⋅2^{−𝑡/30.22}−1137⋅6.022141793×10²³⋅ln2.0.000005 (15)(16)(17)(18)

3.2. Biological Radioactivity

Let 𝑅_𝑏(𝑡) represent the total biological radioactivity, measured in Becquerels, of the marine system at time 𝑡, where 𝑡=0 corresponds to the instance of absorption. Biological radioactivity due to iodine-131 will be denoted by 𝑅_𝐼,𝑏(𝑡), biological radioactivity due to caesium-134 will be denoted by 𝑅_𝐶,𝑏(𝑡), biological radioactivity due to caesium-137 will be denoted by 𝑅_𝑆,𝑏(𝑡), and biological radioactivity due to barium-137m will be denoted by 𝑅_𝐵,𝑏(𝑡).

The total biological radioactivity, that is, the radioactivity of the radionuclides after biological absorption, is the summation of the functions of biological radioactivity of the constituent components, or 𝑅_𝑏(𝑡)=𝑅_𝐼,𝑏(𝑡)+𝑅_𝐶,𝑏(𝑡)+𝑅_𝑆,𝑏(𝑡)+𝑅_𝐵,𝑏(𝑡), (19) where 𝑅_𝐼,𝑏(𝑡)=𝐼₀⋅2^−𝑡/𝑘131⋅6.022141793×10²³⋅ln2𝑘, (20) where 0.328542≤𝑘≤0.377823.

Consider 𝑅_𝐶,𝑏(𝑡)=𝐶₀⋅2^−𝑡/𝑙134⋅6.022141793×10²³⋅ln2𝑙, (21) where 0.054757≤𝑙≤0.164271.

And consider 𝑅_𝑆,𝑏(𝑡)=𝑆₀⋅2^{−𝑡/0.19165}137⋅6.022141793×10²³⋅ln2,𝑅0.19165_𝐵,𝑏(𝑡)=𝑆₀𝜑⋅2^{−𝑡/0.19165}−1137⋅6.022141793×10²³⋅ln2.0.000005 (22)(23)

3.3. Total Radioactivity

The total radioactivity is the summation of physical radioactivity and biological radioactivity, or 𝑅(𝑡)=𝑅_𝑝(𝑡)+𝑎⋅𝑅_𝑏(𝑡+𝑏), (24) where 𝑎 is the percentage of radioactive material absorbed by the organism and 𝑏 is the time required for absorption to occur (see details in Section 5).

3.4. Initial Radioactivity and Initial Mass

According to data released by Tepco [15, 16], the initial concentrations of radioactivity following the release of 11,500 metric tonnes of contaminated water into the Pacific Ocean from 09 April, 2011, to 11 April, 2011, are 310,000 Bq/L of iodine-131, 230,000 Bq/L of caesium-134, and 230,000 Bq/L of caesium-137. Given that the contaminants are trace particles, their volumes are negligible in the computation that total volume of solution released is 1.15×10⁷ L of water, under the assumption that the density of water is 1 kg/L. Multiplication of the initial concentration of radioactivity by the total volume released produces the initial radioactivity of the iodine-131: 𝑅_𝐼(0)=3.565×10¹²Bq, (25) the initial radioactivity of caesium-134: 𝑅_𝐶(0)=2.645×10¹²Bq, (26) and the initial radioactivity of caesium-137: 𝑅_𝑆(0)=2.645×10612Bq. (27)

Substitution of the initial radioactivity into functions (15), (16), and (17), respectively, yields the initial mass of iodine-131: 𝐼₀=2.46271×10⁻¹¹grams, (28) the initial mass of caesium-134: 𝐶₀=1.75354×10⁻⁹grams, (29) and the initial mass of caesium-137 𝑆₀=2.6234×10⁻⁸grams. (30)

3.5. Equations for Radioactivity

Substitution of the initial masses of iodine-131, caesium-134, and caesium-137 into (15), (16), (17), and (18), respectively, yields 𝑅_𝐼(𝑡)=2.46271×10⁻¹¹⋅2^{−𝑡/0.022012}131⋅6.022141793×10²³⋅ln2,𝑅0.022012_𝐶(𝑡)=1.75354×10⁻⁹⋅2^{−𝑡/2.0652}134⋅6.022141793×10²³⋅ln2,𝑅2.0652_𝑆(𝑡)=2.6234×10⁻⁸⋅2^{−𝑡/30.22}137⋅6.022141793×10²³⋅ln2,𝑅30.22_𝐵(𝑡)=2.6234×10⁻⁸𝜑⋅2^{−𝑡/30.22}−1137⋅6.022141793×10²³⋅ln2.0.000005 (31)(32)(33)(34)

Substitution of initial masses of iodine-131, caesium-134, and caesium-137 into (20), (21), (22), and (23), respectively, yields 𝑅_𝐼,𝑏(𝑡)=2.46271×10⁻¹¹⋅2^−𝑡/𝑘131⋅6.022141793×10²³⋅ln2𝑘, (35) where 0.328542≤𝑘≤0.377823, 𝑅_𝐶,𝑏(𝑡)=1.75354×10⁻⁹⋅2^−𝑡/𝑙134⋅6.022141793×10²³⋅ln2𝑙, (36) where 0.054757≤𝑙≤0.164271, 𝑅_𝑆,𝑏(𝑡)=2.6234×10⁻⁸⋅2^{−𝑡/.19165}137⋅6.022141793×10²³⋅ln2,𝑅0.19165_𝐵,𝑏(𝑡)=2.6234×10⁻⁸𝜑⋅2^{−𝑡/0.19165}−1137⋅6.022141793×10²³⋅ln2.0.000005 (37)(38)

4. Numerical Method for Radioactive Pollution Spread Using Diffusion Equation with Advection and Decay

The concentration of radionuclides in seawater is given by partial differential diffusion equation: 𝜕𝐶𝜕𝜕𝑡=𝐷²𝐶𝜕𝑥²+𝜕²𝐶𝜕𝑦²+𝜕²𝐶𝜕𝑧²−𝑣_𝑥𝜕𝐶𝜕𝑥−𝑣_𝑦𝜕𝐶𝜕𝑦−𝜆𝐶, (39) where 𝐶(𝑥,𝑦,𝑧,,𝑡) is the concentration of the radionuclides at location (𝑥,𝑦,𝑧) in the Pacific Ocean and time 𝑡. 𝐷 is the diffusivity coefficient of each radioactive particle in seawater, 𝑣_𝑥 and 𝑣_𝑦 are the velocities of the Kuroshio current in the 𝑥- and 𝑦-directions, respectively, and 𝜆 is the decay coefficient. We assume that there is no advection in the 𝑧-direction.

To find 𝐷 for particular matter diffusing in liquid, we use the Einstein-Stokes equation [17], 𝐷=𝐾_𝐵⋅𝑇⋅(1/6𝜋𝜂𝑟), where 𝐾_𝐵=1.3896488m² kg/s² K is the Boltzmann constant, 𝑇 is the absolute temperature in Kelvins, 𝜂 is the dynamic viscosity of the diffusion medium, and 𝑟 is the radius of the diffusing particle in meters.

As the major current involved in the diffusion of particles from Fukushima is the Kuroshio current, only its effects will be considered for simplicity. The average surface temperature of the Kuroshio current is 24^∘C or 297.15 K [18]. Thus, 𝑇=293.15K. Since the average salinity of the Kuroshio current is 34.5 ppt [18], the dynamic viscosity of seawater is approximately 1.077×10³ kg/ms [19].

Since the hydrodynamic radius is essentially the same for all isotopes of an element [20], the hydrodynamic radius of iodine-131 is 140 pm or 140×10⁻¹² m [21]. Thus, the diffusion coefficient for iodine-131 in seawater is 𝐷_𝐼=1.69015×10⁻⁶ m²/s.

The hydrodynamic radius of caesium-134 and caesium-137 is 260 pm or 260×10⁻¹² m [21]. Thus, the diffusion coefficient for caesium-134 and caesium-137 in seawater is 𝐷_𝐶=8.92909×10⁻⁷ m²/s.

The hydrodynamic radius of barium-137m is 215 pm or 215×10⁻¹²m [21]. Thus, the diffusion coefficient for barium-137m in seawater is 𝐷_𝐵=1.06586×10⁻⁶m²/s.

We use Table 2 to represent the diffusivity coefficients.

Table 2 

The average velocity of the Kuroshio current is 0.75 m/s with an angle of approximately 𝜃=𝜋/6 to the northeast [18]. Thus, the velocity in the 𝑥-direction is 𝑣_𝑥=0.75cos𝜃 and the velocity in the 𝑦-direction is 𝑣_𝑥=0.75sin𝜃.

The decay coefficient 𝜆 is the half-life constant in (1), (3), (7), and (10).

The initial value for the concentration of radionuclides from the Fukushima Daiichi Nuclear Power Plant is the initial mass of the radionuclides divided by the initial volume of water, that is, 2.43584×10⁻¹⁵ g/L.

The decay of the contaminants is governed by (39), but because of the closed nature of the system, Neumann boundary conditions are used: 𝜕𝐶𝜕𝑥_𝑥=0=0𝜕𝐶𝜕𝑦_𝑦=0=0𝜕𝐶𝜕𝑧_𝑧=0=0𝜕𝐶𝜕𝑥_𝑥=𝐿=0𝜕𝐶𝜕𝑦_𝑦=𝐿=0𝜕𝐶𝜕𝑧_𝑧=𝐿=0, (40) where 𝐿 is the dimension of the cube imposed on the Pacific Ocean, discussed further.

To solve the diffusion-advection-decay PDE, the Crank-Nicolson numerical scheme is employed. It is the average of the forward finite difference method and of the backward finite difference method with the index shifted to agree with the forward difference method. Essentially, the scheme uses the relation, as set out by the two finite difference discretizations of the diffusion-advection-decay PDE, between the values of the concentration at the six surrounding points in three-dimensional space at one time step and the values of the concentration at the same six surrounding points in the next time step, to approximate the value of the concentration at the central point over time. The resulting finite difference scheme is desirable over either of its constituents because of its unconditional stability [22]. The directions of the velocities reflect the choice of origin to be the bottom-most southeast corner of the cube, for invertibility purposes.

The 𝑛×𝑛×𝑛 cube is reshaped into an 𝑛³×1 column vector by iterating first over 𝑥 then over 𝑦 and finally over 𝑧. The following matrix equation results:𝐀⋅𝐂^𝑙+1_{𝑖,𝑗,𝑘}=𝐁⋅𝐂^𝑙_{𝑖,𝑗,𝑘}⟹𝐂^𝑙+1_{𝑖,𝑗,𝑘}=𝐀⁻¹𝐁⋅𝐂^𝑙_{𝑖,𝑗,𝑘}, (41) where 𝑙 is the time step, 𝑖 is the space step in the 𝑥-direction, 𝑗 is the space step in the 𝑦-direction, 𝑘 is the space step in the 𝑧-direction, 𝐂 is the 𝑛³×1 column vector of concentration values, and 𝐀 and 𝐁 are 𝑛³×𝑛³ square matrix of coefficients. The invertibility of 𝐀, a strictly diagonally dominant matrix, is guaranteed by the Levy-Desplanques Theorem [23] and the reflection is described in Section 2.

The boundary conditions are incorporated into the matrix using the centered difference method at each boundary, that is, at each of the six faces of the cube. For example, 𝜕𝐂𝜕𝑥_𝑥=0≈𝐂^𝑙_0,𝑗,𝑘−𝐂^𝑙_2,𝑗,𝑘2ℎ=0⟹𝐂^𝑙_0,𝑗,𝑘=𝐂^𝑙_2,𝑗,𝑘,𝜕𝐂𝜕𝑥_𝑥=𝐿≈𝐂^𝑙_{𝑛−1,𝑗,𝑘}−𝐂^𝑙_{𝑛,𝑗,𝑘}2ℎ=0⟹𝐂^𝑙_{𝑛,𝑗,𝑘}=𝐂^𝑙_{𝑛−1,𝑗,𝑘}, (42) where ℎ is the space step. Let {𝑎_𝑖}^𝑛3₁ be the sequence of row of the column vector 𝐂. Given the setup of the cube, the 𝑥=0 boundary occurs where 𝑎_𝑖mod𝑛=1, 𝑥=𝐿 occurs where 𝑎_𝑖mod𝑛=0, 𝑦=0 occurs where 1≤𝑎_𝑖mod𝑛²≤𝑛, 𝑦=𝐿 occurs where 𝑛²−𝑛≤𝑎_𝑖mod𝑛²≤0, 𝑧=0 occurs where 1≤𝑎_𝑖mod𝑛³≤𝑛², and 𝑧=𝐿 occurs where 𝑛³−𝑛²≤𝑎_𝑖mod𝑛³≤0.

The matrix equation is programmed into MatLab and iterated with a time step consistent with the stability requirement [10] as follows:𝑙𝐷⋅ℎ²≤12. (43)

By applying the previous method in MatLab, we obtained the graphs of Iodine-131 (see Figure 1). The graphs of caesium-134, caesium-137, and barium-137m can be obtained similarly.

Figure 1 

5. Ecological Processes: The Bioconcentration Factor, the Biomagnification Factor, and Biouptake Delay Factor

The total radioactivity, according to (24), is 𝑅(𝑡)=𝑅_𝑝(𝑡)+𝑎⋅𝑅_𝑏(𝑡+𝑏), (44) where 𝑎 is the percentage of radioactive material absorbed by the organism and 𝑏 is the time required for absorption to occur.

The biological radiation to which an organism is exposed occurs by the means of biouptake, which is, simply, the flux process whereby particulate matter diffuses from the source toward the organism, and the organism absorbs the particulate matter along its surface and then internalizes it about its volume [24, 25].

The parameter 𝑎 is a measure of the difference in concentration of radionuclides between the organism and the surrounding seawater.

The initial differential concentration between the organism and the surrounding seawater is assumed to be the ratio between the naturally occurring concentration of iodine-131, caesium-134, and caseium-137 in the Pacific Ocean and the initial concentration of radionuclides from the Fukushima Daiichi NPP. However, neither iodine-131, caesium-134, nor caesium-137 is naturally occurring, so the initial differential concentration is 0∶2.43584×10⁻¹⁵.

The bioconcentration factor is the equilibrium concentration between the organism and the surrounding marine ecosystem. That is, the differential concentration will tend toward the bioconcentration factor.

The bioconcentration factor for iodine-131 for fishes in any marine ecosystem is 40 for muscles, 800 for ovaries, and 110,000 for thyroid tissue [26]. The bioconcentration factor for isotopes of caesium in turbid waters for piscivorous fish is 3000/𝐾, where 𝐾 is the stable potassium concentration of water in ppm. The bioconcentration factor for isotopes of Caesium in turbid waters for nonpiscivorous fish is 1000/𝐾 where 𝐾 is the stable potassium concentration of water in ppm [26]. The stable potassium concentration of seawater is 𝐾=0.000390ppm [26]. Thus, the bioconcentration factor for isotopes of caesium for piscivorous fish is 7,692,310 and for nonpiscivorous fish is 2,564,100.

The time required to reach this equilibrium is called the biouptake delay factor. The biouptake delay factor (BUD), measured in years, accounts for the lag time experienced by larger organisms for changes in the concentration of any particulate matter in marine ecosystems. It is generally [27] given by (𝑡)=1BUD(𝑡)𝑤(𝑡)12,for𝑤<100g(𝑡)1200,for𝑤≥100g. (45)

The biomagnification factor (BMF) accounts for the differences in biouptake due to diet. It is given [27] by BMF(𝑡)=25,for𝑤(𝑡)<10g, (46) and if 𝑤(𝑡)≥10g, then BMF(𝑡)=250𝑤(𝑡),forplanktivores2,forbenthivores3,foromnivores4,forpiscivores. (47)

That is, the concentration of iodine-131, denoted by 𝐵_𝐼, in the thyroid (the most susceptible to radiation damage) of a large marine organism will tend toward 𝐵_𝐼=110,000⋅𝜙(𝑥,𝑡)as𝑡⟶[BUD(𝑡)⋅BMF(𝑡)]. (48)

Additionally, the concentration of both caesium-134 and caesium-137 in large piscivorous fish, denoted by 𝐵_𝑝, will tend toward 𝐵_𝑝=7,692,310⋅𝜙(𝑥,𝑡)as𝑡⟶[BUD(𝑡)⋅BMF(𝑡)] (49) or, for large nonpiscivorous fish, denoted by 𝐵_𝑛, toward 𝐵_𝑛=2,564,100⋅𝜙(𝑥,𝑡)as𝑡⟶[BUD(𝑡)⋅BMF(𝑡)]. (50)

Thus, the parameter 𝑎 of (44) is 𝐵_𝐼 for iodine-131 and 𝐵_𝑝 or 𝐵_𝑛 for caesium-134, caesium-137, and barium-137m. The parameter 𝑏 is the term BUD(𝑡)⋅BMF(𝑡).

6. Exposure Dose for Radioactive Material Released from Fukushima

Because radioactivity is a point-source phenomenon, its intensity is inversely proportional to the square of the distance of a target. Additionally, the measure of exposure takes into account the possible discrete energy quanta emitted with each decay.

Let 𝑋(𝑟,𝑡) represent the exposure to the marine environment at a radius of 𝑟 cm away from the Fukushima Daiichi Nuclear Power Plant and at time 𝑡, where 𝑡=0 corresponds to the instance of contamination.

The measurement for direct exposure measured in roentgen to radiation at a point 𝑟 cm away at time 𝑡 is given by the formula 𝑋(𝑟,𝑡)=5.263×10⁻⁶(𝑡)⋅∑⋅𝑅^𝑛𝑖=1𝑦_𝑖⋅𝐸_𝑖⋅𝜇_𝑒𝑛/𝜌_𝑖𝑟², (51) where 𝑅_𝑝(𝑡) is the physical radioactivity of a nuclear isotope measured in Becquerels, 𝑦_𝑖 is the fractional yield of a photon with energy 𝐸_𝑖 measured in MeV, and (𝜇_𝑒𝑛/𝜌)_𝑖 is the mass energy absorption coefficient measured in cm²/g [28].

Direct exposure due to iodine-131 will be denoted by 𝑥_𝐼(𝑟,𝑡), direct exposure due to caesium-134 will be denoted by 𝑥_𝐶(𝑟,𝑡), direct exposure due to caesium-137 will be denoted by 𝑥_𝑆(𝑟,𝑡), and direct exposure due to barium-137m will be denoted by 𝑥_𝐵(𝑟,𝑡).

Using the appropriate (𝜇_𝑒𝑛/𝜌)_𝑖 for water [29] and substitution of the appropriate information [30, 31] for the respective radionuclides produces 𝑥_𝐼(𝑟,𝑡)=5.263×10⁻⁶⋅𝑅_𝐼(𝑡)⋅0.011143𝑟²,𝑥_𝐶(𝑟,𝑡)=5.263×10⁻⁶⋅𝑅_𝐶(𝑡)⋅0.016217𝑟²,𝑥_𝑆(𝑟,𝑡)=5.263×10⁻⁶⋅𝑅_𝑆(𝑡)⋅0.001217𝑟²,𝑥_𝐵(𝑟,𝑡)=5.263×10⁻⁶⋅𝑅_𝐵(𝑡)⋅0.001217𝑟². (52)

The cumulative direct exposure, 𝑋_𝑑(𝑟,𝑡), after conversion from roentgen to sieverts [32], an individual organism experiences 𝑟 cm away is the summation of the radiation exposure from each of the radionuclides, (31) through (34), for some interval of time 0 to time 𝑡; that is, 𝑥_𝑑(𝑟,𝑡)=0.11933^𝑡₀𝑥_𝐼(𝑟,𝑡)+𝑥_𝐶(𝑟,𝑡)+𝑥_𝑆(𝑟,𝑡)+𝑥_𝐵(𝑟,𝑡)𝑑𝑡. (53)

After the organism has been directly exposed to radiation, the absorbed dose continues to damage its tissues at a rate concomitant not with the physical half-life of the radionuclides, but with the biological half-life.

Biological exposure due to iodine-131 will be denoted by 𝑥_𝐼,𝑏(𝑟,𝑡), biological exposure due to caesium-134 will be denoted by 𝑥_𝐶,𝑏(𝑟,𝑡), biological exposure due to caesium-137 will be denoted by 𝑥_𝑆,𝑏(𝑟,𝑡), and biological exposure due to barium-137m will be denoted by 𝑥_𝐵,𝑏(𝑟,𝑡).

Therefore; 𝑥_𝐼,𝑏(𝑟,𝑡)=5.263×10⁻⁶⋅𝑅_𝐼,𝑏(𝑡)⋅0.011143𝑟²,𝑥_𝐶,𝑏(𝑟,𝑡)=5.263×10⁻⁶⋅𝑅_𝐶,𝑏(𝑡)⋅0.016217𝑟²,𝑥_𝑆,𝑏(𝑟,𝑡)=5.263×10⁻⁶⋅𝑅_𝑆,𝑏(𝑡)⋅0.001217𝑟²,𝑥_𝐵,𝑏(𝑟,𝑡)=5.263×10⁻⁶⋅𝑅_𝐵,𝑏(𝑡)⋅0.001217𝑟². (54)

The cumulative biological exposure, 𝑋_𝑏(𝑟,𝑡), after conversion from roentgen to sieverts, an individual organism experiences 𝑟 cm away is the summation of the biological exposure from each of the radionuclides, (35) through (38), for some interval of time 𝑏 to time 𝑡; that is, 𝑋_𝑏(𝑟,𝑡)=0.11933^𝑡_𝑏𝑥_𝐼,𝑏(𝑟,𝑡)+𝑥_𝐶,𝑏(𝑟,𝑡)+𝑥_𝑆,𝑏(𝑟,𝑡)+𝑥_𝐵,𝑏(𝑟,𝑡)𝑑𝑡, (55) where 𝑏 is the time required for absorption to occur.

The total exposure, 𝑋(𝑟,𝑡), an organism experiences is 𝑥(𝑟,𝑡)=𝑋_𝑑(𝑟,𝑡)+𝑋_𝑏(𝑟,𝑡). (56)

To calculate the average exposure from the maximal nonlethal exposure 4 Sv to the minimal lethal exposure 6 Sv [33, 34], the mean value theorem is used. The average exposure is given by 𝑥(𝑡)=1𝑅(𝑡)/4𝜌−𝑅(𝑡)/6𝜌^{𝑟=𝑅(𝑡)/4𝜌𝑟=𝑅(𝑡)/6𝜌}𝑥(𝑟,𝑡)𝑑𝑟. (57)

7. Fish Mortality Rates Resulting from Exposure to Radioactive Material Released from Fukushima

Information of mortality rates correlated to radiation exposure dose is limited to its effects on humans [33, 35]. Unit conversion from roentgen to sievert, which measures the equivalent dose of radiation for any biological tissue, allows for the assumption that human mortality due to radiation exposure is proportional to fish mortality due to radiation exposure by a weighing factor. One can use the von Bertalanffy model for the growth of an individual fish 𝑑𝑤𝑑𝑡=𝛼𝑤^2/3−𝛽𝑤, (58) where 𝛼 is the nutrient intake constant, 𝛽 is the respiration constant, and 𝑤=𝑤(𝑡) is the weight of a fish at time 𝑡, and the assumption that human mortality rates and fish mortality rates correlated to radiation exposure are proportional by weight, to compute a proportionality constant. The von Bertalanffy growth model assumes that two factors contribute to growth: the intake of nutrients uniformly over the surface of the organism and respiration proportional to its volume [36].

To solve the growth model (𝑑𝑤/𝑑𝑡)+𝛽𝑤=𝛼𝑤^2/3, let 𝑣=𝑤/𝑤^2/3=𝑤^1/3. Then 𝑤(𝑡)=(𝛼/𝛽)³⋅[1−(𝛽/𝛼)𝑒^{−(𝛽/3)𝑡}]³, where (𝛼/𝛽)³ is the maximum mass of a fish.

Given that the average mass of an adult human is 70 kg [37] and assuming that the mortality rates of fish due to exposure to radioactivity are proportional by mass to the mortality rates of humans due to exposure to radioactivity, the proportionality constant 𝜌 for an average fish is 𝜌=70lim_𝑡→∞1/𝑡∫^𝑡₀𝑤(𝑡)=𝑑𝑡70𝛼/𝛽³. (59)

Since the excess relative risk of exposure to radiation is 15% [34, 38, 39], the mortality rate 𝑑 for a fish exposed to 𝑥(𝑡) sieverts of radiation is thus given by 𝑑(𝑡)=0.15⋅𝜌⋅𝑥(𝑡). (60)

8. Incorporation of Mortality into Fish Population Model

8.1. Overview of MatLab Program

The numerical method explained in Section 4 produces a square matrix whose entries correspond to the concentration level at a given point in the three-dimensional cube and at a given time, resulting in a four-dimensional structure. Upon each entry, an operation was performed, using (24), to determine the total radioactivity at that given point at that given time. For each spatial 3-vector, the average exposure was calculated using (57), which produced a column matrix. For each entry of the exposure values, the corresponding mortality rate was calculated using (60). This column matrix became the 𝑑(𝑡) used to determine the modified population model.

8.2. Modified Population Model

If the current fish population 𝑁(𝑡) in the Pacific Ocean, assuming no carrying capacity, is a function of its rate of growth [36] that is, the interaction between birth rate 𝐵, mortality rate 𝑀, and fishing effort 𝐹, 𝑑𝑁(𝑡)=()(𝑡),𝑑𝑡𝐵−𝑀−𝐹⋅𝑁 (61) then the standard population model would be 𝑁(𝑡)=𝑁₀⋅𝑒^{(𝐵−𝑀−𝐹)𝑡}. (62) With the inclusion of the modified death rate, 𝑑(𝑡), the ordinary differential equation is 𝑑𝑁(𝑡)=((𝑡(𝑡)𝑑𝑡𝐵−𝑀−𝐹−𝑑))⋅𝑁 (63) whose solution is 𝑁(𝑥,𝑡)=𝑁₀⋅𝑒^{(𝐵−𝑀−𝐹)𝑡−∫𝑑(𝑡)𝑑𝑡}, (64) where 𝑁₀ is the initial fish population.

9. Results

9.1. Pacific Bluefin Tuna

The most reliable sources indicate that the virgin biomass of the Pacific Bluefin Tuna population is approximately 1,500,000 metric tonnes, or 1,500,000,000 kilograms [40–43].

Sources also indicate that the typical mortality rate for the adult Pacific Bluefin Tuna is 𝑀=0.27 per year. The birth rate 𝐵 is between 0.104 and 0.195 [44].

The maximum mass of a Pacific Bluefin Tuna is (𝛼/𝛽)³=450kg. Consequently, 𝛼=√³450⋅𝛽. A Pacific Bluefin Tuna gains, on average, 1 kg of mass in 4 months [45]. Substitution of these values into (58) yields 𝛽4501−𝛼𝑒^−𝛽/93−𝛽4501−𝛼³=1. (65) Expansion, combination of like terms, factorization, and logarithmation yield only one real solution 𝛽≈30.3267kg/yr. Thus, 𝛼≈232.396kg/yr.

Substitution of these values into (59) under the initial condition 𝑡₁=0 yields the proportionality constant: 𝜌=70𝛼/𝛽³=70.450 (66)

Since the maximum age of a Pacific Bluefin Tuna is 15 years [46], the average weight of a Pacific Bluefin Tuna from year 0 to year 15 is given by finding the average value of (58): 115⋅₀¹⁵𝑤(𝑡)𝑑𝑡=449.636kg, (67) and thus the initial population of Pacific Bluefin Tuna 𝑁₀=1,500,000,000/449.636=3,336,030.

The best estimate for the commercial catch of Pacific Bluefin Tuna in 2010 is 82,543 metric tonnes, or 82,543,000 kg [22]. Since the average mass of a Pacific Bluefin Tuna is 449.636 kg, the number of fish caught in 2010 is 82,543,000/449.636=183,577. The percentage of fish caught is 183,577/3,336,030=0.055029%. The fishing effort 𝐹≈0.055029.

Since the weight function for Pacific Bluefin Tuna is given by 𝑤(𝑡)=𝛼𝛽³⋅𝛽1−𝛼𝑒^{−(𝛽/3)𝑡3}, (68) where 𝛽≈30.3267kg/yr and 𝛼≈232.396kg/yr, the biological uptake delay factor for Pacific Bluefin Tuna is (𝑡)=1BUD(𝑡)12,for𝑤<100g,𝛼/𝛽³⋅1−𝛽/𝛼𝑒^{−(𝛽/3)𝑡3}(𝑡)1200,for𝑤≥100g. (69)

Given that the Pacific Bluefin Tuna is a piscivore, its biomagnification factor is given by (𝑡)=(𝑡)BMF25,for𝑤<10g,4,forpiscivores. (70)

Therefore, the biological delay parameter 𝑏 is (𝑡)(𝑡)=BUD⋅BMF25(𝑡)112,for𝑤<10g,3,for10≤𝑤(𝑡)<100g,𝛼/𝛽³⋅1−𝛽/𝛼𝑒^{−(𝛽/3)𝑡3}(𝑡)300,for𝑤≥100g. (71)

The model for the population of the Pacific Bluefin Tuna, after substitution of the appropriate values into (64), becomes 𝑁(𝑡)=3,336,030⋅𝑒^{−(0.27+𝑑(𝑡)−𝐵+0.055029)}, (72) where 𝐵 is the birth rate parameter such that 𝑑(𝑡)=()⋅0.104≤𝐵≤0.195,0.964286𝑥(𝑡). (73)(74)

9.2. Pacific Pink Salmon

The maximum weight (𝛼/𝛽)³=6.8kg. The minimum population doubling time is between 1.4 and 4.4 years [23, 47].

The best estimate for the current biomass of Pink Salmon in the Northern Pacific Ocean is 4,250,000 metric tonnes [48], or 4,250,000,000 kg.

While the nutrient intake 𝛼≈11.772g/day and the metabolic rate 𝛽≈2.291g/day [49], both cannot simultaneously be true, given the constraint (𝛼/𝛽)³=6.8kg. Thus, to approximate the initial population of Pacific Pink Tuna, 𝑁₀, 𝛼 and 𝛽 shall be parameters.

When 𝛼=11.772g/day or, after conversion to the appropriate units 𝛼=4.29972kg/year, the average weight of a Pink Salmon, given a lifespan of 3 years [49], is given by finding the average value of (58): 1⋅3⋅365.25₀^3⋅365.25𝑤(𝑡)𝑑𝑡=5.66677kg. (75) When 𝛽=2.291g/day, or after conversion to the appropriate units 𝛽=0.836788kg/year, the average weight of a Pink Salmon is 1⋅3⋅365.25₀^3⋅365.25𝑤(𝑡)𝑑𝑡=2.08268kg. (76) The initial Pink Salmon population then ranges from 𝑁₀=4250000000/5.66677=749,986,000 to 𝑁₀=4250000000/2.08268=2,040,640,000.

The best estimate for the rate of growth of the Pacific Pink Salmon is 268,000,000 fish/year [23, 47]. The birth rate parameter is 268,000,000/2,040,640,000≤𝐵≤268,000,000/749,986,000 or as a percentage 0.131331≤𝐵≤0.35734.

The mortality rate of Pacific Pink Salmon ranges from 83.6% to 98.7% [50]. The mortality rate parameter is 0.836≤𝑀≤0.987.

The commercial catch for Pacific Pink Salmon is approximately 22,500,000 fish/year [51]. The fishing effort parameter is 22,500,000/2,040,640,000≤𝐹≤22,500,000/749,986,000 or as a percentage 0.011026≤𝐹≤0.030001.

Substitution of the nutrient intake and respiration parameters into (59) yields the proportionality constant: 𝜌=70𝛼/𝛽³=70.6.8 (77)

Since the weight function for Pacific Pink Salmon is given by 𝑤(𝑡)=𝛼/𝛽³⋅𝛽1−𝛼𝑒^{−(𝛽/3)𝑡3}, (78) where 𝛽=0.836788kg/year and 𝛼=4.29972kg/year, the biological uptake delay factor for Pacific Bluefin Tuna is (𝑡)=1BUD(𝑡)12,for𝑤<100g,𝛼/𝛽³⋅1−𝛽/𝛼𝑒^{−(𝛽/3)𝑡3}(𝑡)1200,for𝑤≥100g. (79)

Given that the Pacific Pink Salmon is a piscivore, its biomagnification factor is given by (𝑡)=(𝑡)BMF25,for𝑤<10g,3,forpiscivores. (80)

Therefore, the biological delay parameter 𝑏 is (𝑡)(𝑡)=BUD⋅BMF25(𝑡)112,for𝑤<10g,4,for10≤𝑤(𝑡)<100g,𝛼/𝛽³⋅1−𝛽/𝛼𝑒^{−(𝛽/3)𝑡3}(𝑡)400,for𝑤≥100g. (81)

The model for the population of the Pacific Pink Salmon, after substitution of the appropriate values into (64), becomes 𝑁(𝑡)=𝑁₀⋅𝑒^{−(𝑀+𝑑(𝑡)−𝐵+𝐹)(𝑡)}, (82) where 749,986,000≤𝑁₀𝑑(𝑡)≤2,040,640,000,0.836≤𝑀≤0.987,=0.15⋅70⋅6.8𝑥(𝑡),0.131331≤𝐵≤0.35734,0.011026≤𝐹≤0.030001. (83)(84)

9.3. Sensitivity Analysis

The nature of the available data necessitates an analysis of the sensitivity of the model to its parameters. The instantiation of the model with regard to the Pacific Pink Salmon fish population involves the use of more parameters than does the model of the Pacific Bluefin Tuna fish population, and so the sensitivity analysis of the former is presented. Single-factor analysis revealed that there was no change in the results of the model as a result of the parameterization of the average weight of a Pink Salmon as well as the biological half-lives. The average weight is negligible when compared to the large population of fish, and the effects due to biological half-life are spread throughout the system thin enough to be miniscule. Minuscule effects on the results were seen as a result of the parameterization of fishing rate (𝑎±0.001 change yielding a 0.22% difference), initial population (𝑎±1 change yielding a 7.2×10⁻⁸ difference), and death rate (𝑎±0.1 change yielding a 2.2% difference). Because our model changes the death rate due to effects from the radioactive effluent, sensitivity in that regard is desirable. The model was most sensitive to parameterization of birth rate (𝑎±0.1 change yielding a 2.4% difference). However, these large percent differences occurred when both the baseline solution and the test solution became sufficiently small, such that even small absolute differences produce nominally large percent differences. This is presented in Figure 2.

Figure 2 

9.4. Conclusions and Implications

A qualitative analysis of the results shows that the Pacific Bluefin Tuna will experience a steeper population decline in the short term compared to its expected population decline (see Figures 3 and 4). This decline will reach its peak a few years after the event, after which the population will return to the expected population.

Figure 3 

Figure 4 

The Pacific Pink Salmon, on the other hand, will simply decline at a faster pace than the expected population decline. That is, the radioactive effluent will result in a marked and lasting decrease in population over time (see Figures 5 and 6).

Figure 5 

Figure 6 

One explanation for the seeming incongruity in results for the Bluefin Tuna and for the Pink Salmon is the difference in average weight: Bluefin Tuna tend to be much larger than their Pink Salmon brethren, which, according the weight-proportionality assumption, will contribute significantly to the observed population effects. Although it has been experimentally shown that the effects of exposure to radiation are proportional to weight in humans, it has not been shown that those results are generalizable to other species. Experimental confirmation is needed.

Before the conclusions are subjected to social analysis, the model’s limitations must be considered. Because the distribution of fish was assumed to be uniform, the model does not reflect the migratory nature of both fish species. This does not, however, entirely negate the validity of the simulation: over a sufficiently long period of time, the short-term movement of the fish throughout the Pacific Ocean becomes negligible. Some schools of fish may be nearer to the source than others at the outset and therefore be threatened by a greater excess relative risk, while other fish may be significantly further away from the source, mitigating the effects due to radiation. The results of the model presented herein forgo the intricacies of these complicated migratory patterns and instead opt to consider the fish species’ population on average. Further research may wish to verify the validity of this simplification.

Whereas many studies into the current state of the Pacific Bluefin Tuna population are grim, the results of the simulation offer hope for the future of the species [52]. On the assumption that the ban on fishing was absolute and long standing, the rate of decline of the population, even with the effects of the radioactive effluent considered, was below the rate of decline of the population without the effects of the effluent considered. This suggests that the decline of the species may be forestalled and possibly even prevented. This, however, would require coordinated efforts on behalf of policy makers to protect the Pacific Bluefin Tuna lest the species disappear in their entirety. The results of the simulation suggest that measures beyond sustainable fishing may be required, including but not limited to periodic temporary bans on fishing.

The very high rate of decline of the Pacific Pink Salmon indicates that live specimens may contain relatively high levels of radioactivity. Continued monitoring of the Pacific Pink Salmon, as well as all marine species, for radioactivity will be critical to the avoidance of health problems for humans. Because the species migrates throughout freshwater rivers and tributaries of British Columbia, Alaska, and the Pacific Northwest of the United States, inland areas are also at risk of exposure to, at the very least, low-level radioactivity.

Moreover, the rapid rate of decline of the Pacific Pink Salmon, in conjunction with rapidly deteriorating conditions, seems to necessitate drastic action [53]. Work beyond sustainability is needed to protect the species from possible extinction [54]. Policy makers may want to consider aggressive tactics for the repopulation of the Pacific Pink Salmon.

Because of the deleterious effects on the marine environment predicted by the simulation and to prevent any future radioactive pollution of marine environments, it seems reasonable to suggest that any new nuclear power plants be constructed sufficiently far from coastal waterways so as to mitigate the absorption of any radioactive contaminants into the biosphere. This, however, would pose a risk to the environment near the nuclear power plant without the capacity of an ocean to diffuse the radioactivity. The costs and benefits of this approach should be further studied.

Acknowledgment

The authors are thankful to NSF HRD-0902132 (LSAMP) for the support to do this research.