An inverse modeling method to assess the source term of the Fukushima Nuclear Power Plant accident using gamma dose rate observations

The Chernobyl nuclear accident, and more recently the Fukushima accident, highlighted that the largest source of error on consequences assessment is the source term, including the time evolution of the release rate and its distribution between radioisotopes. Inverse modeling methods, which combine environmental measurements and atmospheric dispersion models, have proven efficient in assessing source term due to an accidental situation (Gudiksen, 1989; Krysta and Bocquet, 2007; Stohl et al., 2012a; Winiarek et al., 2012). Most existing approaches are designed to use air sampling measurements (Winiarek et al., 2012) and some of them also use deposition measurements (Stohl et al., 2012a; Winiarek et al., 2014). Some studies have been performed to use dose rate measurements (Duranova et al., 1999; Astrup et al., 2004; Drews et al., 2004; Tsiouri et al., 2012) but none of the developed methods were carried out to assess the complex source term of a real accident situation like the Fukushima accident. However, dose rate measurements are generated by the most widespread measurement system, and in the event of a nuclear accident, these data constitute the main source of measurements of the plume and radioactive fallout during releases. This paper proposes a method to use dose rate measurements as part of an inverse modeling approach to assess source terms. The method is proven efficient and reliable when applied to the accident at the Fukushima Daiichi Nuclear Power Plant (FD-NPP). The emissions for the eight main isotopes 133Xe, 134Cs, 136Cs, 137Cs, 137mBa, 131I, 132I and 132Te have been assessed. Accordingly, 105.9 PBq of 131I, 35.8 PBq of132I, 15.5 PBq of137Cs and 12 134 PBq of noble gases were released. The events at FD-NPP (such as venting, explosions, etc.) known to have caused atmospheric releases are well identified in the retrieved source term. The estimated source term is validated by comparing simulations of atmospheric dispersion and deposition with environmental observations. In total, it was found that for 80 % of the measurements, simulated and observed dose rates agreed within a factor of 2. Changes in dose rates over time have been overall properly reconstructed, especially in the most contaminated areas to the northwest and south of the FD-NPP. A comparison with observed atmospheric activity concentration and surface deposition shows that the emissions of caesiums and 131I are realistic but that132I and132Te are probably underestimated and noble gases are likely overestimated. Finally, an important outcome of this study is that the method proved to be perfectly suited to emergency management and could contribute to improve emergency response in the event of a nuclear accident.


Introduction
The Great East Japan Earthquake followed by a tsunami that occurred in Japan in March 2011 led to the most significant nuclear accident since the one that occurred at Chernobyl in 1986.To assess the consequences of such accidents on humans and the environment, it is necessary to know the temporal evolution of the release quantities of all emitted radionuclides.However, previous accidental situations have shown that estimating releases may be a difficult task due to strong uncertainties.For example, the source term related to the Chernobyl accident remains uncertain (Devell et al., 1995).In the case of the accident at FD-NPP, incomplete knowledge of the state of the facility made the assessment of releases particularly difficult.
To reduce uncertainty and improve assessment of atmospheric releases, one solution consists of estimating the source term by combining measurements of the environment with atmospheric dispersion models.The following measurements can be carried out in the event of a nuclear accident: atmospheric activity concentrations, surface activities and gamma dose rates.Measurements of activity concentration provide the concentration of each radionuclide in the air.While the measurement is accurate, it is, however, averaged over a period of time, usually on the basis of several hours or even daily.Measurement of surface activity quantifies fallout for each radionuclide.The measurement is integrated over periods of time and in the case of Fukushima, for example, they are daily at best.Lastly, ambient gamma dose rate measurements sum the direct contribution of the plume (plume radiation) and the gamma radiation emitted by radionuclides that fell to the ground (deposited radiation) through dry and wet deposition processes.The data interpretation is complex since the signal provides no information about isotopic composition or the respective contributions of the plume and deposition.Nevertheless, networks of gamma dose rate stations are generally dense and the devices offer a high temporal resolution.
Among the methods developed to estimate source term from a nuclear accident by using environmental measurements, certain approaches use all types of measurements.This is the case for simple approaches for source estimation such as that proposed by Chino et al. (2011).These methods consider each monitoring location point independently of the others.On an empirical basis, they estimate emissions that may explain the signal observed at each monitoring location.It is the sum of all emissions that constitutes the source term.The disadvantage of these approaches results from the algorithm not taking into account errors from the measurements, from the dispersion model, and from the meteorological data that force the model.While these errors cannot be overlooked, the measurement points can no longer be independently considered from one another and agreement between simulations and all available observations must be optimized.
Methods based on inverse modeling have the advantage of taking into account the various errors.Since the Chernobyl disaster, inverse methods have been developed to estimate source term from measurements of atmospheric activity concentration and surface activity (Gudiksen et al., 1989;Krysta and Bocquet, 2007;Davoine and Bocquet, 2007;Winiarek et al., 2011;Stohl et al., 2012a).Assimilation methods have been developed to enhance the prediction of gamma dose field by using synthetic dose rate data (Astrup et al., 2004).No source inversion has been done with real gamma dose rate data.
In the case of the Fukushima accident, simple methods for source estimation and inverse methods have been applied.The estimates of emissions obtained using simple methods were first proposed by Chino et al. (2011) and later reassessed by Katata et al. (2012) and Terada et al. (2012).Their approach is based on the use of activity concentrations and surface activity measurements and part of the dose rate measurements.The use of dose rate measurements requires assuming the isotopic composition of the release.
The Mathieu et al. (2012) and Korsakissok et al. (2013) studies are also based on a simple method for source estimation (IRSN, 2012).The isotopic composition of the release is based on the standard radionuclide inventory of a reactor.It evolves over time depending on the quantities that have already been emitted and the volatility of the radionuclides (D.Corbin and J. Denis, personal communication, 2012).Emission peaks and released quantities are assessed based on the detailed analysis of events that have occurred at the facility and measurement of ambient dose rates.Stohl et al. (2012a) and Winiarek et al. (2012) propose estimates of source terms obtained by inverse modeling.In both cases, measurements of activity concentration are used.Stohl et al. (2012a) and Winiarek et al. (2014) also use measurements of surface activity.Winiarek et al. (2012) demonstrates that it is not possible to obtain complete source term using only this very limited number of measurements available in the Honshu island, which are mainly distributed south of FD-NPP.Stohl et al. (2012a) and Winiarek et al. (2012) used a prior knowledge of the emissions (first guess) to compensate for the lack of measurements.Inverse modeling was then used to improve the first guess.
For some accidental situations, it may be tricky to quickly estimate a reliable first guess.In such cases, it is relevant to estimate the source term by using environmental measurements without having to construct a specific a priori source term.Given the configuration of measurement networks in countries with nuclear facilities, it is highly likely that, as in the case of the Fukushima accident, measurements of activity concentration and surface activity are insufficient and provide data which are not acquired with adequate frequency.It is thus crucial to be able to use a network of dose rate stations, which are usually dense.This paper proposes a method based on inverse modeling techniques to estimate source term by using dose rate measurements.The approach can be applied with or without a specific first guess.The method is built to be used in the management of an accidental situation.
Section 2 describes the characteristics of the dose rate signal and the method used.In Section 3, the method is applied to the Fukushima accident and the source term is estimated based on dose rate measurements.The realism of (2) shows the signal behavior when the plume moves further away.
(3) shows the dose rate signal due to the deposition.
the assessed source term is discussed in Section 4. First, it is compared with source terms that have already been published.The source term is then used to simulate atmospheric dispersion of the FD-NPP releases and the results are compared with dose rate measurements, activity concentration and surface activity.Section 5 provides conclusions and perspectives about this study.

