© SFRP, 2020

1 Introduction

Upon the occurrence of an atmospheric accidental release from a nuclear power plant, the emergency situation has to deal with many unknowns regarding the kinetics of the nuclear event, the potential radioactive atmospheric releases and their amounts over time, the meteorology which may vary in the next days. All these factors influence the key decisions to be taken in real time by the authorities in charge, dealing essentially with the radiological protection of the population that may require sheltering or evacuation. Subsequently, part of the population will be potentially exposed to low doses ionizing radiation (IR) for a long period of time (Bourguignon et al., 2017).

Feedback from nuclear accidents, e.g., Fukushima, has shown that the first decisions regarding the management of the population are critical in terms of psychological effects and decisive to gain or lose confidence of the public (Lochard et al., 2019). Therefore, the best scientific support is needed for helping authorities to make the best possible decisions.

In case of an accidental situation involving radioactive material, Technical Support Organizations (TSO) such as the Institute for Radiation Protection and Nuclear Safety (IRSN) have to provide in support of the public authorities a scientific estimate of the releases and of their consequences on human health and environment. To produce a rapid and reliable expertise, the Technical Emergency Centre (TEC) of IRSN relies on an organisation based on the use of expert methods with specific tools. The C3X platform (Tombette et al., 2014) is used to calculate the atmospheric transport of the released radionuclides at local and regional scales and to assess the potential consequences on the population and the environment. In the context of emergency, decision-making relies on comprehensive information about the situation and its possible evolution. In particular, the source term has to be assessed. Currently, the method to estimate the radionuclide release is a bottom-up approach that requires information on the reactor of the power plant. This approach is essential since it allows us to obtain a prognosis of the release and thus to anticipate the consequences and to recommend protective actions.

However, past accidental situations such as in Chernobyl or in Fukushima-Daichii have shown that it can be difficult to rapidly obtain information on the status of damaged reactors. Under such circumstances, the assessment of the source term may be very uncertain. To overcome this difficulty, one way is to assess the source term using inverse modeling technique which consists of coupling atmospheric transport model with environmental observations. Following the Fukushima accident, IRSN developed an inverse modeling method based on a variational approach. This method has been applied to the Fukushima accident using dose rate measurements (Saunier et al., 2013) and ¹³⁷Cs air concentration measurements (Winiarek et al., 2012, 2014; Saunier et al., 2016).

The challenge is now to make this research tool operational to ensure that it can be encapsulated in the C3X platform and usable in an emergency context. More precisely, the challenge is to define a methodology to guide the user in applying good practices depending on the situation. The previous applications were done to reconstruct the Fukushima source term after the end of the releases when a lot of measurements were available (Saunier et al., 2013). One of the issues of the operationalization is henceforth to use the inverse modeling tool in real-time when the first significant radioactive levels of measurements are reported. National exercises are set up to test the ability of the crisis organization to respond to a crisis situation, they represent an ideal and a realistic application case. Therefore, the inverse modeling method has been used during national exercises. The measurements taken into account in the inverse modeling method are ambient gamma dose rate. In the present paper, the implementation of the method based on the use of ambient gamma dose rate is described and its ability to reconstruct a release in real-time is assessed.

2 Inverse modeling methodology using dose rate measurements

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 dose rate data interpretation is complex since the signal provides no direct information about isotopic composition or the respective contributions of the plume and deposit.

Inverse modeling methods based on dose rate measurements are therefore more difficult to implement than those based on air concentration measurements. Among the existing methods, the way to integrate information on the source term composition is the main challenge. Zhang et al. (2017) have developed an Ensemble Kalman filter where information on isotopic ratios is included within error covariance matrix associated with the source term. However, solving the inverse problem can lead to unrealistic source terms due to ill-conditioned matrix. This difficulty is overcome by adding null artificial measurements but this makes the algorithm more expensive. Moreover, this approach requires to manually determine coefficients of the error covariance matrix associated to the source term, which makes it less relevant for operational use. Similarly, Tichý et al. (2018) have developed a Bayesian method where isotopic ratios are also incorporated into the error covariance matrix associated to the source term. This approach is successfully applied to reconstruct the source term on the synthetic test case where measurements are simulated by the dispersion model. However, to date, none of these approaches has been applied to real accidental releases. The method used in this paper is based on a variational approach which allows to easily take into account information related to the composition of the release. This method was applied and validated to reconstruct the Fukushima source term (Saunier et al., 2013). In this method, a first step consists of reducing the size of the inverse problem by only selecting the main radionuclides that contribute to the dose rate.

