Long-term prediction of [Formula: see text]Cs in Lake Onuma on Mt. Akagi after the Fukushima accident using fractional diffusion model.

The Fukushima Daiichi Nuclear Power Plant accident also contaminates lakes in Japan. Especially in closed lakes, there is a problem of prolonged low-level [Formula: see text]Cs contamination because the activity concentration of [Formula: see text]Cs declines sharply immediately after the accident, but then begins to decrease slowly. In this paper, we derived a long-term prediction formula based on the fractional diffusion model (FDM) for the temporal variation in [Formula: see text]Cs activity concentrations of the water in Lake Onuma on Mt. Akagi, one of the closed lakes, and of pond smelt (Hypomesus nipponensis), a typical fish species inhabiting in the lake. The formula reproduced well the measured [Formula: see text]Cs activity concentration of the lake water and pond smelt for 5.4 years after the accident. Next, we performed long-term prediction for 10,000 days using this formula and compared it with the prediction results of the two-component decay function model (TDM), which is the most common model. The results suggest that the FDM prediction will lead to a longer period of contamination with low-level [Formula: see text]Cs than the TDM prediction.

Introduction

During the Fukushima Daiichi Nuclear Power Plant (FDNPP) accident, aquatic ecosystems are contaminated with radioactive Cs[1-6]. At present, the activity concentration of Cs in river water is quite low[7-9]. However, in closed lakes, there is a problem of long-tailed temporal change in low-level radioactive contamination because the activity concentration declines steeply immediately after the accident, but then turns to a gradual decrease[6,10]. Such a problem is also observed in closed lakes and small lakes after the Chernobyl accident[11,12]. Lake Onuma on Mt. Akagi in Gunma Prefecture located about 190 km from the FDNPP is a semi-closed system with limited amount of inflow and runoff water and an average water residence time of 2.3 years[10]. In August 2011, Cs plus Cs activity concentration of pond smelt (Hypomesus nipponensis), a typical fish species inhabiting in Lake Onuma exceeded 500 Bq  the standard limits of radiocesium of foods in Japan at that time. As a result, the future of radioactive contamination in closed lakes became a matter of public concern. Therefore, the Cs activity concentrations of the lake water and pond smelt in Lake Onuma were measured for a long period of 5.4 years after the accident[10].

The main processes governing Cs activity concentration after the accident can be categorized into an initial phase and mid-, to long-term phases, corresponding to the post-accident phase. The most important processes in the initial phase are the atmospheric fallout of Cs on the water surface, subsequent sorption and fixation by sediments, and their settling on the bottom[13-15]. Besides pure fitting, a parameter for a given lake and radionuclide in the diffusion model can be estimated using characteristics with clear meaning Cs deposition, effective diffusion coefficient, and distribution coefficient[11], and thus predictions of Cs long-term dynamics in a lake can be provided immediately after the accident using previously measured values of the characteristics, literature data or expert judgments. Studies in the mid-, to long-term phases after the accident include models using sum of exponential functions[13,16,17] and process-level diffusion models[11,18,19]. The two-component decay function model (TDM)[16,17], which is the sum of two exponential functions, was used to fit the measured Cs activity concentration and showed good agreement with the measured values[6,10,12]. On the other hand, by using the diffusion models, Bulgakov et al.[11] and Konoplev et al.[18,19] succeeded in showing the temporal trend of radiocesium in the mid- and long-term phases after the accidents at Chernobyl and Fukushima.

Various factors are considered to be responsible for the temporal variation in Cs activity concentration of lake water, including sedimentation on the lake sediment due to adsorption on inorganic particles, elution from the lake sediment[20,21], circulation in the lake due to convection[22-25], turbulent mixing[26], absorption by plankton and other organisms, and many other factors are possible. A model that describes such complex diffusion is the fractional diffusion model (FDM)[27,28], which includes a fractional order differentiation[29] with respect to time. The FDM has been used in the study of actual mass transfer to explain phenomena that are different from classical diffusion, especially those that diffuse slower than the classical diffusion. Some of the earliest engineering applications of the FDM include studies of concentration of an electroactive species at the surface of an electrode[30] and transient photocurrent in inorganic and organic amorphous materials[31]. Other examples include applications to fluids in porous media[32] and sorption and convection of ions in heterogeneous media[33,34] and studies on turbulence[35-37]. Based on these backgrounds, we will try to reproduce the temporal variation in Cs activity concentration of the water in Lake Onuma using the FDM.