Inverse modeling of accidental release by using dose rate measurements
This section presents the generic method developed to estimate the source term due to a nuclear accident by using dose rate measurements.The implementation of the method to the Fukushima accident is described in Sect.3.

Dose rate measurements
The variable "ambient dose rate" includes various information.It is the sum of the contributions of all gamma-emitting radionuclides, whether airborne or fallout: with D ambient the ambient dose rate; D plume n the contribution of radionuclide n to dose rate from plume radiation; D deposit n the contribution of radionuclide n to dose rate due to deposited radiation and M the number of gamma emitters.The temporal evolution of the ambient dose rate signal measured when a contaminated air mass passes over a monitoring station is given in TEPCO ( 2012).An idealized situation is shown in Fig. 1.
When the plume blows above the monitoring station, the dose rate level rises and reaches a peak (e.g.event 1 on Fig. 1).This signal is due to radioactivity in the air and the deposition that gradually forms in the vicinity of the device.
As the plume moves further away, the dose rate decreases (e.g.event 2 on Fig. 1).After the plume has passed, the dose rate level is higher than it was before the plume (e.g.event 3 on Fig. 1).The difference results from the deposits that formed during the plume passage.At this point, the temporal evolution of dose rate is mainly due to radioactive decay of the radionuclides that constitute the deposit.The dose rate decreases over time and the signal slope results from the isotopic composition of the deposit.

Inverse modeling of accidental release using dose rate measurements
The objective of the method is to estimate the release rate for each radionuclide n using the measurements of ambient dose rate given that the number M of radionuclides contributing to the dose rate signal may be significant (more than a hundred in the case of a release from a nuclear power plant).We consider that the source location is known.The reconstruction of the source term can be expressed as the inverse problem described by the source-receptor equation below.The equation relates the measurements to the source term using an atmospheric dispersion model: In our case, µ ∈ R d is the vector of dose rate measurements where d is the number of measurements.σ ∈ R N is the source term vector where N is the number of radionuclides multiplied by the number of time steps.σ includes the temporal evolution of the release rates for each radionuclide.ε ∈ R d is the error vector.It represents a combined modelrepresentativeness and instrumental error.H ∈ R d×N is the Jacobian matrix that contains the results of the dispersion model at the measurement locations.
Since the Jacobian matrix is ill-conditioned, the source receptor equation constitutes an ill-posed inverse problem (Enting, 2002) and cannot be solved directly.To compensate for this problem, a Tikhonov regularization is added to ensure the existence and unicity of the solution (Tikhonov and Glasko, 1965).Resolution of the inverse problem can then be summed up by minimizing the following cost function: When the dose rate network is sufficiently dense and welldistributed, it is possible to obtain a good solution without having to construct a specific a priori and σ b = 0 is used.
Solving the inverse problem by using dose rate measurements and a σ b = 0 requires to reduce the size of the problem, i.e., reduce the number of parameters to solve and limit the solution space.To achieve this objective, we propose a method, which requires proceeding by steps.
-The first step is to define the a priori information required; O. Saunier et al.: Fukushima Nuclear Power Plant accident -the second step is to identify the potential release periods; -the third step is to estimate the release rates during the periods identified during the second step.

Step 1: defining the a priori information
For an accident involving a nuclear power plant, the main part of the dose rate signal is due to few radionuclides.The relevant radionuclides can be identified by using other types of measurements made in the environment or by using the core inventory of the damaged facility.The number of unknown parameters is decreased by reducing the composition of the source term to this list of radionuclides.
To solve the inverse problem, the respective contributions of the principal radionuclides must be distinguished using only information contained in the dose rate measurements.The temporal evolution of signal due to the radioactive decay of the deposit contains indirect information on the isotopic composition of the emissions.However, it is necessary that the half-lives of the selected radionuclides be sufficiently different so that the inversion process can distinguish their respective contribution to the dose rate signal.On a physical point of view, the various radionuclides are released simultaneously in proportions that depend on their physicochemical properties and the core inventory, among other things.The inverse problem can be constrained by imposing that the radionuclides be released in realistic proportions.For example, the release rates of radionuclides i in relation with radionuclide 1 can be made to respect at each time step t, the following proportions: with a i and b i the limit values of the isotopic relationships.a i and b i can be estimated using the core inventory of the damaged reactor and the measurements of activity concentration and surface activity.The constraints may be more or less flexible depending on the emissions knowledge.