Indeed, in case of an accident from a nuclear power plant, the main part of the dose rate signal is usually due to few radionuclides N. Knowledge related to the core inventory of the damaged facility may provide information about isotopic composition of the release. IRSN accidental database in which source terms are pre-calculated can be also exploited to determine the main radionuclides. Moreover, the relevant radionuclides can be set by using other types of environmental measurements such as spectrometry from air samples, although this information ask for longer delays.

In our approach, the composition of the source term σ is therefore reduced to this list of the main radionuclides N. The contributions of the main radionuclides N have to 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. We use a variational approach to assess the source term σ consisting in the minimization of a cost function J(σ). The first term of J(σ) measures the differences between the model predictions Hσ and the real dose rate measurements μ: with: H is the Jacobian matrix where each column represents the dispersion model’s response to a unitary release emitted for one radionuclide whose release rate is to be estimated. The gamma dose rate calculation at each element of H matrix is assessed from the activity concentrations and surface activities simulated by atmospheric transport model. This is done with the C3X platform’s ConsX model. With a low number of dose rate observations, the minimization of the J(σ) may fail by providing many solutions. A second term is therefore added in the cost function to regularize the inverse problem. It measures the differences between a priori (or background) source term σ_b and the updated source term σ.

R = E[εε^T] is the error covariance matrix related to the measurements and model. The vector ε is the observation error aggregating instrumental and modeling errors and B = E[(σ−σ_b)(σ−σ_b)^T] is the background error covariance matrix. Simple parametrization for R and B matrixes are used. It is assumed that they are diagonal and the error variance is the same for all diagonal elements of each matrix (homoscedasticity property): The parameter determines the scale of the fluctuations in the source term and has to be determined during the inversion process. Therefore, the cost function J(σ) takes the form: (1)The third term takes into account information about the isotopic composition of the release. On a physical point of view, the various radionuclides are released simultaneously in proportions that depend on their physicochemical properties and the core inventory. Therefore, the release rates of radionuclide j in relation with radionuclide j + 1 are made to respect the following proportions: with a_j and b_j the limit values of the isotopic relationship. They can be estimated using the core inventory of the damaged reactor or using pre-calculated source terms from IRSN accidental database. The L-BFGS-B algorithm (Liu and Nocedal, 1989) is used to minimize the cost function (1). It is a limited-memory quasi-Newton algorithm for bound constrained optimization. The positivity constraint on the source vector is enforced and no upper bound is used.

3 Application: Cruas national exercise

The inverse method was applied during a national exercise considering an accidental scenario at nuclear power plant located in Cruas in the southeast of France. The scenario involves a failure of a steam generator tube and leads to a continuous release emitted in the environment between 08:00 and 16:00. The accident involved rather small releases since it does not require actions to protect the population. However, the radioactive materials released into the environment are significant enough to be detected by dose rate monitoring network. The source term consists of 39 radionuclides but only few radionuclides contribute significantly to the dose rate signal (¹³⁴Cs, ¹³⁶Cs, ¹³⁷Cs, ^137mBa, ¹³¹I, ¹³³I, ¹³⁵I and noble gas). By their very nature, noble gases are the first radionuclides emitted in the early hours of the accident. It means that the ratio between noble gas and other radionuclides is very high during the early stages of the accident and then gradually decreases. The total released amount is 1.3 × 10¹³ Bq including 8.7 × 10¹² Bq of noble gases (⁸⁸Kr, ¹³³Xe, ¹³⁵Xe), 1.8 × 10¹² Bq of cesium and 2.2 × 10¹² Bq of iodine. The maximum release rate occurred between 12:00 and 13:00. The temporal resolution of the source term varies between 5 minutes and 15 minutes.

