Towards the operational application of inverse modelling for the - PowerPoint PPT Presentation

Towards the operational application of inverse modelling for the source identification and plume forecast of an accidental release of radionuclides Victor Winiarek (victor.winiarek@cerea.enpc.fr) Julius Vira, Marc Bocquet, Mikhail Sofiev, Olivier Saunier Universit´ e Paris-Est CEREA, Ecole des Ponts Paris-Tech / EDF R&D Finnish Meteorological Institute IRSN, Institute of Radiation Protection and Nuclear Safety V. Winiarek HARMO13, Paris, France, 1-4 June 2010 1 / 19

Problem definition and objectives Context In an event of an accidental release of radionuclides, authorities need to know (if possible in advance) the impacted area. Numerical models are used to forecast the radioactive plume. The performance of these tools are mainly forced by the knowledge of the source field. Data assimilation methods (such as inverse modelling) have shown, at least at an academic level, good skills to help in this matter. Objectives To propose data assimilation methods to help forecasting radionuclides plume simple enough to be implemented in an operational context (Which assumptions ? Which simplifications ?) and to be understood by operators but still efficient. V. Winiarek HARMO13, Paris, France, 1-4 June 2010 2 / 19

Outline Sequential semi-automatic data assimilation system 1 Bayesian tests for the identification of the release site 2 V. Winiarek HARMO13, Paris, France, 1-4 June 2010 3 / 19

Sequential semi-automatic data assimilation system Outline Sequential semi-automatic data assimilation system 1 Bayesian tests for the identification of the release site 2 V. Winiarek HARMO13, Paris, France, 1-4 June 2010 4 / 19

Sequential semi-automatic data assimilation system Framework of the study France : 20 nuclear facilities - Monitoring network of 100 stations (Saunier et al . 2009) Finland : 6 sites - Monitoring network of 255 stations (“uljas” network) V. Winiarek HARMO13, Paris, France, 1-4 June 2010 5 / 19

Sequential semi-automatic data assimilation system Accidental dispersion synthetic experiment Source term Hypothetical fast core meltdown, without hull 3e+09 breach (nuclear power plant) Release rate (Bq/s) 2e+09 Dispersion of caesium 137. Intentional release 24 hours after the start of 1e+09 the accident → “double-peak” temporal profile. 0 12 24 36 48 60 72 Time (hours after the start of the release) Transport simulation and perturbed observations The networks are supposed to monitor the activity concentrations of 137 Cs (actually, most of them measure γ -dose). Transport simulated with Polair3D or SILAM → computation of synthetic observations µ synth . . Lognormal perturbations of synthetic observations : µ perturb . ∼ exp( N (0 , 0 . 5)) µ synth . (1) . i i V. Winiarek HARMO13, Paris, France, 1-4 June 2010 6 / 19

Sequential semi-automatic data assimilation system Method (1/2) : Inverse modelling Inverse modelling Source-receptor relationship: µ = H σ + ε Assumption : the location of the accident is known → Unknown : temporal profile of the source σ ( N imp emission rates ; N imp of the order of 10 2 ) Many more observations (though very noisy) than unknown → direct computation of the Jacobian matrix H ∈ R N obs × N imp (column by column) The estimated reconstructed source is given by : � − 1 H T R − 1 H H T R − 1 µ � σ = (2) Prior modelling of errors Gaussian uniform R = χ I d . Gaussian relative R = diag ( χ 1 , χ 2 ,..., χ d ), with √ χ i ≃ 0 . 5 µ i (+ threshold) V. Winiarek HARMO13, Paris, France, 1-4 June 2010 7 / 19

Sequential semi-automatic data assimilation system Method (2/2) : Data assimilation scheme Reconstruction of the source (analysis) Measurements in the interval [ t a − ∆ t a , t a ] are collected (those in [ t 0 , t a − ∆ t a ] were already available). This allows to build the measurement vectors up to t a : µ a . One prolongates the source-receptor matrix H at t a , by prolongating or computing all the elementary solutions from t a − ∆ t a to t a . Then one computes an estimate of the source term σ a . Forecast A forecast is performed from t a to t f , using the transport model. Forecast driven by the best estimation of the source up to t a ( σ a ), then by a guess of the source term from t a to t f (usually persistence hypothesis). V. Winiarek HARMO13, Paris, France, 1-4 June 2010 8 / 19

