Long-term prediction of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{137}$$\end{document}137Cs 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{137}$$\end{document}137Cs contamination because the activity concentration of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{137}$$\end{document}137Cs 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{137}$$\end{document}137Cs 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{137}$$\end{document}137Cs 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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{137}$$\end{document}137Cs than the TDM prediction.

www.nature.com/scientificreports/ 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 137 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][23][24][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][36][37] . Based on these backgrounds, we will try to reproduce the temporal variation in 137 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 137 Cs after that period. Since there are less measured values of 137 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 137 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 137 Cs activity concentration of the water and pond smelt in Lake Onuma based on the FDM, (2) conduct long-term prediction of 137 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 137 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 ( 137 Cs) and the waiting time between jumps are represented by a probability density function (PDF) ψ(x, t) . We consider the jump length PDF (x) and the waiting time PDF w(t) are independent variables, and one finds the separation of variable ψ(x, t) = (x)w(t).We describe the CTRW model using ψ(x, t) by the following master equation 39 : where η(x, t) 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 η(k, t) . We Laplace transform η(k, t) from time t to s and obtain η (k, s) . In this study, the Fourier transform and the inverse Fourier transform are defined, respectively, as and 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 ψ(x, t) in Eq. (1) can be written as follows: where ˜ (k) and ŵ(s) are the jump length PDF in the Fourier space and waiting time PDF in the Laplace space, respectively.
. www.nature.com/scientificreports/ Let �(t) be the survival probability that the waiting time on a site exceed t: The Laplace transform of �(t) corresponds to Equation (1) leads to the PDF of C(x, t) for 137 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 (x) behaves as a Gaussian jump length PDF: where σ 2 is given by We Fourier transform (x) and expanded to the order of k 2 , Eq. (10) can be approximated by the following formula: Suzuki et al. 10 reported the decay processes of 137 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 σ 2 τ α k 2 s α of the higher order 41