Early in the morning, the wind on Cruas power plant came from the north then it gradually shifted toward the north–northeast in the afternoon. Around twenty dose rate stations are located in the vicinity of the Cruas power plant (Fig. 1).

During the exercise, only five stations reported an increase of the dose rate due to the wind mostly blowing in one direction. Among these five stations, three are located less than 1 km from the power plant and the other stations are located 10 km south–southwest of the power plant. Therefore, there is a fundamental issue in assessing the relevance of the inversion method in a situation where the number of the measurements is small.

Fig. 1 Monitoring dose rate stations located around the Cruas nuclear power plant.

3.1 Computation of synthetic measurements

Since the releases from the Cruas power plant are fictitious, the “true” dose rate measurements μ_t are calculated using atmospheric transport model. At a local scale (from several hundred meters up to several tens of kilometers), the Gaussian puff model pX is used (Soulhac and Didier, 2008; Korsakissok et al., 2013) to simulate the atmospheric dispersion of the radioactive plume. The pX model is part of the C3X operational platform. Although Gaussian models are not suitable to describe the atmospheric dispersion in the very short range of a facility (0–1 km) due to their assumptions (flat terrain, point source...), they show a good performance within a few kilometers from the source and are widely used in a crisis context due to their simple approach and low computational burden. In this study, the Doury (1976) standard deviation laws are used. In the C3X platform, two types of meteorological fields can be considered:

• homogeneous in space but time-dependent data based on observations;

• tridimensional meteorological and time-dependent data provided by Meteo France.

In this study, the meteorological measurements provided by the Cruas power plant are used to generate the “true” dose rate measurements. The meteorological measurements are available every 30 minutes, temporal resolution of the “true” measurements is 10 minutes and calculated using the following formula: where σ_t is the “true” source term and no perturbation is added to synthetic measurements μ_t.

3.2 Real-time source term reconstruction

The aim consists of estimating the source term of atmospheric emissions between 08:00 and 16:00. The temporal resolution of the source term is assumed to be 10 minutes in accordance with the temporal resolution of the measurements. The different steps of the source term reconstruction are described below. The main objective is to focus on the ability of the inversion method to reconstruct the “true” source term (ST-T) using a small number of measurements. We have therefore chosen to make simplifying assumptions about meteorological data, pX dispersion model parameters and source term composition.

3.2.1 H matrix calculation

The H matrix is computed column by column (Winiarek et al., 2011) using the forward pX model. The matrix is updated every 30 minutes taking into account the changes in meteorological conditions. The computational time related to this step is not prohibitive since it requires only several seconds. To avoid additional model errors, the meteorological data and the physical parameters of pX are similar to those used to generate synthetic measurements.

3.2.2 A priori information

We assume that we have no knowledge about the source term which results in σ_b = 0. It means that the source term reconstruction only depends on the measurements. If a release event is not detected by a dose rate station, it may not be assessed by inverse modeling (σ = 0). In the case of sparse measurements, it is expected that several release periods will be omitted. Another appropriate choice would be to consider a non-zero a priori source term based for instance on the analysis of the status of the reactor.

3.2.3 Isotopic composition assumptions

For this accident scenario, the main part of the dose rate signal is due to few radionuclides: ¹³⁴Cs, ¹³⁶Cs, ¹³⁷Cs, ^137mBa, ¹³¹I, ¹³³I, ¹³⁵I and noble gases. As in initial approach, it is assumed that the main contributors are known. By definition, noble gases do not deposit and they only contribute to the radioactive plume. In our inversion process, it is therefore not possible to distinguish two isotopes in the form of noble gases. For that reason, all of the noble gases emitted during the accident are grouped and are estimated as ⁸⁸Kr emissions. Moreover, secular equilibrium is supposed for ^137mBa and his parent ¹³⁷Cs: (2) The half-lives of ¹³⁷Cs, ¹³⁴Cs, and ¹³⁶Cs (respectively 30.17 years, 2.06 years and 13.2 days) are very different but are too long as compared to the duration of the release to have separate impact on the dose rate signal due to deposition. The same is true for ¹³¹I, ¹³³I and ¹³⁵I. Therefore, constant isotopic ratios can be assumed for , , , . We considered that the following constant isotopic ratios are known: (3)

