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

The contamination of the Pacific Ocean by the radioactive pollutants released from the FukushimaDaiichi 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 twomarine species: the Pacific BluefinTuna (Thunnus orientalis) and the Pacific Pink Salmon (Oncorhynchus gorbuscha).


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][2][3][4][5][6][7][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 farreaching implications for decisions related to the location of nuclear power plants as well as to fishing policy.

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 biological decay of iodine-131 is 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, where  0 = (0) and  is measured in years.Similarly, the biological decay of caesium-134 is given by 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 must be solved simultaneously.The solution to (5) is 30.22 , where  0 =  (0) , where  is measured in years.
Here, we use Table 1 to present the previous information.

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 V = (/  ) ⋅   ⋅ (ln 2/ 1/2 ), where  is the mass of the radioisotope in grams,   is the atomic mass of the radioisotope,   = 6.022141793 × 10 23 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 where 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).

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: 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, V  and V  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.3896488 m 2 kg/s 2 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.15K [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 3 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 −12 m [21].Thus, the diffusion coefficient for iodine-131 in seawater is   = 1.69015 × 10 −6 m 2 /s.
We use Table 2 to represent the diffusivity coefficients.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 V  = 0.75 cos  and the velocity in the -direction is V  = 0.75 sin .
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 −15 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: 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  3 × 1 column vector by iterating first over  then over  and finally over .The following matrix equation results: 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, C is the  3 × 1 column vector of concentration values, and A and B are  3 ×  3 square matrix of coefficients.The invertibility of A, 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, where ℎ is the space step.Let {  }  3 1 be the sequence of row of the column vector C. Given the setup of the cube, the  = 0 boundary occurs where   mod  = 1,  =  occurs where   mod  = 0,  = 0 occurs where 1 ≤   mod  2 ≤ ,  =  occurs where  2 −  ≤   mod  2 ≤ 0,  = 0 occurs where 1 ≤   mod  3 ≤  2 , and  =  occurs where  3 −  2 ≤   mod  3 ≤ 0.
The matrix equation is programmed into MatLab and iterated with a time step consistent with the stability requirement [10] as follows: 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.

Ecological Processes: The Bioconcentration
Factor, the Biomagnification Factor, and Biouptake Delay Factor The total radioactivity, according to (24), is 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 −15 .
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.000390 ppm [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 for  () < 100 g  () 1200 , for  () ≥ 100 g.
The biomagnification factor (BMF) accounts for the differences in biouptake due to diet.It is given [27] by BMF () = 25, for  () < 10 g, (46) and if () ≥ 10 g, then , for planktivores 2, for benthivores 3, for omnivores 4, for piscivores. ( 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 Additionally, the concentration of both caesium-134 and caesium-137 in large piscivorous fish, denoted by   , will tend toward or, for large nonpiscivorous fish, denoted by   , toward 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().

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 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 2 /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] 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, +  (, ) +   (, )] .(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 ( 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, where  is the time required for absorption to occur.The total exposure, (, ), an organism experiences is 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 (, ) .(57)

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 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].
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 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)

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.

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 , then the standard population model would be With the inclusion of the modified death rate, (), the ordinary differential equation is whose solution is where  0 is the initial fish population.
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 (/) 3 = 450 kg.Consequently,  = 3  √ 450 ⋅ .A Pacific Bluefin Tuna gains, on average, 1 kg of mass in 4 months [45].Substitution of these values into (58) yields Expansion, combination of like terms, factorization, and logarithmation yield only one real solution  ≈ 30.3267 kg/yr.Thus,  ≈ 232.396 kg/yr.Substitution of these values into (59) under the initial condition  1 = 0 yields the proportionality constant: 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): and thus the initial population of Pacific Bluefin Tuna  0 = 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.636kg, 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 where  ≈ 30.3267 kg/yr and  ≈ 232.396 kg/yr, the biological uptake delay factor for Pacific Bluefin Tuna is for  () < 100 g, 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.772 g/day and the metabolic rate  ≈ 2.291 g/day [49], both cannot simultaneously be true, given the constraint (/) 3 = 6.8 kg.Thus, to approximate the initial population of Pacific Pink Tuna,  0 ,  and  shall be parameters.
When  = 11.772g/day or, after conversion to the appropriate units  = 4.29972 kg/year, the average weight of a Pink Salmon, given a lifespan of 3 years [49], is given by finding the average value of (58): The initial Pink Salmon population then ranges from  0 = 4250000000/5.66677= 749,986,000 to  0 = 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.
Substitution of the nutrient intake and respiration parameters into (59) yields the proportionality constant: (81) The model for the population of the Pacific Pink Salmon, after substitution of the appropriate values into (64), becomes 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 −8 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.

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.
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).
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 weightproportionality 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 shortterm 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.

Table 2 I
1.69015 × 10 −6 m 2 /s 8.92909 × 10 −7 m 2 /s 1.06586 × 10 −6 m 2 /s 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