Step 2: identification of potential release periods
The purpose of this second step is to determine when the radionuclides responsible for the increase of dose rate during the passage of a plume over the monitoring stations were emitted.Inverse modeling techniques can be used to solve this problem following the subsequent process.
It is assumed that the release is composed of a single radionuclide which acts like a passive tracer.The vector of measurements, µ plume , is built with dose rate measurements as follows.An algorithm is used to analyze the slope in the dose rate signal and the peaks corresponding to the detection of the passage of plume are extracted.An example is given in Fig. 2. The dose rate measured at the Tsukuriya location during the Fukushima accident is shown in black and "the extraction of the plume component" is in red.Extraction process is performed on data from all the monitoring locations and the "plume component extraction" constitutes the vector of measurements, µ plume .
A specific first guess is not required and σ b = 0.The cost function of Eq. (3) becomes the following: Minimization of J 1 provides an estimate of vector σ 1 .The potential release periods are encapsulated in vector σ plume based on vector σ 1 according to the following rules: with t ∈ [0; N t ], N t is the number of time steps and σ plume ∈ R N t .The second step of the proposed method brings the number of unknown parameters to the number of time steps corresponding to potential release periods multiplied by the number of radionuclides selected.It remains to estimate the release rates during the potential emission periods.

Step 3: estimating source term
The purpose of the third step is to solve the inverse problem, i.e., determine the source term, σ , composed of the temporal evolution of the release rates for each radionuclide.The release periods must be among those identified during the second step.The M selec radionuclides that constitute the source term are selected during the first step as being the main contributors to the dose rate signal.Finally, the source term, σ , must comply with the constraints on isotopic relationships defined during the first step.The cost function of Eq. (3) becomes µ contains all the dose rate measurements.M selec − 1 is the number of constraints involving the isotopic composition of the source term.r i is the function that constrains the isotopic relationship of radionuclides i and in relation with radionuclide 1.It is defined as a barrier-protection function in order to comply with the constraints: The M selec − 1 constraints are used to force the method to release the different radionuclides simultaneously.

Towards an implementation of the method
In summary, the implementation of the method requires a set of a priori information: -the location of the source; -the selected radionuclides identified as responsible of the main part of the dose rate signal; -the limit values a i and b i of the isotopic ratios (Eq.4).
When a nonzero first guess is used assuming it includes all the release periods, the second step of the method becomes unnecessary and can be avoided.When no first guess is used, the release events reconstructed by the inverse method are only those that have been measured.Before being applied to a real accident situation, the method was tested and validated on a synthetic case (results not shown).The implementation on the FD-NPP accident is presented in the next section.

Application: the Fukushima accident
The objective consists of estimating the Fukushima source term of atmospheric radiation emissions at hourly intervals between 11 March 2011 at 9 a.m.(JST) and 27 March at 6 p.m. (JST), for a total of 381 time intervals.