By the end, the inverse problem to solve consists in the assessment of release rates of ¹³⁷Cs, ¹³¹I, ⁸⁸Kr which are assessed by minimizing the cost function (1). For instance, the number of unknown source parameters reconstructed by inverse modeling at 16:00 will be equal to 48 × 3 = 144 considering a source term of 10 min temporal resolution. The release rates of ¹³⁴Cs, ^137mBa, ¹³⁶Cs, ¹³³I, ¹³⁵I are then determined in a simple way using equation (2) and equation (3). An analysis of the ST-T indicates that the bounds of the isotopic ratios between ¹³⁷Cs, ¹³¹I, ⁸⁸Kr are the following: In our case, larger bounds are willingly chosen in order to evaluate the ability of the inversion method to converge towards an optimal solution:

4 Results and discussion

The time calculation required to converge to a minimum of the cost function J(σ) is only few minutes. Thanks to the flexibility of the inversion method, an update of the source term has been carried out every 30 minutes between 08:00 and 16:00.

4.1 Determination of λ parameter using L-curve method

When the number of observations used to reconstruct a release event is low, the inverse problem (1) may not be sufficiently constrained. This will result in source term magnitude being very sensitive to the value of λ parameter. In that case, strong values of λ lead to source term tempered too much whereas small values of λ yield to unrealistic source term in which the norm may tend to infinity. In practice, the choice of the λ parameter therefore involves a compromise between reducing the error value related to the cost function and increasing the amount of radionuclides emitted. In this study, the optimal choice of λ is determined using L-curve method (Hansen et al., 1989). It consists in plotting the source term norm versus the norm of the first term of J(σ) (residues) and identifying the maximum curve point in the graph. For every source term assessment (10:00, 13:00 and 16:00), an optimal value of λ is determined applying the L-curve technique. The initial λ values are in the range of 10⁻⁹ to 10⁻³. Figure 2 shows the evolution of the residues of J(σ) as a function of the ¹³⁷Cs source term assessed at 10:00, 13:00 and 16:00. Figure 2 emphasizes that the maximum curve point is reached for inversion performed with λ = 10⁻⁴ at 10:00 and for inversion performed with λ = 10⁻⁵ at 13:00 and 16:00. The results set out below are based on the optimal values of λ.

Fig. 2 Log-log plot of the norm of the 137Cs source term as a function of the norm of residues related to the cost function J. Inversion procedure is performed at 10:00 (blue), 13:00 (green) and 16:00 (red) with λ included in the range [10−9; 10−3].

4.2 Source term reconstruction

The ¹³⁷Cs, ¹³¹I and ⁸⁸Kr source terms reconstructed by inverse modeling are plotted in Figure 3 at three different hours (10:00, 13:00 and 16:00). The source term is validated by comparing the environmental dose rate measurements with the atmospheric dispersion simulation performed by forcing the pX model with the inverted source terms. The objective of the model-to-data comparisons is to check whether all of the release events are reconstructed properly and in the right time scale, and whether the quantities and the isotopic composition of the releases are consistent with the ST-T source term. The Figure 4 shows comparison of the observed and simulated dose rates for five monitoring stations: CS0KRS921MA, CS0KRS801MA, CS0KRS802MA, CS0KRS810MA and CS-PPI-22.

Fig. 3 137Cs, 131I and 88Kr release rates obtained by inverse modeling (left: ST-1, estimated at 10:00, middle: ST-2, estimated at 13:00 and right: ST-3, estimated at 16:00). The whole period from 08:00 to 16:00 is shown on all figures.

Fig. 4 Comparisons of the dose rate observations (black dots) with the simulated dose rate computed with reconstructed source term at different times (red: 10:00, green: 13:00, blue: 16:00).

4.2.1 Source term reconstruction at 10:00