The actual measurement period of the water in Lake Onuma is the one of the longest for a closed lake in Japan, 5.4 years after the FDNPP accident. However, there are less measured values of Cs after that period. Since there are less measured values of Cs activity concentration of the water and pond smelt in Lake Onuma since 5.4 years after the accident and not publicly available, the long-term temporal changes in Cs activity concentration of the lake water and pond smelt over 10,000 days since the accident have not been clarified.

The main goals of the present study are to (1) derive a long-term prediction formula for the temporal variation in Cs activity concentration of the water and pond smelt in Lake Onuma based on the FDM, (2) conduct long-term prediction of Cs activity concentration of the water and pond smelt in Lake Onuma for 10,000 days using this prediction formula, (3) clarify the difference between the predictions of the FDM and of the TDM. Such a prediction formula will enable long-term prediction of radioactive contamination in closed lakes and provide long-term prospects to the residents living in the vicinity, thereby relieving their psychological anxiety and preventing the spread of rumors.

Mathematical modeling

Roche et al.[38] used the continuous time random walk (CTRW)[39] to analyze the mass transport within the sediments of Lake DePue, a backwater lake of the Illinois River. The results showed that the CTRW could reproduce the concentration in the sediments, especially from the sediment interior at the lake bottom to the water interface. On the other hand, the present study examined the applicability of the CTRW up to the upper water above the sediment. Based on the study of Metzler et al.[40], we generalized diffusion of Cs from the CTRW to the fractional diffusion model (FDM), and derived the fractional diffusion equation (FDE). The derivation of the FDE is described below.

In the CTRW, the jump length of a particle (Cs) and the waiting time between jumps are represented by a probability density function (PDF) . We consider the jump length PDF and the waiting time PDF w(t) are independent variables, and one finds the separation of variable .We describe the CTRW model using by the following master equation[39]:where is the PDF when the particle arrives at time t, position x, and the second term on the right-hand side corresponds to the initial condition of the random walk.

We Fourier transform Eq. (1) from position x to k and obtain . We Laplace transform from time t to s and obtain . In this study, the Fourier transform and the inverse Fourier transform are defined, respectively, asand

The Fourier-Laplace transform of Eq. (1) leads to the solution

Now applying the method of separation of variables and Fourier-Laplace representation, the PDF in Eq. (1) can be written as follows:where and are the jump length PDF in the Fourier space and waiting time PDF in the Laplace space, respectively.

Let be the survival probability that the waiting time on a site exceed t:

The Laplace transform of corresponds to

Equation (1) leads to the PDF of C(x, t) for Cs at x at time t is given by

Combining Eqs. (4), (5) and (7), we get the Fourier-Laplace transform of the PDF C(x, t):

Let us consider a situation in which behaves as a Gaussian jump length PDF:where is given by

We Fourier transform and expanded to the order of , Eq. (10) can be approximated by the following formula:

Suzuki et al.[10] reported the decay processes of Cs activity concentration of the lake water in Lake Onuma were well suited by not the single-component decay function model (SDM) but the two-component decay function model (TDM) between 234–1981 days from March, 15, 2011. This implies that there is a long-tailed (heavy-tailed) distribution with respect to time t. Accordingly, we assume that the PDF of waiting time w(t) behaves asymptotically as a power law, and w(t) is described as follows[31,41,42]:where

The Laplace transform of Eq. (13) is given by the following expression for a small s (corresponding to a large time t)[41]:

Substituting Eqs. (12) and (15) into the Eq. (9) gives:

We can neglect the mixed term of the higher order[41] and obtainwhere

In order to derive FDE from Eq. (17), we employ the fractional derivative of an arbitrary order defined by so-called Riemann–Liouville fractional integral as follows:where n is positive integer and . Employing the integral rule for the Laplace transform of fractional derivative[29]we inverse Laplace transform from s to t in Eq. (17):and inverse Fourier transform from k to x in Eq. (21) leads to the following expression for C(x, t):

Let us apply the differential operator to Eq. (22), we finally obtain the FDE:

Our next step is to find the solution of Eq. (23). We look for the solution in the form:with expansion coefficients dep ending on the initial conditions and use the separation of variables ansatz:

Substituting Eq. (25) into Eq. (23), one obtainswhere Eq. (26) can be rewritten as follows:

From Eq. (27), we get the following decoupled set for the equation for the spatial and the temporal eigenfunctions using the constant value :and

The general solution of Eq. (28) is given as follows:where prefactors and should be determined by the boundary conditions.

The solution to Eq. (29) is given through the Mittag-Leffler function [43,44] :and written as follows:where

The full solution of the FDE (Eq. (23)) is given by the sum over all eigenfunctions, i.e., by

Since our purpose is to make long-term predictions, the fundamental mode in Eq. (34) becomes dominant for a sufficiently large time t. Therefore, C(x, t) can be approximately written as follows:where . Now, let us consider averaging Eq. (35) over space (distance x) to make time-dependent function, therefore we get an equation based on the FDM as follows:where is the radioactive decay constant per day of Cs and

When we calculated using Eq. (36), it becomes numerically unstable as the value of the series expansion term k in Eq. (31) increased. In order to avoid this phenomenon, the formula proposed by Gorenflo et al.[45] may be used for the calculation of the Mittag-Leffler function. That is, if and , then can be rewritten as:

To calculate the Eq. (36), replace t in Eq. (38) with . Equation (36) is thus

We calculated the integral in Eq. (39) by using SciPy[46] which is an open-source scientific computing library for the Python programming language.

The physical meaning of the four parameters is explained below. The parameter A is the initial average activity concentration of Cs (Bq ) in the whole lake. The parameter A is averaged with respect to the water depth, as defined in Eq. (37). The parameter represents the speed of diffusion: a large value of represents slow diffusion, while a small value of represents fast diffusion. Specifically, it is expressed in Eq. (14) as the average value of the waiting time. The parameter is given by Eq. (33). The parameter on the right hand side of Eq. (33) is related to the water depth of the lake. On the other hand, on the right hand side of Eq. (33) gives the variance of the jump lengths of a random walker, i.e. the variance related to how far it travels per unit time when it diffuses from the bottom of the lake (Eq. (11)). The parameter is proportional to and thus represents how fast individual Cs moves through the lake water. The vorticity diffusion coefficient () in the lake water is given in the literature from observations in lakes all over the world and is defined by Eq. (18). The parameter is the most important parameter and its value determines the overall behaviour of the diffusion. The value indicates the extent to which the diffusion differs from the classical behaviour. This can only be determined by fitting.

Results and discussion

In the following discussion, Cs activity concentration refers to the total Cs activity concentration (sum of dissolved and particulate Cs activity concentrations). The total Cs activity concentration is the average concentration of measurements at 0 m, 8 m and 15 m in the water column at the center of the lake[10]. Although total Cs activity concentration is very much variable and dependent on suspended matter concentration in the lake water, our main concern is the radiation safety of the lake water. That is why we use the total activity concentration in the present study.

The prediction formulae for Cs activity concentration used for comparison with the measured value are the fractional diffusion model (FDM) given by Eq. (39) and the two-component decay function model[10] (TDM) written as follows:where , , and are fitting parameters.

In order to determine the fitting parameters of Eqs. (39) and (40), the prediction formulae for Cs activity concentration of the water in Lake Onuma were fitted to the measured values. For all prediction formulae, we used the Levenberg–Marquardt algorithm to find the fitting parameters so that the squared relative error , defined by the following equation is minimized.where is the ith measurement value and is the fitting value at time .

The FDM reproduced well the temporal variation in the measured Cs activity concentration (Bq ) of the water in Lake Onuma collected in 2011–2016 (see Fig. 1). The solid line (FDM) in Fig. 1 is the result of fitting Eq. (39) with the unknown parameters (A, , and ) to the measured values. The values of the fitting parameters were , , and . The solid red circles indicate the measured Cs activity concentrations listed in the supplementary data[10]. The value of for the FDM was 0.207, and the of the TDM was 0.208 (Table 1). In order to reproduce the decay process of Cs activity concentration, the TDM required two exponential functions, while the FDM required only the Mittag-Leffler function. Since the fitting formula based on the FDM reproduced well the temporal change in the measured Cs activity concentration in Lake Onuma, Eq. (39) was used as a prediction formula for the long-term prediction of Cs activity concentration.