Model simulations
We use the Eulerian ldX model to simulate the radionuclide dispersion.This model is part of IRSN's (French Institute for Radiation protection and Nuclear Safety) C3X operational platform (Isnard, 2006).It is based on the Polair3D chemistry transport model (Boutahar et al., 2004), which is itself part of the Polyphemus system, and has been validated on the Algeciras incident as well as on the Chernobyl accident (Quelo et al., 2007).It takes into account the dry deposition and the wet deposition as well as the decay products and the radioactive decay.The dry deposition is modeled using a simple scheme with a constant deposition velocity: The wet scavenging is parameterized in the form s = 0 p 0 , where 0 = 5 × 10 −5 h (mm −1 s −1 ) and p 0 is the rain intensity in mm h −1 (Baklanov and Sørensen, 2001).The vertical diffusion follows the Louis scheme (Louis, 1979) in stable meteorological situations, and the Troen-Mahrt scheme (Troen and Mahrt, 1986) if the atmosphere is unstable.
The Fukushima accident simulations are performed by forcing the model with three-hourly operational meteorological data from the European Center for Medium-Range Weather Forecasts (ECMWF).The spatial resolution of these data is 0.125 • × 0.125 • .At this resolution, Japan's complex orography cannot be fully resolved, and Mathieu et al. (2012) and Korsakissok et al. (2013) state that some situations are poorly reproduced by the forecast model.
The spatial domain of the ldX simulations is [135 , which encompasses Japan.The time resolution is one hour and the model provide instantaneous outputs.The spatial resolution is that of the meteorological data (0.125 • × 0.125 • ), and 11 vertical levels between 0 and 3400 m are considered.The release height is taken for the first two levels of the model, i.e. between 0 and 160 m.
Lastly, the dose rate at each measurement point is computed from the activity concentrations and surface activities simulated by ldX.This is done with the C3X platform's ConsX model.
The Jacobian matrix, H, is computed under the approach proposed by Abida et al. (2009) and subsequently adopted by Winiarek et al. (2011).Each column of H represents the dispersion model's response to a unitary release emitted for each radionuclide whose release rate is to be estimated.In total, 381 simulations are carried out in order to construct the H matrix.

Step 1: defining the a priori information
A total of 57 dose rate monitoring stations located in Japan are considered in the inversion process (see Fig. 3).Only approximately thirty stations recorded the temporal evolution of the dose rate at least hourly from the beginning to the end of the accident.The other stations' data are incomplete, as their measurements generally begin after 15 March.
The large majority of the available measurements were obtained by the automatic detectors in the SPEEDI (System for Prediction of Environment Emergency Dose Information) network.Other measurements were taken for each prefecture by the Japanese Ministry of Education, Culture, Sport, Science and Technology (MEXT) (MEXT, 2011a).Other measurements are available in different prefectures (Tochigi Prefecture, 2011;Miyagi Prefecture, 2011;KEK, 2011;Fukushima Prefecture, 2011a, b).Most of those observations are now gathered in the International Atomic Energy Agency Fukushima Monitoring database (IAEA, 2012).
The measurements include those of the Chiba monitoring station (Nagaoka et al., 2012), which is the only station that measured the dose rate of each radionuclide.This station provides considerable useful information, as it detected many releases during the accident.The radionuclides measured are therefore fairly representative of the radionuclides released during the FD-NPP accident.The measurements indicate that most of the dose rate signal is due to 134 Cs, 137 Cs, 136 Cs, 137m Ba, 131 I, 132 I, 132 Te, and 133 Xe.We therefore try to assess the release rate of these eight radionuclides.
All of the noble gases emitted during the accident are grouped and are estimated as 133 Xe emissions.By definition, noble gases do not deposit and they only contribute to the radioactive plume.By their very nature, noble gases are quickly emitted, and it is safe to assume that they were no longer present in the emissions after a certain date.In this study, we assume that all of the noble gases were emitted before 21 March.
The elements of the 137m Ba/ 137 Cs pair are in secular equilibrium. 132I (half-life 2.3 h) is the son of 132 Te (half-life 3.26 d).More than one day after the shutdown (i.e. a period of time greater than ten times the half-life period of the 132 I element) the 132 Te/ 132 I pair can be considered in equilibrium.At each instant (t), the release rates of the two pairs obey the following rules: The half-lives of 137 Cs, 134 Cs, 131 I, 136 Cs, and 132 Te (respectively, 30.17 yr, 2.06 yr, 8.02 days, 13.2 days and 3.26 days) are sufficiently different to have separate impact on the deposition dose rate signal.The simulations do not last long enough (simulation run of 15 days) to detect the impact of the decrease of 134 Cs and 137 Cs concentrations on the dose rate signal.The time window is indeed too short so that 137 Cs and 134 Cs can be differentiated by using the deposition dose rate measurements.On the other hand, an analysis of the activity concentration measurements for the whole of Japan has shown that the ratio between the 137 Cs and the 134 Cs was constant over the time, and approximately equal to 0.94.The following isotopic ratio was therefore required in order to calculate the emissions: Ultimately, it is only necessary to explicitly assess the release rate of five radionuclides.Finally, flexible constraints have been defined for the isotopic ratios by comparing a typical inventory of a nuclear plant core of the same type as the FD-NPP with the activity concentration and dose rate measurements at Chiba (Corbin and Denis, personal communication).In the cost function, we specified that the emissions must meet at each instant (t), the following constraints:

Step 2: identification of potential release periods
The method's second step, consisting of determining the potential release periods, is applied.The error covariance matrices B 1 and R 1 are chosen by using a simple parameterization (Winiarek et al., 2011).It is assumed that the two matrices are diagonal, and the error variance is the same for all diagonal elements of each matrix (homoscedasticity property): m 1 , k 1 are the diagonal elements of the matrixes B 1 and R 1 .I d denotes the identity matrix.
Based on these assumptions, the cost function J 1 to be minimized is reduced to the following expression: , where: A threshold can be applied on σ 1 to reduce the CPU time of step 3 and not have any impact on the source term assessment.In our application, we used σ 1 (t) > 10 8 for that purpose.
As the parameter λ 1 determines the scale of the source term's fluctuations, it defines how much importance should be given to minimizing the source term's norm.The errors relating to the measurements and the model are contained in the parameterization of λ 1 .This parameter can be automatically determined by means of hyperparameter techniques 46.9 % 114 (Davoine and Bocquet, 2007;Krysta et al., 2008;Winiarek et al., 2012).In our case, it is set directly by the user.
The L-BFGS-B algorithm (Liu and Nocedal, 1989) is used to minimize the cost function J 1 for a wide range of values of λ 1 , of which Table 1 presents a sample.L-BFGS-B is a limited-memory quasi-Newton algorithm for boundconstrained optimization.The positivity of the source vector is assumed and no upper bound is used.The second column of the table shows the relative error computed using the values of the cost function J 1 before and after minimization.The relative error is sensitive to the values of λ 1 : the smaller the value of λ 1 , the lower the relative error.The value of the parameter λ 1 has a smaller influence upon the number of estimated potential release instants.Choosing the value of λ 1 involves a compromise between reducing the error value and increasing the number of potential release instants.As a result, λ 1 = 10 −12 was chosen as the value.The second step therefore enabled us to reduce the 381 initial instants to 120 potential release instants.

Step 3: estimating source term
The third step is applied in order to assess the release rates from the dose rate measurements, assuming σ b = 0.The same error covariance matrix parameterization is used as in step 2. The cost function J 2 is then written as follows: The homoscedasticity property of the parameterization leads to give more weight on the high values than on the low values, as it is the absolute value of the model-to-measurement difference that is considered in the cost function to be minimized.
The sensitivity to the parameter λ 2 has been tested for a wide range of values (0 ≤ λ 2 ≤ 10 −8 ).The results obtained are relatively insensitive to the value chosen.The method converges towards the same value of J 2 and the same source term.Several reasons help to reduce the influence of the regularization term.First, many observations are used in the inversion process.Second, steps 1, 2 as well as the use of a barrier-protection function for step 3 reduce the solution space of the inverse problem.Third, the zero initial condition contributes to the regularization of the problem.In practice, this allows the use of very low values of λ 2 .The results shown hereafter are those obtained from λ 2 = 10 −12 .

Source term
The quantities emitted in the source term estimated by our inverse modeling approach (hereafter called retrieved ST) are shown in Table 2 and compared with the other assessments.The retrieved ST's total emissions in 136 Cs, 137 Cs, 131 I, 132 I and 133 Xe are 3.8 × 10 15 Bq, 1.55 × 10 16 Bq, 1.059 × 10 17 Bq, 3.58 × 10 16 and 1.2134 × 10 19 Bq, respectively, which is consistent with the other estimates.It should be noted that any events not detected by the dose rate monitoring stations, such as releases transported directly to the Pacific Ocean, could not be reconstructed correctly.The estimated quantities therefore provide a lower bound of FD-NPP emissions.However, like the results of Stohl et al. (2012a), the 133 Xe quantities obtained are slightly higher than the noble gases inventory of the three damaged reactors.The xenon releases are probably overestimated.As suggested by Stohl et al. (2012b), the releases of short-lived radionuclides as well as the 133 I emissions, not taking into account in the retrieved ST, may explain the probable overestimation of xenon releases.
Figure 4 represents the temporal evolution of 137 Cs release rates resulting from various estimates.Until 20 March, the retrieved ST's release periods are consistent with the other assessments.Table 3 compares the chronology of the events listed for the nuclear facility (TEPCO, 2012; Corbin and Denis, 2012) with our release periods.It shows that our approach makes possible to find all nuclear facility events except two ventings on 17 and 18 March, when the resulting releases were not detected in Japan because they were transported directly towards the Pacific Ocean.
From 20 to 27 March, the retrieved ST shows additional release peaks compared to the other estimates.These releases may not be related to specific events in the nuclear facility, but the measurements clearly indicate that releases occurred.The variations in the water injection intensity of reactors 2 and 3, combined with the pressure variations and the fact that smoke was observed on the site could explain these emissions.
The method is therefore effective in reconstructing the nuclear facility events.This result is very important because the source term is obtained without requiring a prior knowledge of the emissions.
Figure 4 shows that the different estimated release rates vary considerably, illustrating how difficult it is to reconstruct   3).
the source term.All of the source term estimation methods inherit many uncertainties arising from the number of measurements used, the quality of the meteorological data, the realism of the first guess if the method uses one, and the quality of the atmospheric dispersion and deposition simulations.
The main limiting criteria affecting our method are the number of monitoring stations that detected each release event, and the quality of the meteorological data.Although a release event can be reconstructed if the plume is detected by only one or two detectors, the inverse problem is less constrained than if the release were detected by more monitoring stations.Furthermore, errors in the wind or precipitation data can bias the atmospheric dispersion and deposition simulations, thereby preventing the signal observed at the monitoring stations from being reproduced.In this case, the source term estimated by the inverse modeling approach is a compromise between the simulations and the observations.The quality index associated with each event (5th column of Table 3) specifies the origin of the uncertainties that might bias the release estimations.The index varies between 0 and 3.An index of 0 corresponds to events that are not reconstructed.An index of 1 corresponds to events that have been reconstructed using too few measurements in number, too distant or too small to have sufficient weight in the source term estimation.An index of 2 corresponds to situations in which there are sufficient measurements but the meteorological conditions are poorly reproduced, possibly resulting in inconsistencies between the simulations and the observations.Lastly, an index of 3 indicates events that have been reconstructed properly.Very few events have an index of 1 or less.The weather forecasts are often insufficiently accurate, however, and many episodes have an index of 2. It is why the release rates of the retrieved ST must be validated by comparing it with the environmental measurements.This is the subject of the following section.