As shown in Figure 3 left column, the comparison of the inverted source term calculated at 10:00 (ST-1) and the “true” source term (ST-T) highlights some discrepancies in terms of magnitude of the release rates. Moreover, the radionuclide composition of the ST-1 source term is not very realistic. Indeed, the release rates of noble gases assessed by inverse modeling are underestimated whereas the releases rates in ¹³⁷Cs and ¹³¹I are overestimated. Since only one station (the CS0KRS810MA station) detects a rise in the dose rate between 08:00 and 10:00, it is therefore not surprising that the inversion procedure has difficulties to properly reconstruct the radionuclide composition of the source term between 08:00 and 10:00. It is reasonable to suppose that the isotopic composition reconstruction would have been more reliable if more stations have detected the plume. Moreover, the method is not always able to distinguish between the respective contributions of the noble gases and those of the radionuclides with short half-lives when the time window of the source reconstruction is very short (two hours for ST-1). However, as shown in Figure 4, the observed dose rate is well enough reproduced at CS0KRS810MA station. In particular, the main peak occurring at 08:20 is well simulated using the ST-1 source term. It highlights the ability of the inversion method to reproduce observed dose rate measurements even if the radionuclide composition of the ST-1 source term can be improved.

4.2.2 Source term reconstruction at 13:00

The Figure 3 (central column) shows that the source term assessed at 13:00 (ST-2) is rather consistent with the ST-T. In particular, the magnitude, the duration of the main peak and the radionuclide composition are well reproduced. We notice that four dose rate stations detect the plume between 10:00 and 13:00 which explains a better relevance of the source term reconstruction. Moreover, the releases that occurred between 08:00 and 10:00 changed slightly in comparison with the previous estimation done at 10:00. The release rates of ¹³⁷Cs and ¹³¹I are lower whereas the release rates of noble gases are higher. Observations between 10:00 and 12:00 at CS0KRS810MA station (Fig. 4) enabled to better constrain the radionuclide composition of the reconstructed source term.

At 12:00, we observe a significant increase of the dose rate at CS0KRS921MA, CS0KRS801MA and PPI-22 stations (Fig. 4). In general, the peak is very well reproduced by the simulation. From 13:00, the simulated dose rate decreases at PPI-22, CS0KRS921MA and CS0KRS801MA stations but remains higher than before the detection of plume which is consistent with observations.

4.2.3 Source term reconstruction at 16:00

The last assessment of the source term (ST-3) is carried out at 16:00. Between 13:00 and 14:00, the releases of the ST-T are constant but lower than during the previous hour. The observed dose rates ranged from 0.1 μSv/h to 2.5 μSv/h during this period.

These releases are well reconstructed by inverse modeling since five dose rate stations detect the plume. After 14:00, the release rates of the ST-T become higher and remain constant between 14:30 and 16:00. Due to a slight change in the wind direction, the dose rate on CS0KRS802MA station increases from 14:00 whereas it starts to rise from 15:00 at CS0KRS921MA station. We observe that the inverse modeling method is able to reproduce every peak even if overestimation is noticed at both stations. Between 08:00 and 13:00, we also notice that the ST-3 is very similar compared to the previous assessments ST-1 and ST-2. On average, the simulations performed with the inverted ST-3 are in good agreement with the dose rate observations. For 95% of the measurements, observed and modeled values agreed within a factor of 2. The retrieved total emissions of ST-3 in cesium, iodine and noble gases are respectively 1.1 × 10¹² Bq, 1.2 × 10¹² Bq and 5.6 × 10¹¹. The quantities are therefore underestimated compared to the ST-T emissions. The comparison between ST-3 and ST-T for noble gas emissions, though, should be treated with caution since all of the noble gases emitted during the accident are grouped and are estimated as ⁸⁸Kr emissions equivalent in the inversion procedure. For instance, ¹³³Xe accounts for 80% of the total amount of noble gases included in ST-T source term but ¹³³Xe contributes very little to the dose rate contrary to ⁸⁸Kr which is one of the main contributors to the dose rate. It is therefore logical that the amount of noble gases obtained by inversion method be lower than that of ST-T source term. On the other hand, the iodine and cesium emissions assessed by inverse modeling are very consistent with those of ST-T.