Sequential semi-automatic data assimilation system Results (1/2) : Statistical indicators Reconstruction of the source � i =1 ([ σ ] i − [ σ t ] i ) 2 N ∑ N 1 rmse = [ σ − < σ > ] i [ σ t − < σ t > ] i ρ = ∑ N i =1 �� �� � j =1 [ σ − < σ > ] 2 j =1 [ σ t − < σ t > ] 2 ∑ N ∑ N j j Forecast 1 0.8 ∑ min([ c ] h , [ c t ] h ) Figure of merit 0.6 h ∈ S figure of merit = ∑ max([ c ] h , [ c t ] h ) 0.4 h ∈ S Median figure of merit 0.2 Average figure of merit 0 0 12 24 36 48 60 Time (hours after the start of the release) V. Winiarek HARMO13, Paris, France, 1-4 June 2010 9 / 19

Sequential semi-automatic data assimilation system Results (2/2) : Plume forecasts 25.00 t = 1 h t = 3 h t = 9 h t = 24 h 50 N 10.00 � 2.00 45 N � 1.00 40 N � 0.50 50 N � 0.20 45 N � 0.10 0.01 40 N � 0 � 10 � E 0 � 10 � E 0 � 10 � E 0 � 10 � E The source term is quickly (after 1-2 hours of observations) well-estimated. The forecast of the radioactive plume is of good quality. But in an operational context, the average behaviour of the system is not sufficient → One must pay attention to the cases where the system fails. V. Winiarek HARMO13, Paris, France, 1-4 June 2010 10 / 19

Sequential semi-automatic data assimilation system Fail cases : causes and solutions Analysis of the fail cases Some power plants are located on the shores, near the frontiers or far from the monitoring networks → Some accidents are not well-observed → H T R − 1 H is severely ill-conditioned and the inversion step is not achieved. One alternative solution : Regularisation Use of a background term : for example a True source Retrieved source without regularisation (Polair3D and SILAM) 3e+09 Retrieved source with regularisation (Polair3D) Gaussian assumption for the source term dis- Retrieved source with regularisation (SILAM) tribution, with B the background error co- Release rate (Bq/s) 2e+09 variance matrix, leads to a new estimation of the source : 1e+09 σ = ( H T R − 1 H + B − 1 ) − 1 H T R − 1 µ (3) 0 0 12 24 36 48 60 72 84 96 Time (hours after the start of the release) But... But an important part of the information is still lost → An international network would be the best solution. V. Winiarek HARMO13, Paris, France, 1-4 June 2010 11 / 19

Bayesian tests for the identification of the release site Outline Sequential semi-automatic data assimilation system 1 Bayesian tests for the identification of the release site 2 V. Winiarek HARMO13, Paris, France, 1-4 June 2010 12 / 19

Bayesian tests for the identification of the release site Objectives and general principles If the release site is unknown ? The operational data assimilation system needs to be complemented with more sophisticated methods that do not assume that the release localisation is known (parametrical or non-parametrical methods, progressive reduction of candidates group) or with statistical tools which indicate the probability of a power plant to be responsible for the accident, knowing the measurements. Bayesian tests Such statistical tests are based on Bayesian inference theory � p ( µ ) = p ( σ ) p ( µ | σ ) d σ , (4) and differ from each other by the assumptions made on the source prior p ( σ ) V. Winiarek HARMO13, Paris, France, 1-4 June 2010 13 / 19

Bayesian tests for the identification of the release site Gaussian prior (1/2) Principles If p ( σ ) follows a Gaussian multivariate distribution, one obtains � − 1 µ ) H i BH T 2 µ T � p i ( µ ) = exp( − 1 i + R (5) | H i BH T 1 i + R | 2 p i ( µ ) represents the likelihood of the dataset µ provided the source prior statistics are Gaussian, and that the source is located at site i. H i being the Jacobian matrix of site i (a submatrix of H ). V. Winiarek HARMO13, Paris, France, 1-4 June 2010 14 / 19

Bayesian tests for the identification of the release site Gaussian prior (2/2) Results 100 (a) The (hypothetical) accident occurs in 80 Sosnovy Bor power plant Likelihood (%) 60 Sosnovy Bor Forsmark, Oskarhamn and Kalinin After 3 hours, the likelihood is about 75% Loviisa and Olkiluoto 40 for this power plant 20 After 5 hours, 100% 0 1 3 6 Time (in hours after the start of the release) Range of validity 100 Belleville Dampierre Fessenheim 80 This indicator strongly depends on B . Likelihood (%) 60 But in this case, there is a range of 40 validity of four orders of magnitude. 20 0 1 2 4 6 8 10 12 14 10 10 10 10 10 10 10 Root mean square background error V. Winiarek HARMO13, Paris, France, 1-4 June 2010 15 / 19