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

SLIDE 1

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

SLIDE 2

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

SLIDE 3

Outline

1

Sequential semi-automatic data assimilation system

2

Bayesian tests for the identification of the release site

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 3 / 19

SLIDE 4

Sequential semi-automatic data assimilation system

Outline

1

Sequential semi-automatic data assimilation system

2

Bayesian tests for the identification of the release site

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 4 / 19

SLIDE 5

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

SLIDE 6

Sequential semi-automatic data assimilation system

Accidental dispersion synthetic experiment

Source term Hypothetical fast core meltdown, without hull breach (nuclear power plant) Dispersion of caesium 137. Intentional release 24 hours after the start of the accident → “double-peak” temporal profile.

12 24 36 48 60 72 Time (hours after the start of the release) 1e+09 2e+09 3e+09 Release rate (Bq/s)

Transport simulation and perturbed observations The networks are supposed to monitor the activity concentrations of 137Cs (actually, most of them measure γ-dose). Transport simulated with Polair3D or SILAM → computation of synthetic

bservations µsynth..

Lognormal perturbations of synthetic observations : µperturb.

i

∼ exp(N(0,0.5))µsynth.

i

. (1)

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 6 / 19

SLIDE 7

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

f the source σ (Nimp emission rates ; Nimp of the order of 102)

Many more observations (though very noisy) than unknown → direct computation of the Jacobian matrix H ∈ RNobs×Nimp (column by column) The estimated reconstructed source is given by : σ =

HTR−1H

−1 HTR−1µ (2) Prior modelling of errors Gaussian uniform R = χId. Gaussian relative R = diag(χ1,χ2,...,χd), with √χi ≃ 0.5µi (+ threshold)

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 7 / 19

SLIDE 8

Sequential semi-automatic data assimilation system

Method (2/2) : Data assimilation scheme

Reconstruction of the source (analysis) Measurements in the interval [ta −∆ta,ta] are collected (those in [t0,ta −∆ta] were already available). This allows to build the measurement vectors up to ta : µa. One prolongates the source-receptor matrix H at ta, by prolongating or computing all the elementary solutions from ta −∆ta to ta. Then one computes an estimate of the source term σ a. Forecast A forecast is performed from ta to tf, using the transport model. Forecast driven by the best estimation of the source up to ta (σ a), then by a guess of the source term from ta to tf (usually persistence hypothesis).

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 8 / 19

SLIDE 9

Sequential semi-automatic data assimilation system

Results (1/2) : Statistical indicators

Reconstruction of the source rmse =

1

N ∑N i=1 ([σ]i −[σ t]i)2

ρ = ∑N

i=1 [σ−<σ>]i [σt−<σt>]i

∑N

j=1[σ−<σ>]2 j

∑N

j=1[σt−<σt>]2 j

Forecast

figure of merit =

∑

h∈S

min([c]h,[ct]h)

∑

h∈S

max([c]h,[ct]h)

12 24 36 48 60

Time (hours after the start of the release)

0.2 0.4 0.6 0.8 1

Figure of merit

Median figure of merit Average figure of merit

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 9 / 19

SLIDE 10

Sequential semi-automatic data assimilation system

Results (2/2) : Plume forecasts

0.01 0.10 0.20 0.50 1.00 2.00 10.00 25.00

40

N

45

N

50

N

t = 1 h

40

N

45

N

50

N
10
E

t = 3 h

10
E

t = 9 h

10
E

t = 24 h

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

SLIDE 11

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 → HTR−1H is severely ill-conditioned and the inversion step is not achieved. One alternative solution : Regularisation Use of a background term : for example a Gaussian assumption for the source term dis- tribution, with B the background error co- variance matrix, leads to a new estimation

f the source :

σ = (HTR−1H+B−1)−1HTR−1µ (3)

12 24 36 48 60 72 84 96

Time (hours after the start of the release)

1e+09 2e+09 3e+09

Release rate (Bq/s)

True source Retrieved source without regularisation (Polair3D and SILAM) Retrieved source with regularisation (Polair3D) Retrieved source with regularisation (SILAM)

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

SLIDE 12

Bayesian tests for the identification of the release site

Outline

1

Sequential semi-automatic data assimilation system

2

Bayesian tests for the identification of the release site

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 12 / 19

SLIDE 13

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)

r 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

SLIDE 14

Bayesian tests for the identification of the release site

Gaussian prior (1/2)

Principles If p(σ) follows a Gaussian multivariate distribution, one obtains pi(µ) = exp(− 1

2µT

HiBHT

i +R

−1 µ) |HiBHT

i +R|

1 2