4.3 Radionuclide composition assessment

Since only dose rate measurements were taken into account to reconstruct the source term, comparing the simulations with another variable such as ¹³⁷Cs total deposit is a relevant way of validating the ST-3 source term. The ST-3 source term was therefore used to simulate the ¹³⁷Cs total deposit at 16:00 using the pX model. Simple models of deposition were chosen by considering a constant deposition velocity of v_dep = 2.10⁻³ m.s⁻¹. For wet scavenging, the parameterization used is Λ_s = Λ₀p₀, where Λ₀ = 5 × 10⁻⁵ h/(mm.s) and p₀ the rainfall intensity in millimeters per hour. Figure 5 shows the ¹³⁷Cs total deposit simulated at 16:00 using the ST-T source term (reference simulation) and the ST-3 source term.

We first notice that the deposit simulated to the east of the Cruas power plant is higher when the ST-3 source term is used. Indeed, the area where deposits are above 100 Bq.m⁻² extends up to 10 km south-southwest of the power plant using ST-3 source term, whereas it does not exceed 4 km for the reference simulation. These deposits result from releases occurring between 08:00 and 10:00, period during which the ¹³⁷Cs quantities in the ST-3 source term are larger than those of the ST-T source term (Fig. 3). The other deposition areas reconstructed using the ST-3 source term are fully consistent with the reference simulation, in particular the zone where deposits values are higher than 10,000 Bq.m⁻².

Fig. 5 Total 137Cs deposition simulated at the end of the exercise (left) using ”true” source term (reference simulation) (right) using ST-3 source term.

5 Conclusion

In this study, a variational inverse method was used during the Cruas national exercise to assess radioactive releases in the environment. In particular, the L-curve technique was applied in order to efficiently determine the λ parameter ensuring that the inverse problem is well posed. We demonstrated that the inverse modeling method is efficient in reconstructing the main release events despite a small number of measurements. Other inverse modeling methods (Tichy et al., 2017; Zhang et al., 2017) are able to efficiently reconstruct reliable source term on similar synthetic cases but their relevance has not been assessed when a small number of measurements is available.

The comparison between “true” and simulated measurements shows a very satisfactory agreement. Moreover, the calculation time required to assess the inverted source term are fully compatible with a real-time use. For a first real-time use, we made simplifying assumptions since no perturbation is added to “true” measurements and model error coming from the meteorological fields and atmospheric transport model is zero. The only source of error explaining the differences between the inverted source term (ST-1, ST-2 and ST-3) and the ST-T comes from the assumptions about the isotopic composition of the release. However, we showed that the isotopic composition of the inverted source term could be improved. It is explained by the low number of measurements and by the very short time window of the source reconstruction. In the latter case, the inversion process is not able to distinguish two radionuclides from the slope of the dose rate signal. Recently, the relevance of the inversion method has also been evaluated during another exercises using realistic 3-D meteorological fields and by taking into account measurement and dispersion modeling errors.

References

All Figures

	Fig. 1 Monitoring dose rate stations located around the Cruas nuclear power plant.
In the text

	Fig. 2 Log-log plot of the norm of the 137Cs source term as a function of the norm of residues related to the cost function J. Inversion procedure is performed at 10:00 (blue), 13:00 (green) and 16:00 (red) with λ included in the range [10−9; 10−3].
In the text

	Fig. 3 137Cs, 131I and 88Kr release rates obtained by inverse modeling (left: ST-1, estimated at 10:00, middle: ST-2, estimated at 13:00 and right: ST-3, estimated at 16:00). The whole period from 08:00 to 16:00 is shown on all figures.
In the text

	Fig. 4 Comparisons of the dose rate observations (black dots) with the simulated dose rate computed with reconstructed source term at different times (red: 10:00, green: 13:00, blue: 16:00).
In the text

	Fig. 5 Total 137Cs deposition simulated at the end of the exercise (left) using ”true” source term (reference simulation) (right) using ST-3 source term.
In the text