Figure 1

Change in the Cs activity concentration of water in Lake Onuma using the fractional diffusion model (FDM). The solid red circles indicate the measured Cs activity concentration found in the supplementary data[10].

Table 1

Values of the squared relative errors obtained as a result of fitting with each prediction models.

Model	ϵ^2
Fractional diffusion model (FDM)	0.207
Two-component decay function model (TDM)	0.208

On the other hand, a prediction formula based on the Bulgakov model[11] is as follows:where is a fitting parameter, is the radioactive Cs decay constant per day. Furthermore, a comparison between the FDM (Eq. (39)) and the Bulgakov model (Eq. (42)) showed that the FDM was a generalization of the Bulgakov model. Bulgakov et al.[11] derived the Cs activity concentration in the lake water as follows:where is the activity concentration of Cs in the lake water at time t, is the radionuclide deposition density, and is the average depth of the water column. The variable z in Eq. (43) is expressed as follows:here, is the effective diffusion coefficient in sediments, is the dimensionless distribution coefficient. Bulgakov et al.[11] derived Eq. (42) by approximating in Eq. (43) as follows:

The following relation holds between the Mittag-Leffler function and the function in Eq. (43)[47]:

Substituting in Eq. (46) into Eq. (43), we obtain the following equation:

Thus, the FDM is a generalization of the Bulgakov model[11]. In other words, in the Mittag-Leffler function in the Bulgakov model is 1/2, whereas the value of in the FDM (Eqs. (36) and (39)) is in the range of . For large time t the following relation holds[48]:

Therefore, Eq. (42) is reasonable for sufficiently large t, but when t becomes small, the approximation does not hold. Furthermore, in the FDM, the parameter of the Mittag-Leffler function has a degree of freedom . Therefore, the following discussion of the change in the Cs activity concentration of the lake water over time will be compared between the FDM and the TDM.

We fitted the FDM and the TDM to the measured values 234–1527 days after March 15, 2011, predicted the 137Cs activity concentration at 1527–1981 days, and then compared the predicted values and the measured values between the FDM and the TDM (see Fig. 2). The values of for the FDM and the TDM were 0.215 and 0.251, respectively and slightly better agreement was found for the FDM than for the TDM.

Figure 2

Comparison of blind test between the fractional diffusion model (FDM) and the two-component decay function model (TDM) using the measured data (234–1527 days after March 15, 2011) for the Cs activity concentrations of water in Lake Onuma. The solid red circles indicate the measured Cs activity concentration found in the supplementary data[10].

The fitting formula based on the TDM[10] consisted of the sum of two Mittag-Leffler functions with = 1. Noting that , Eq. (40) can be written using the two Mittag-Leffler functions as follows:

The FDM can reproduce the long-tailed distribution with respect to time t with a single Mittag-Leffler function. The TDM loses the degrees of freedom of because the value of is fixed to one. Therefore, at least two Mittag-Leffler functions must be used to reproduce the apparent long-tailed distribution.

A comparison of the FDM and the TDM predictions for Cs activity concentrations in the lake water up to 5000 days after the accident (see Fig. 3) suggested that the FDM prediction indicates a prolonged low-level radioactive contamination after 2000 days compared to the TDM prediction. The squared relative errors of the FDM and the TDM were 0.207 and 0.208, respectively, and there were no large differences in the prediction results up to 3000 days after the accident. The initial radioactivity of the lake water was 380 Bq  in the TDM prediction, while it was 1256 Bq  in the FDM prediction. Most of the Cs in Lake Onuma was remained within the uppermost layer (5 cm) of the deposits on the lake bottom. The inventory of Cs in the sediment has been estimated at 20 kBq per unit area () of the lake bottom[49,50]. We multiply the value, 20 kBq , by the bottom area of the lake[51], then divide it by the total volume of the lake water[51]. Thereby we obtain the estimation of the initial activity concentration 1960 Bq . Comparing initial value of our data (1960 Bq ) with the result of the FDM (A = 1256 Bq ), the predictions of the FDM at was in relatively good agreement. On the other hand, after 2000 days, the TDM prediction showed an exponential decrease in Cs activity concentration, while the FDM prediction showed a gradual decrease in Cs activity concentration. The difference between the TDM and the FDM is due to the dominance of the second term () in Eq. (40) after 2000 days. On the other hand, in the FDM prediction Eq. (48) holds when the time t becomes large and the Cs activity concentration follows the power law of t, so that the decrease of Cs becomes gradual.