Comparison with observations
The source term is validated by comparing the environmental measurements with the atmospheric dispersion simulation obtained by forcing the model with the retrieved ST.The objective of the model to data comparisons is to check whether all of the release events are reconstructed properly and in the right timescale, and whether the quantities and the isotopic composition of the releases are accurate.
Of the source terms considered in Table 2, only the retrieved ST and the Mathieu et al. (2012) ST afford the release rate of the required radionuclides to simulate dose rate.The Mathieu et al. (2012) ST was assessed using the same dose rate measurements and the same meteorological data as those we used in our study.Mathieu et al. (2012) and Korsakissok et al. (2013) have shown that the Mathieu et al. ( 2012) ST was realistic and highly comparable to the observations.If our inversion method is relevant, Table 3. Main facility events that caused environmental releases (TEPCO, 2012;Corbin and Denis, personal communication).Comparison with the release periods identified by means of inverse modeling.The 4th column indicates the number of detectors used to reconstruct the source term.The 5th column provides a reconstruction quality index.Index of 0: events not reconstructed.Index of 1: events reconstructed using too few measurements in number, too distant or too small.Index of 2: events reconstructed using sufficient measurements but using meteorological conditions poorly reproduced.Index of 3: events properly reconstructed.The 6th column indicates the release trajectory.