(5) pi(µ) represents the likelihood of the dataset µ provided the source prior statistics are Gaussian, and that the source is located at site i. Hi being the Jacobian matrix of site i (a submatrix of H).

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 14 / 19

SLIDE 15

Bayesian tests for the identification of the release site

Gaussian prior (2/2)

Results The (hypothetical) accident occurs in Sosnovy Bor power plant After 3 hours, the likelihood is about 75% for this power plant After 5 hours, 100%

1 3 6

Time (in hours after the start of the release)

20 40 60 80 100

Likelihood (%)

Sosnovy Bor Forsmark, Oskarhamn and Kalinin Loviisa and Olkiluoto

(a)

Range of validity This indicator strongly depends on B. But in this case, there is a range of validity of four orders of magnitude.

1 10

2

10

4

10

6

10

8

10

12

10

14

Root mean square background error

20 40 60 80 100

Likelihood (%)

Belleville Dampierre Fessenheim

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 15 / 19

SLIDE 16

Bayesian tests for the identification of the release site

Non-gaussian prior (1/2)

Principles Main assumption : The temporal profile σ b of the source rates is known, but not their real magnitude. Thus, the source is assumed to be of the form σ = λσ b Different prior p(λ) are assumed (gamma distribution, semi-gaussian) leading to different estimations of pi(µ). Example : with a semi-gaussian prior p(λ) =

2

π e− λ2

2θ

pi(µ) = e

(µTR−1Hi σb)2 2(θ−1+σT bHT i R−1Hi σb)

θ −1 +σT

bHT i R−1Hiσ b

×  1+Φ   µTR−1Hiσ b

2(θ −1 +σT

bHT i R−1Hiσ b)

    (6) where Φ(u) =

2 √π

u

0 e−x2dx is the error function.

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 16 / 19

SLIDE 17

Bayesian tests for the identification of the release site

Non-Gaussian prior (2/2)

Global performance When σ b is close from the true source, average results are excellent (2-3 hours to obtain at least 80% for the correct site). When σ b is different from the true source (wrong shape, wrong magnitude), the average performance still reaches 80%, but it does so later (typically 5 hours later). Critical example Fictitious accident in Kalinin power plant, with a poorly observed plume. The Gaussian test failed. The reconstruction of the source, knowing the release site was difficult (need of regularisation). The non-Gaussian tests managed to gradually identify the responsible site.

1 3 6 9 12 15 18

Time (in hours after the start of the release)

20 40 60 80 100

Likelihood (%)

Kalinin Forsmark Oskarhamn Sosnovy Bor Loviisa and Olkiluoto

(c)

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 17 / 19

SLIDE 18

Conclusion

A semi-automatic data assimilation system, using inverse modelling techniques, has been proposed to forecast a radioactive plume in an event of an accidental release from a nuclear power plant. Very good average performances (source quickly well-estimated, plume accurately forecasted). Fail situations have been identified and some complementary solutions have been proposed (regularisation, international network). In the case where the release site is unknown, statistical tests have been proposed and implemented to help identifying the responsible site. Excellent results in average and even good results in critical situations.

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 18 / 19

SLIDE 19

Selection of CEREA’s publications on the topic

Reconstruction of an atmospheric tracer source using the principle of maximum entropy. II: Applications. Bocquet M., Q. J. R. Meteorol. Soc., 131, part B (610), p. 2209-2223, 2005. Source reconstruction of an accidental radionuclide release at European scale. Krysta M. and M. Bocquet, Q. J. R. Meteorol. Soc., 133, 529-544, 2007. High resolution reconstruction of a tracer dispersion event: application to ETEX. Bocquet M., Q. J. R. Meteorol. Soc., 133, 1013-1026, 2007. Inverse modelling-based reconstruction of the Chernobyl source term available for long-range transport. Davoine X. and M. Bocquet, Atmos. Chem. Phys., 7, 1549-1564, 2007. Probing ETEX-II data set with inverse modelling. Krysta M., M. Bocquet, and J. Brandt, Atmos. Chem. Phys., 8, 3963-3971, 2008. Model reduction via principal component truncation for the optimal design of atmospheric monitoring networks. Saunier O., M. Bocquet, A. Mathieu and O. Isnard, Atmos. Env., 43, 4940-4950, 2009. Targeting of observations for accidental atmospheric release monitoring. Abida R. and M. Bocquet. Atmos. Env., 43, 6312-6327, 2009. Towards the operational estimation of a radiological plume using data assimilation after an accidental atmospherical release. Winiarek V., J. Vira, M. Bocquet, M. Sofiev and O. Saunier. J. Env. Radioactivity (submitted).

V. Winiarek

HARMO13, Paris, France, 1-4 June 2010 19 / 19