Figure 3

Comparison between the fractional diffusion model (FDM) and the two-component decay function model (TDM) for changes in the Cs activity concentrations of water in Lake Onuma over 5000 days from March 15, 2011. The solid red circles indicate the measured Cs activity concentration found in the supplementary data[10].

To clarify the characteristics between the FDM and the TDM, we plotted the double logarithmic plot of both prediction results up to 10,000 days after the accident (see Fig. 4), and the FDM prediction suggested that the low-level Cs contamination was prolonged than the TDM prediction. The TDM prediction shows a rapid decrease in Cs activity concentration after 4000 days after the accident, while the FDM prediction values decreased linearly in double logarithmic plot until 10,000 days. Hence the prediction of the FDM decreased in inverse proportion to the power of time (), suggesting that the contamination of low-level Cs is prolonged. Although the measured data for Lake Onuma is only available until 1981 days after March 15, 2011, the comparison of the FDM with the measured data for Svyatoe lake in the Bryansk region of Russia[11], which is a closed lake, suggested that the FDM could be applied to the prediction of Lake Onuma 2000–5000 days after the Fukushima accident. Bulgakov et al.[11] compared the measured Cs activity concentrations of the water in Svyatoe lake during the period of about 2400–4900 days (about 7–13 years) after the Chernobyl accident with the results of calculations using Eq. (42) and showed a good agreement. Therefore, in order to verify the validity of the prediction of Cs activity concentration in Lake Onuma after 2000–5000 days after March 15, 2011 using the FDM, the fitting parameters obtained from Lake Onuma (, , ) were substituted into Eq. (39) and the fitting parameter A was fitted to the measured data of the water in Svyatoe lake[11] from 2400 to 4900 days and compared with the measured value (see Fig. 5). Both the predicted values based on the FDM and the values measured by Bulgakov et al.[11] decreased with the power law of time, and the two values agreed well. Therefore, the prediction of the FDM is considered to be valid for this period in Lake Onuma. The validity of the prediction for 5000–10,000 days after the accident needs to be evaluated by future measurements.

Figure 4

Double logarithmic plot for comparison between the fractional diffusion model (FDM) and the two-component decay function model (TDM) for changes in the Cs activity concentrations of water in Lake Onuma over 10,000 days from March 15, 2011. The solid red circles indicate the measured Cs activity concentration found in the supplementary data[10].

Figure 5

Double logarithmic plot for comparison between the measured values by Bulgakov et al.[11] and the fractional diffusion model (FDM) for changes in the Cs activity concentrations of water in Svyatoe lake in the Bryansk region of Russia days after the Chernobyl accident 2400–4900 days (7–13 years) . The solid line indicates the FDM using , , in Eq. (39) and red circles indicate the measured Cs activity concentration shown in Fig. 3 in Bulgakov et al.[11].

On the other hand, the correlation coefficient between the Cs activity concentration (Bq ) of the water in Lake Onuma and the Cs activity concentration (Bq  wet weight) of pond smelt (Hypomesus nipponensis), a typical fish species inhabiting the lake, was 0.95, indicating a strong positive correlation (Fig. 6). The activity concentration of Cs of pond smelt decreased rapidly from August 2011 (160 days after March 15, 2011) to September 2012 (552 days), but showed a gradual decreasing trend after October 2012 (592 days)[10]. This trend is consistent with the change in the Cs activity concentration of the lake water in Lake Onuma. If the correlation coefficient between the Cs activity concentration of the lake water and the Cs activity concentration of pond smelt is close to one, the Cs activity concentration of pond smelt can be predicted by applying the FDM. Therefore, the correlation coefficient between the Cs activity concentration of the lake water and the Cs activity concentration of pond smelt was calculated, and a strong correlation (correlation coefficient = 0.95) was found.

Figure 6

Correlations of the Cs activity concentrations between lake water and pond smelt in Lake Onuma. The calculated value of correlation coefficient was 0.95.