Comparison with the dose rate measurements
The agreement obtained between the simulations and the dose rate observations is summarized in Table 4. Fac2 represents the proportion of the simulated dose rates that is within a factor of 2 of the observed values.The relative bias is given by where dm is the simulated dose rate, do is the observed dose rate, and d is the number of dose rate measurements.On average, the simulations performed with the retrieved ST are in good agreement with the observations.For 80 % of the measurements, observed and modeled values agreed within a factor of 2. The 0.4 bias value is low.The proportion of the observations where simulated and observed dose rates agreed within a factor of 2 is improved by 20 % and the bias is reduced by 29 % compared with the values obtained with the Mathieu et al. (2012) ST.
In order to analyze the results in greater detail, Japan was then divided into four zones located to the north, the northwest, the west and the south of the FD-NPP.Figure 3 shows the positioning of the zones.The zones are defined on the basis of the dose rate signals' and the detected release events' consistency: -The North zone contains few monitoring stations and the observations are only available from 17 March onwards, with the exception of the Minami Soma station, which was operational from 12 March onwards.
-The North-West zone is the most highly contaminated area, and includes the cities of Aizuwakamatsu, Fukushima, Yonezawa and Nihonmatsu.
-The West zone is located at the edge of the most heavily contaminated areas.The area is crossed from north to south by a range of mountains, with the towns of Tadami and Minamiaizu to the west and Shirakawa to the east.
-The South zone is the second important area, as it has been affected by many release events and includes major Japanese cities, including Tokyo.It consists of the stations located in the Ibaraki prefecture (approximately 100 km to the south of the nuclear facility), on the one hand, and stations further to the south-west, such as Saitama, Chiba and Tokyo, on the other hand.
Table 4 shows the scores achieved by the simulations for each zone.For the two most "emblematic" zones (North-West and South), the retrieved ST provides simulations that are highly consistent with the measurements and greatly improves upon the Mathieu et al. (2012) ST scores.The scores are significantly lower for the West zone (fac2 = 45 % and = 0.75), even though the retrieved ST improves them.In the North zone, the simulations are moderately consistent with the measurements.The fairly good proportion of observations where simulated and observed dose rates agreed within a factor of 2 (fac2 = 64 %) do not compensate the significant bias of 0.7 that the retrieved ST fails to improve compared with the results obtained with the Mathieu et al. (2012) ST.
Figure 5 shows examples, for each zone, of the comparison of the observed and simulated dose rates for monitoring stations for which the modeled values are consistent with the measurements.Some examples of model-measurement inconsistencies are shown in Fig. 6 for the North and West zones.
Column 7 of Table 3 indicates which zones are affected by each release event.Event # 6 stands out because it affected all four zones in turn.This is the episode responsible both for the main contamination of the Japanese territory and for the moderately good scores obtained in the North and West zones.The releases firstly drifted southwards, then westwards and north-westwards where the plume was washedout.Next, they were carried northwards then southwards again, and then westwards.The second drift westwards was accompanied by wet scavenging.
From 15 March onwards, the simulations significantly underestimate the dose rates of the North zone monitoring stations close to the Pacific coast (e.g.Fig. 6a) and overestimate those of the stations further inland (e.g.Fig. 6b).Design the source term to obtain simulated values that are more consistent with the measurements for the monitoring stations near to the coast inevitably has the opposite result for the inland stations (or vice versa), as the wind data do not allow the simulated plume to drift sufficiently towards the north.
The comparisons for the eastern part of the West zone (upstream of the central mountain range crossing the island of Honshu from north to south) show that the dose rate signal is generally well simulated with the retrieved ST.The results are much less satisfactory in the case of the stations located on the foothills of the range, however.The fact that the results for the stations located on both sides of the mountain area (see Fig. 6, Minamiaizu and Shirakawa stations) are incompatible is due to the distance and time lag affecting the precipitation and wind data, which is probably because the effects of the complex terrain are not fully taken into account.Release event # 6 can nevertheless be considered as being properly reconstructed by the retrieved ST, as the stations in the North-West and South zones (see Fig. 5) constrain it sufficiently and marginalize the impact of the few incompatible stations in the West and North zones.
Figure 5 shows that, on average, the dose rate signal's temporal evolution is accurately reproduced by the simulations.Most of the peaks are simulated at the correct times, and their intensity is perfectly consistent with the observations.Some secondary peaks are occasionally poorly reproduced, but such events are marginal.As a result, these comparisons validate the source term in terms of the release quantities and the temporal evolution of the emissions, leaving us to assess the method's suitability for reconstructing the releases' isotopic composition.This is done by comparing the temporal evolution of the dose rate signal due to deposition irradiation when the simulated and observed values are of a similar magnitude.The signal's slope is mainly due to the radioactive decay of the deposited species.The first-order simulated and observed slopes are highly consistent.Nevertheless, we can see that the simulations' slope is slightly less pronounced than that of the actual observations.These findings suggest that the retrieved ST's composition is generally realistic but perhaps slightly underestimates the proportion of radionuclides with the shortest half-lives, such as 132 I and 132 Te.In order to verify this hypothesis and complete the validation of the source term, the simulations were compared with the activity concentration and surface activity observations.

Comparison with the atmospheric activity concentrations measurements
The activity concentration measurements were not used to estimate the Mathieu et al. (2012) ST and the retrieved ST.As a result, comparing the simulations with these observations is an excellent way of validating the source terms, albeit only partially, as most of the few measurements available are from the south zone (MEXT, 2011b).Figure 7 provides an example of a comparison obtained with the As regards the released quantities and the emissions' isotopic composition, the degree to which the modeled value is consistent with the measurements is quite encouraging.For the main episode on 15 and 16 March, the Tokai monitoring station detected the main peak over a longer period in the simulations than in the actual measurements.The simulation reproduces correctly the main peak on 15 March.For the rest of the episode, the simulation results are fairly consistent with the measurements in the case of the caesiums.The 136 Cs activity concentrations are slightly overestimated by a factor of 2 to 3, and the simulation results are most consistent with the measurements in the case of 137 Cs.The 131 I activity concentrations are overestimated but remain within a factor of 2 of the observations.The results are less conclusive in the case of the 132 Te, as the activity concentrations are underestimated by a factor 3 to 4 in the simulations.The 132 Te/ 132 I pair has a relatively short half-life of 3.26 days.If several plumes are detected above the same monitoring station for a short period, the dose rate due to deposition between two successive plume passes is not measured for a long enough period to enable the inverse method to accurately estimate the releases' isotopic composition.In this context, the radioactive elements with relatively short half-lives, such as those of the 132 Te/ 132 I pair, are considered as noble gases by the inversion process and the method will tend to emphasize the 133 Xe in the release rather than the 132 Te/ 132 I pair.Two plumes are detected above the Tokai monitoring station on 15 and 16 March.The underestimated 132 Te value should therefore be compensated for by overestimating the 133 Xe value.This hypothesis cannot be checked, as we do not have simultaneous noble gas and 132 Te measurements, but all of the comparisons tend to support it.
The simulation made using the Mathieu et al. (2012) ST tends to overestimate the activity concentrations, resulting in the simulated dose rate being overestimated after 15 March due to deposition during the event.
Table 5.Comparison of the observed and simulated daily depositions using the retrieved ST and the Mathieu et al. (2012) ST.The statistics are the proportion of the measurements for which observed and modeled values agreed within a factor of 2, 5 and 10.