and obtain
where 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: www.nature.com/scientificreports/ where n is positive integer and n − 1 ≤ α < n . 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 ∂/∂t 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 a n depending on the initial conditions and use the separation of variables ansatz: Substituting Eq. (25)  where prefactors A c n and A sn should be determined by the boundary conditions. The solution to Eq. (29) is given through the Mittag-Leffler function E α (z) 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 www.nature.com/scientificreports/ Since our purpose is to make long-term predictions, the fundamental mode (n = 1) in Eq. (34) becomes dominant for a sufficiently large time t. Therefore, C(x, t) can be approximately written as follows: where ξ = ξ 1 . 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 d is the radioactive decay constant per day of 137 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 0 < α < 1 and −t α < 0 , then E α (−t α ) can be rewritten as: To calculate the Eq. (36), replace t in Eq. (38) with ξ 1/α (t/τ ) . 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 137 Cs (Bq m −3 ) 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, σ 2 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 σ 2 and thus represents how fast individual 137 Cs moves through the lake water. The vorticity diffusion coefficient ( K α ) 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, 137 Cs activity concentration refers to the total 137 Cs activity concentration (sum of dissolved and particulate 137 Cs activity concentrations). The total 137 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 137 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 137 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 Q 1 , k 1 , Q 2 and k 2 are fitting parameters.
In order to determine the fitting parameters of Eqs. (39) and (40), the prediction formulae for 137 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 ǫ 2 , defined by the following equation is minimized.
{A c n cos(µ n x) + A sn sin(µ n x)} T n (0) E α −ξ n t τ α .  10 . The value of ǫ 2 for the FDM was 0.207, and the ǫ 2 of the TDM was 0.208 (Table 1). In order to reproduce the decay process of 137 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 137 Cs activity concentration in Lake Onuma, Eq. (39) was used as a prediction formula for the long-term prediction of 137 Cs activity concentration.
On the other hand, a prediction formula based on the Bulgakov model 11 is as follows: where A B is a fitting parameter, d is the radioactive 137 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 137 Cs activity concentration in the lake water as follows: where C B (t) is the activity concentration of 137 Cs in the lake water at time t, σ d is the radionuclide deposition density, and h w is the average depth of the water column. The variable z in Eq. (43) is expressed as follows: here, D E is the effective diffusion coefficient in sediments, K d is the dimensionless distribution coefficient. Bulgakov et al. 11 derived Eq. (42) by approximating erfc(z) in Eq. (43) as follows: The following relation holds between the Mittag-Leffler function and the function exp(z 2 ) erfc(z) in Eq. (43) 47 :    (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 0 < α < 1 . 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 (0 < α < 1) . Therefore, the following discussion of the change in the 137 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 137 Cs 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 ǫ 2 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.
The fitting formula based on the TDM 10 consisted of the sum of two Mittag-Leffler functions with α = 1. Noting that E 1 (±z) = exp(±z) , 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 137 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 m −3 in the TDM prediction, while it was 1256 Bq m −3 in the FDM prediction. Most of the 137 Cs in Lake Onuma was remained within the uppermost layer (5 cm) of the deposits on the lake bottom. The inventory of 137 Cs in the sediment has been estimated at 20 kBq per unit area ( m 2 ) of the lake bottom 49,50 . We multiply the value, 20 kBq m −2 , by the bottom area of the lake 51 , then divide it by the total volume of the lake water 51   www.nature.com/scientificreports/ 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 137 Cs contamination was prolonged than the TDM prediction. The TDM prediction shows a rapid decrease in 137 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 ( t −α ), suggesting that the contamination of low-level 137 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 137 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 137 Cs activity concentration in Lake Onuma after 2000-5000 days after March 15, 2011 using the FDM, the fitting parameters obtained from Lake Onuma ( α = 0.62 , ξ = 5.00 , τ = 693 ) 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.  www.nature.com/scientificreports/ On the other hand, the correlation coefficient between the 137 Cs activity concentration (Bq m −3 ) of the water in Lake Onuma and the 137 Cs activity concentration (Bq kg −1 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 137 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 137 Cs activity concentration of the lake water in Lake Onuma. If the correlation coefficient between the 137 Cs activity concentration of the lake water and the 137 Cs activity concentration of pond smelt is close to one, the 137 Cs activity concentration of pond smelt can be predicted by applying the FDM. Therefore, the correlation coefficient between the 137 Cs activity concentration of the lake water and the 137 Cs activity concentration of pond smelt was calculated, and a strong correlation (correlation coefficient = 0.95) was found.  www.nature.com/scientificreports/ We compared the fitting results of 137 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, ǫ 2 = 5.60 ) and the TDM fitting ( ǫ 2 = 4.96 ) agreed with the measured values. In the following discussion, we compare the predictions of the FDM and the TDM.
As shown in Fig. 8, there was no clear difference between the FDM and the TDM in the prediction of 137 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.
As shown in Fig. 9, double logarithmic plot for the predictions of 137 Cs activity concentration of pond smelt up to 10,000 days showed that the FDM prediction resulted in a longer period of low-level 137 Cs contamination than the TDM prediction. In the FDM prediction, when time t ≫ τ , the 137 Cs activity concentration decreased according to the power law t −α , whose index was α = 0.76 . This indicated that decrease of the low-level 137 Cs activity concentration of pond smelt was prolonged because the 137 Cs activity concentration decreased gradually according to the power law of time. In other words, the results suggested that the low-level 137 Cs contamination of pond smelt was also persistent due to the influence of 137 Cs activity concentration in the lake water.

Conclusion
In conclusion, we showed that using the FDM for long-term prediction of 137 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 137 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 t −α . Using the FDM based formula (Mittag-Leffler function), it was possible to fit the temporal variation in the measured 137 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 137 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 137 Cs activity concentration of the lake water and the 137 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 137 Cs activity concentration of pond smelt showed a good agreement, we performed long-term prediction of 137 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 137 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 137 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 .