We compared the fitting results of Cs activity concentrations of pond smelt in the lake water using the FDM and the TDM, and found that both the FDM fitting (Fig. 7, ) and the TDM fitting () agreed with the measured values. In the following discussion, we compare the predictions of the FDM and the TDM.

Figure 7

Change in the Cs activity concentration of pond smelt in Lake Onuma using the fractional diffusion model (FDM). The solid red circles indicate the measured Cs activity concentration in pond smelt found in the supplementary data[10].

As shown in Fig. 8, there was no clear difference between the FDM and the TDM in the prediction of Cs activity concentration of pond smelt from 2000 to 4000 days, but the FDM prediction showed a gradual decrease after 4000 days. The reason for this difference is that the TDM decreases exponentially, while the FDM decays with the power law of time.

Figure 8

Comparison between the fractional diffusion model (FDM) and the two-component decay function model (TDM) for changes in the Cs activity concentrations of pond smelt in Lake Onuma over 5000 days from March 15, 2011. The solid red circles indicate the measured Cs activity concentration found in the supplementary data[10].

As shown in Fig. 9, double logarithmic plot for the predictions of Cs activity concentration of pond smelt up to 10,000 days showed that the FDM prediction resulted in a longer period of low-level Cs contamination than the TDM prediction. In the FDM prediction, when time t , the Cs activity concentration decreased according to the power law , whose index was . This indicated that decrease of the low-level Cs activity concentration of pond smelt was prolonged because the Cs activity concentration decreased gradually according to the power law of time. In other words, the results suggested that the low-level Cs contamination of pond smelt was also persistent due to the influence of Cs activity concentration in the lake water.

Figure 9

Double logarithmic plot for comparison between the fractional diffusion model (FDM) and the two-component decay function model (TDM) for changes in the Cs activity concentrations of pond smelt in Lake Onuma over 10,000 days from March 15, 2011. The solid red circles indicate the measured Cs activity concentration found in the supplementary data[10].

Conclusion

In conclusion, we showed that using the FDM for long-term prediction of Cs activity concentration in Lake Onuma, a closed lake, prolongs the low-level radioactive contamination compared to the TDM, a conventional fitting model. We derived the Mittag-Leffler function as a solution of the fractional diffusion equation based on the FDM in order to reproduce the gradual decrease of Cs activity concentration in closed lakes. This function has the property of exponential function when time t is small, and becomes long-tailed distribution when time t becomes large because it follows the power law . Using the FDM based formula (Mittag-Leffler function), it was possible to fit the temporal variation in the measured Cs activity concentration in the lake water without using formula consisting of two component exponential functions. The long-term prediction 5000 days after the accident using the fitting parameters to the measured data suggested that the FDM would result in a prolonged contamination with low activity concentrations of Cs compared to the TDM. Although the measured data for Lake Onuma is available until 1981 days after the accident, the comparison of the FDM calculation with the measured data for Svyatoe lake in Russia (about 2400–4900 days after the Chernobyl accident), which is a closed lake, suggested that the FDM could be applied to the prediction of Lake Onuma 2000–5000 days after the Fukushima accident. The validity of the prediction for 5000–10,000 days after the Fukushima accident needs to be evaluated by measurement in the future. The correlation coefficient between the Cs activity concentration of the lake water and the Cs activity concentration of pond smelt was calculated, and a strong positive correlation (correlation coefficient = 0.95) was found. Because the fitting with the measured Cs activity concentration of pond smelt showed a good agreement, we performed long-term prediction of Cs activity concentration in pond smelt using the FDM. A comparison of the FDM prediction with the TDM prediction suggested that the FDM prediction showed prolonged low-level Cs contamination. The prediction can be improved by applying the sum of three or more exponential functions instead of the sum of two exponential functions as in the TDM. Since no monitoring has been carried out since 2000 days after the accident, it makes sense to use the prediction formulae to give the residents living in the vicinity in Lake Onuma an idea of the future prospects of radioactive contamination. On the other hand, the monitoring of radiocesium is a research topic that should be continued in the future. To be summarized, a combination of the exponential functions to describe Cs dynamics in a lake for initial/short-term phase and the power law (diffusional) approach for mid-, and long-term is the merit of the proposed fitting model. The disadvantage of the FDM is that the fitting parameters of the FDM are less related to physicochemical processes in a lake and their characteristics than the models of Smith et al.[13] and Konoplev et al.[18,19].