Radionuclide
Retrieved ST Mathieu et al. (2012) ST Fac2 (%) Fac5 (%) Fac10 (%) Fac2 (%) Fac5 (%) Fac10 (%) 136  No comparison can be done for the second episode between 20 and 23 March because the simulated event is slightly late, undoubtedly explaining why the activity concentrations obtained with the two source terms are underestimated.
All of the available measurements have been taken into consideration in the validation process.The scores measuring the agreement between the observed and the modeled values are given Table 5.The proportions of the data where simulations and observations agreed within a factor of 2 are low especially for the 131 I but the proportions of the data where simulations and observations agreed within a factor of 5 and 10 are acceptable and similar for all the radionuclides.No distinction between gaseous and aerosol-bound iodine have been made in the simulations.The failure to take the complex behavior of iodine into consideration could be partly responsible for the lower performance with respect to iodine.The results obtained with the retrieved ST are slightly better than those obtained with the Mathieu et al. (2012) ST.
The comparisons between the modeled and observed atmospheric activity concentrations partially validate the isotopic composition of the releases measured in the South zone.They suggest that the released quantities of 137 Cs, 136 Cs and 131 I in the retrieved ST are realistic.The 132 I and 132 Te releases are underestimated.Indeed, 72 % of the data where the observed and modeled 132 Te values agreed within a factor of 5 are underestimated.The 132 Te underestimation is balanced by a probable overestimation of the 133 Xe released quantities.Tsuruta et al. (2013) presented a new set of activity concentration measurements in the Fukushima Prefecture (not public data) which should help to improve the validation of the isotopic composition of the retrieved ST.

Comparison with the surface activity measurements
The assessment of the release's isotopic composition is completed by comparing the simulation results with the surface activity measurements.Like the activity concentration measurements, these were not used to estimate the source terms; as a result, they provide an effective and independent means of validation.The daily deposition measurements by prefecture for Japan, taken by the MEXT (2011c), are effective only from 18 March.No data could be collected in the Fukushima prefecture.Table 6 gives the proportion of the simulated daily deposition that is within a factor of 2, 5 and 10 of the observed values.For the whole of Japan, the scores are not as good as for activity concentration (the proportion of the simulated deposition that is within a factor of 5 of the observed values is 40 % for the retrieved ST and only 30 % for the Mathieu et al. (2012) ST) but they are significantly higher when considering only the measurements in the South zone.
In the northern zone (Iwate and Yamagata), the deposition amounts increased on 20 March (event # 11) but the event is poorly reconstructed, as the depositions caused by event # 6 are overestimated (see Fig. 6).In the prefecture of Ibaraki (Hitachinaka) within the southern zone, the highest deposition amounts were on 20 and 21 March (event # 11), and on 21 and 22 March towards Utsonomiya, Tokyo, Saitama, Kofu and Chiba.Figure 8 shows the degree to which the simulations are consistent with the observations for 137 Cs and 131 I, for two examples of monitoring stations.The first activity peak is correctly reconstructed, whereas the second is underestimated.In the case of the correctly reconstructed peak, the simulated 137 Cs and 131 I quantities are consistent to within a factor of 2 of the observations, which confirms the comparisons made regarding the activity concentration data and tends to validate these two radionuclides' release quantities.
The only available means of validating the release episodes that have contaminated the other zones (West, North-West, and North) is the total surface activity map.This is the result of several airborne measurement campaigns (MEXT and DOE 2011; US DOE/NNSA, 2011).These measurements, some of which were taken after the accident, are especially useful for studying the deposition of long-lived radionuclides.Figure 9 compares the observed surface activity of 137 Cs with the simulations.
The simulations performed using the two source terms produce different results.The deposition area to the northwest is largely underestimated using the Mathieu et al. (2012) ST.It is better represented using the retrieved ST, but it is too extensive and too far north compared with the observations.This mismatch is consistent with the overestimated dose rate levels in the Yamagata and Yonezawa prefectures.
The north part of the most contaminated pattern is slightly under-estimated by approximately a factor of 2.
To the south (notably the Iwaki region), the simulation performed with the Mathieu et al. (2012) ST results in the total deposition amount being considerably overestimated.The retrieved ST is much more consistent with the observations.Other, less important, deposition areas can be seen, to the west, bordering the Tochigi prefecture, and to the south of the Fukushima prefecture, 100-150 km away from the nuclear facility, upstream of the mountainous areas.These depositions resulting from the releases on 15 March (event # 6) are incorrectly reproduced by the simulations.This explains why the dose rates are considerably underestimated in the western zone.The differences can be seen with the two source terms and are due to insufficiently precise meteorological data.As a result, they do not invalidate the retrieved ST's realism.
In conclusion, the simulations with the retrieved ST are, on average, consistent with the observations in the most heavily contaminated zones (the North-West and South zones).This result helps to validate the 137 Cs emission quantities in the retrieved ST.The same can be said of the 134 Cs comparisons.

Conclusions and prospects
In this paper, we have proposed an effective method of reconstructing the source term of a nuclear accident.The approach is based on inverse modeling techniques, and differs from the existing methods in that it is designed to use dose rate measurements.Dose rate stations are the most common device and, in the event of a nuclear accident, they are the primary source of plume and deposition measurements during releases.It is therefore crucial to be able to use dose rate data in an inverse modeling approach to assess source term.This is a difficult task, as the measurements include the contributions of all gamma emitters on the ground and in the air without distinguishing between them.
The proposed approach consists of three steps.The first step is dedicated to the determination of the a priori information required to solve the problem.The second uses inverse modeling techniques to identify the periods during which releases may have occurred.The third step consists of using inverse modeling to estimate the source term during the potential emission periods.
The method is generic, meaning that it can be applied to any nuclear accident when sufficient dose rate measurements are available.It can be applied with any atmospheric dispersion model.However, the nature of the dose rate measurements leaves too many degrees of freedom for the inverse problem, which is ill-conditioning as a result.The following prerequisites are therefore necessary in order to constrain the problem: -As more than a hundred different radionuclides may be emitted during an accident, those whose contribution to the dose rate signal is the greatest must be selected.
In the event of a nuclear power plant accident, this list can be a priori assumed.
-The isotopic composition of releases can be constrained by requiring that the ratio between the different radionuclides' release rates remain within defined limits.The defined limits may be more or less constraining, depending on the prior knowledge available regarding the releases.
The method has been validated on the Fukushima accident.
We have shown that it is appropriate for reconstructing the source term of a nuclear accident.Thanks to its flexibility (it is not dependent upon Fukushima) and its effectiveness (the calculation can be performed within a few hours), it can easily be used in crisis situations.However, it must be further improved in order to obtain a more reliable isotopic composition of the estimated releases.This is because the method is not always able to distinguish the respective contributions of the noble gases and those of the radionuclides with short half-lives.We plan to overcome this by extending the method to integrate the activity concentrations and surface activities into the inversion process, in addition to the dose rate measurements.Properly used, these measurements should further constrain the releases' isotopic composition.
The approach has been applied to the Fukushima accident, and it provides estimated release rates at hourly intervals between 11 March at 00:00 UTC and 26 March at 21:00 UTC, for all radionuclides which mainly contribute to the dose.There are eight of these radionuclides.They consist of 134 Cs, 136 Cs, 137 Cs, 137m Ba, 131 I, 132 I, 132 Te, and 133 Xe.We use 57 dose rate monitoring stations covering the Japanese territory.The method is applied by using the ldX atmospheric dispersion model and the ConsX dose calculation model from the IRSN C3X operational platform.It is estimated that the emissions amounted to 105.9 PBq for 131 I, 35.8 PBq for 132 I, 15.5 PBq for 137 Cs and 12134 PBq for the noble gases.
The approach succeeded in finding the different release periods correlated with the FD-NPP's events, without requiring prior knowledge of emissions.It has confirmed and validated the release rates emitted during certain events and already estimated in other studies.It has also helped to improve the reconstruction of several release episodes.
In order to validate the temporal evolution of the estimated release rates, the simulations of the atmospheric dispersion and of the retrieved ST's consequences have been compared with the environmental observations.On overall, this shows that the simulations are consistent with the dose rate measurements, since it was found that for about 80 % of the observations, simulated observed dose rates agreed within a factor of 2. The temporal evolution of the dose rate is generally well reconstructed, particularly in the most heavily contaminated zones to the north-west and to the south of the FD-NPP.The results are more mixed in some zones to the north and to the west of the FD-NPP, as the meteorological data do not take the complex orography of this region of Japan sufficiently into account.Concerning the isotopic composition of the releases, the comparison with the activity concentration and surface activity observations suggests that although the estimated caesium and 131 I emissions are realistic, they are underestimated for 132 I and 132 Te and probably overestimated for the noble gas 133 Xe.
We plan to further improve the accuracy of the FD-NPP release estimations by utilizing more detailed meteorological data in the method, and adopting a multi-scale approach that would make it possible to use measurements that can either be close to the source or further from it.This assumes that the measurements' relative importance is taken into account.At present, the absolute value of the model-to-measurement difference is considered in the cost function to be minimized.As a result, more weight is given to the most contaminated zones.By adapting the method to a multi-scale approach, it should be capable both of reconstructing the releases that were transported to the Pacific Ocean, and of utilizing the measurements taken close to the FD-NPP, where the dose rate gradients are very steep.

Fig. 1 .
Fig. 1.Idealized dose rate signal.(1) shows the detection of a plume.(2)shows the signal behavior when the plume moves further away.(3)shows the dose rate signal due to the deposition.
is the error covariance matrix related to the measurements and model.E is the mathematical expectation.The matrix B = E (σ − σ b ) (σ − σ b ) T represents the error covariance matrix of the first guess, σ b .Release rate values are positive or null.

Fig. 2 .
Fig. 2. Extraction of plume component from the dose rate signal in Tsukuriya.The station is located in the Ibaraki prefecture, 139 km south of the facility.Japan Standard Time (JST) is 9 h ahead of Coordinated Universal Time (UTC).

Fig. 3 .
Fig. 3. Location of the dose rate monitoring stations used for the inversion.The color varies depending on the start date of the signal concerned; blue: before 12 March, green: 12-15 March, yellow: 15-20 March, red: after 20 March.The white dashed lines show the demarcation of the 4 geographical areas.The shape of the symbol is used to indicate to which geographical area the monitoring stations belong.The diamonds are for the North zone, the circles for the North-West zone, the down triangles for the West zone and the hexagons are for the South zone.The purple triangle indicates the FD-NPP's position.

Fig. 5 .
Fig. 5. Comparisons of the dose rate observations (black dots) with the simulated dose rate computed with the retrieved ST (red) and with the Mathieu et al. (2012) ST (blue).Example of closely corresponding values in the four zones.

Fig. 6 .
Fig. 6.Comparisons of the dose rate observations (black dots) with the simulated dose rate computed with the retrieved ST (red) and with the Mathieu et al. (2012) ST (blue).Example of poorly corresponding values in the North zone (top) and in the West zone (bottom).

Fig. 7 .
Fig. 7. Measured activity concentrations (white blocks), and activity concentrations simulated with the retrieved ST (red circle) and the Mathieu et al. (2012) ST (blue circle) in Tokai, for 137 Cs (a), 136 Cs (b), 131 I (c), and 132 Te (d).Circles are centered on the respective time interval measurements.

Fig. 8 .
Fig. 8. Daily surface activities of 137 Cs (top) and 131 I (bottom), observed (white blocks) and simulated using the retrieved ST (red circle) and the Mathieu et al. (2012) ST (blue circle).Circles are centered on the respective time interval measurements.

Table 1 .
Evolution in the number of release instants and of the error diminution obtained after minimizing the cost function J 1 according to the value of λ 1 .The results obtained for λ 2 = 10 −12 (selected value) are in bold font.

Table 2 .
Total emissions (in PBq) for the source terms of 133 Xe, 131 I, 132 I, 137 Cs and 136 Cs.

Table 4 .
Mathieu et al. (2012)erved and simulated dose rates using the retrieved ST and theMathieu et al. (2012)ST.The statistics are the proportion of the measurements for which observed and modeled values agreed within a factor of 2 and the bias.They were produced for Japan as a whole and for four zones located to the north, the north-west, the west and the south of the FD-NPP.The results obtained for Japan as a whole are in bold font.

Table 6 .
Mathieu et al. (2012)erved and simulated daily depositions using the retrieved ST and theMathieu et al. (2012)ST.The statistics are the proportion of the measurements for which observed and modeled values agreed within a factor of 2, 5 and 10.