Probing ETEX-II data set with inverse modelling

We give here an account on the results of source inversion of the ETEX-II experiment. Inversion has been performed with the maximum entropy method on the basis of non-zero measurements and in conjunction with a transport model P OLAIR 3D. The discrepancy scaling factor between the true and the reconstructed mass has been es- 5 timated to be equal to 7. The results contrast with the method’s performance on the ETEX-I source. In the latter case its mass has been reconstructed with an accuracy exceeding 80%. The large value of the discrepancy factor for ETEX-II could be ascribed to modelling di ﬃ culties, possibly linked not to the transport model itself but rather to the quality of the measurements.


Introduction
The European Tracer EXperiment (ETEX) was designed to evaluate forecast performance of dispersion models as well as readiness and communication capacity of European institutions in case of an accidental release of harmful substances to the atmosphere. The experiment consisted of two controlled releases of harmless and non-15 depositing inert tracers (perfluorocarbon compounds). Since the release term was controlled and thereby well-known, and since the deposition and chemistry processes were non-existing, the design of the experiment made it an excellent test on models of atmospheric transport and diffusion from point sources at continental scale.
The location of the two releases was in Monterfil, Brittany, France. The first release 20 (ETEX-I) took place on 23 October 1994 and the second release (ETEX-II) some three weeks later, on 14 November 1994. Both experiments were carried out under comparable meteorological conditions with a low pressure system situated north of Scotland and westerly wind over the release site and over Western and Central Europe. In this paper we focus on the second ETEX release which took place between 15:00 UTC Introduction Interactive Discussion romethylcyclopentane) were released with a constant rate. In the first stage after the release a narrow cloud was advected rapidly eastwards (see the modelled plume in Fig. 1). Advection was driven by strong winds preceding a frontal passage which was linked with the low pressure area moving from the North Atlantic towards Scandinavia and Russia. The wind during ETEX-II was, however, much stronger compared to the 5 ETEX-I episode, and the low pressure system over Western Europe was accompanied by heavy rain. The experimental campaign coincided with a modelling exercise in which 28 models took part. Modelling turned out to be more difficult than in the case of the first release (Brandt, 1998). The first ETEX experiment has been treated in many papers in which model results 10 have been compared to observations, see (e.g., Brandt et al., 1998). The ultimate evaluation of the outcome of model intercomparison was presented in the final reports in 1998 (Mosca et al., 1998;Graziani et al., 1998a) where results from 49 models were included. In (Mosca et al., 1998) it was concluded that: "almost all the models show a satisfactory agreement with the measured values".

15
Even though the meteorological conditions were quite similar, and the same measurement sites and procedures were used in both experiments, there was a significant difference in models' performance between ETEX-I and ETEX-II. In the latter case an overall bias of a model validated against the measurements is suggested to be of a factor of 10 (Brandt, 1998). The experiment's website (http://rem.jrc.cec.eu.int/etex/) 20 reports spatial and temporal analyses for a very limited number of time instants and monitoring stations. Figures of merit and time of arrival for some subsets of models are the only statistical indices disclosed by the JRC. In particular, a gap in the quality of model results for the two releases is underlined there. For the figure of merit in space (FMS) at 24 hours after the release, half of the models exceeded 20% but only one of 25 them reached the value of 40%. In the case of the first ETEX release, an FMS above 40% was attained by 4 models out of 28, and there were still models exceeding this value at 36 h and even 60 h after the release.
According to the detailed conclusions from the ETEX symposium in 1997 in Vienna, the following was stated about the ETEX-II release: "The second release showed larger discrepancies between observations and model results where all models significantly overpredicted surface concentrations after 12 h from the start of release. No clear explanation has yet been given for this result" (Nodop, 1997). The eventual evaluation of model intercomparison within ETEX-II was treated in the final report in 1998 (Graziani 5 et al., 1998b) where results from the 28 models were presented. Some of the conclusions in this report stated: "Compared to the first release, much less tracer was found". It was also reported that: "Not only fewer non-zero values were observed, but also the concentration levels were much lower". Furthermore, according to the conclusions: "the cloud is apparently broken up into several parts moving at different speed in different directions. At many sites, the tracer was detected intermittently, rather than in a continuous sequence, as it might be expected". In the end, the report settles that: "The generally poor performance of the long-range dispersion models employed in the realtime simulation of the release are to put in relation with the complexity of the weather situation after the end of the second ETEX experiment".

15
Bad models' performance was considered rather strange among the modelling community because all the models performed satisfactorily with respect to ETEX-I. Since all the models showed a large overestimation with respect to ETEX-II, the reason for this overprediction could be explained by one of the following two hypotheses: 1) there is a systematic error in all the model results e.g. due to the complexity of the weather 20 situation which could result in a strong vertical mixing, or 2) there is a systematic error in the measurements which results in a significant negative bias (i.e. all the measurements are too small). The first hypothesis was concluded as a fact in (Graziani et al., 1998b). In this study, an inverse modelling procedure that was evaluated against ETEX-I with very good results (Bocquet, 2007a)  We tackle in this work the reconstruction of the source of ETEX-II in order to check to what extent the retrieved temporal profile and the emitted mass agree with the true release. Source retrieval problem is addressed with the principle of maximum entropy on the mean. The method has been presented in (Bocquet, 2005a,c), and successfully 5 applied to inversion of the source of ETEX-I (Bocquet, 2007a) as well as the temporal and vertical profile of the source of the Chernobyl accident (Davoine and Bocquet, 2007).

Measurements
There are 2248 valid measurements collected in the experiment. In contrast, however, 10 with 969 non-zero values gathered in ETEX-I, only 476 measurements (12% of all the measurements) resulting from the second release are non-null. Solely the latter subset has been used to perform inversion. The reason is our previous experience with the inversion method which made it clear that adding additional zero measurements is inefficient (Bocquet, 2007a).

Inversion method
The retrieval of a source from measurements sums up to an inversion of a sourcemeasurement relationship. For a discrete formulation of the problem the relationship is of the form: Interactive Discussion want to invert and a comparatively small number of measurements d , H has been computed in an adjoint mode.

Duality test
The modelling of the transport phenomenon for the needs of the inversion procedure is ensured by our Eulerian chemistry transport model POLAIR3D. It is equipped with a 5 module devoted to radionuclides which has recently been validated on ETEX-I, Chernobyl and Algeciras releases (Quélo et al., 2007). There are two ways of obtaining an adjoint associated with a discrete direct problem. One can discretise an adjoint of an original continuous model or adjoin an already discretised one. We have chosen here the first option which has the advantage to be simple. It boils down to time reversal of meteorological fields in a direct model. The latter possibility could have also been chosen since POLAIR3D has an extensively tested and employed adjoint module. In that case an agreement between the synthetic measurements obtained via the adjoint model, c ( Whenever both of them are below 0.01 ng m 3 , it has been assumed that they match Interactive Discussion perfectly (negligible concentration is assessed in both cases). The mismatch indicator computed with respect to all the measurements, D=2248, is greater then 20% for 366 measurements (that is to say for 16.3% of the concentrations), whereas 112 measurements (5%) have a mismatch exceeding 50%. We deal in this paper, however, with an inversion of real measurements which are 5 subject to measurement and modelling error, Eq. (1). In comparison to them, the error linked to the way the adjoint solutions have been computed can be regarded as negligible.

Cost function
The inversion method allows to take into account some prior knowledge which is ap-10 propriate to an accidental source. This capacity is of crucial importance due to scarcity of information on an accidental release and hence the necessity of making use of any available pieces of it. Prior information to be exploited refers in particular to positivity and boundedness of the source. Mathematical form of such information is encoded in a cost function arising from, for instance, Bernoulli law. However, any other positive prior 15 law, such as an exponential law, would have served a decent purpose. The minimum of a cost function allows to shift from prior (ν(σ) for the source and ζ (ε) for the errors) to posterior probability density function (p(σ) and q(ε), respectively), and consequently defines the inverted source. Usually such a cost function is defined in the space of the discretised source σ ∈ R N (and errors ε ∈ R d ): Interactive Discussion maximum entropy on the mean approach. Ultimately, one deals with the following cost function (Bocquet, 2005b;Krysta and Bocquet, 2007): where β ∈ R d is the vector conjugated to the measurement vector µ. R is the observation error covariance matrix which parametrises the Gaussian prior law on the 5 errors, ζ (ε). It is chosen diagonal, of the simplest form R=χ I, with I the identity matrix in observation space. The estimated source and errors are then given, in terms of the minimum β of L, by the estimators: Note that m k is a mass scale that is taken constant on all grid-cells, m k ≡m. γ k is the 10 probability that a release occurs in cell k, which is also taken constant, γ k ≡γ, and set to a very low value.

L-curve estimates
The prior probability density functions, ν(σ) and ζ (ε), are defined with the help of some parameters whose values are set arbitrarily according to experience, rather than to 15 some objective criteria. In particular, for Bernoulli law, the emitted mass m is an example of such a parameter. Similarly, Gaussian distribution describing errors is parametrised with an arbitrary value of the standard deviation √ χ . Therefore, in the first place, we establish the optimal values of those parameters. Due to nonlinearity of the dependence of the cost function on m, the analysis has been undertaken with the L-Introduction  (Hansen, 1992), and its application to source inversion is detailed in (Davoine and Bocquet, 2007 Then, having set m=10 (dropping the conventional unit M/λ 2 ), we have applied the

Inversion results for ETEX-II
We present in this section the inversion results obtained for the ETEX-II experiment. Interactive Discussion to 87 kg. On the other hand, inversion on the 1.125 • ×1.125 • grid fails but, as shown in (Bocquet, 2005b), such a failure is bound to occur for a fine grid.

Marginal gain of information provided with by each observation
The assimilation of the data through the source inversion process results in a net gain of information accumulated by the optimality system. What the observations convey to 5 the optimality system is measured by a relative entropy number denoted K σ , Eq. (3). In a general, non-Gaussian case, a single measurement contribution to the gain K σ depends on the solution and hence on the actual values of the other measurements. That is why one investigates the marginal contribution ∂ µ i K σ of an observation µ i , for all the observations. In (Bocquet, 2007b), it was shown how to compute such a 10 marginal contribution. These contributions have also been computed for ETEX-II and are plotted in Fig. 7. Moreover, the Lagrange multipliers β have been pictured. They represent the marginal total gain of information (which can either serve to identify the source or be lost in the deblurring of the errors in the data): ∂ µ i K σ ,ε . They are also a direct measure of how 15 sensitive the optimality system is to each measurement seen as a constraint. Finally, because the prior hypothesis on the errors is Gaussian, they are also directly proportional to the diagnosed errors (ε=Rβ).
From Fig. 7, one learns that the diagnosed errors depend strongly on most measurements in France and Benelux: a change in a measurement results in an increase 20 of the errors, rather than a new piece of information on the source. This contrasts with the analogue result for ETEX-I (Bocquet, 2007b). Nor is it the case for ETEX-II data coming from Central and Eastern Europe. The information that is gained on the source when perturbing the latter measurements is rather strong (again as compared to ETEX-I). This can be explained by the fact that the plume moved eastwards rapidly Interactive Discussion mation got from the Western measurements (at least relatively to the similar ETEX-I analysis). So that, comparing to ETEX-I, the optimality system trusts more the Eastern observations, which were probably found more consistent than the Western ones.

Comparison to ETEX-I
We report here inversion results as have been obtained for the ETEX-I experiment 5 (Bocquet, 2007a), but improved by parameters' estimation. In the first ETEX release 340 kg of PMCH were emitted to the atmosphere. Before source inversion has been addressed, a similar L-curve analysis has been performed. Consequently, an en- In order to enable comparison with the ETEX-II retrievals we show in Fig. 9 a reconstructed profile corresponding to the one given in Fig. 6 and the spatial distribution of the reconstructed mass in Fig. 8, an analogue of Fig. 5.

Remarks
The localisation of the ETEX-II source is rather good. 71% of the recovered mass is attributed to the correct grid-cell, as compared to 65% in the ETEX-I case. This can be explained by the optimality system to deblur some of the data. It is also explained by very high sensitivity of the source to the nearby measurements which makes the We note that the inverse modelling analysis of ETEX-II signals weak consistency in the data. It can be observed in Fig. 10, where the measurements arising from the reconstructed source are largely scattered as a function of the real measurements. Indeed, the inconsistency could have already been remarked in data -simulation comparison. The inverse modelling analysis, however, goes one step further because it is 5 capable of accounting for errors. Since model is taken into account in this study, and the errors are retrieved by the reconstruction method, the detected inconsistencies go beyond statistical ones which could be diagnosed through geostatistical procedures (Dubois et al., 2005). Still, consistency cannot be improved.

10
We have recalled here the results of source retrieval in the case of the ETEX-I release. The reconstruction performed with the method of maximum entropy has been shown to render satisfactory results in this case. The same procedure employed in source retrieval of the ETEX-II release gives much worse results, although an optimal framework of inversion has been much searched for.

15
Similarity of both experimental set-ups let us suppose that it is not the inversion method which is to blame. It does perform well in the former case. One could, however, point several other reasons for the difficulties encountered by the inversion procedure in the latter case.
Firstly, the number of the measurements is limited which has an impact on the attain-20 able source resolution (coarser than for ETEX-I). Secondly, the quality of the inversion is highly dependent on the possibility of incorporating error into the modelling procedure. Because convection is not properly taken into account, merely fronts, a sudden splitting of the plume may have been missed, resulting in an important model error and a loss of mass not reproduced in the modelling. Introduction

Back Close
Full Screen / Esc

Printer-friendly Version
Interactive Discussion from the measurements. This could be either due to poor modelling or poor measurements, or both. Our study rises that level up to 85%, as if the observed concentrations had been diluted by an overall mean factor of 6.8. From the reconstruction point of view the data seem inconsistent, in particular the observations performed in the early evolution of the plume (over France and Benelux).

5
Other clues, independent from this study, plead against the data. Although PMCP concentrations extracted from the collected samples are of good quality, the samples themselves are more difficult to judge. In many stations cloud's presence was intermittent. In particular, several cases are reported in which some stations exhibit cloud passage twice during the sampling period while another station, located close by, misses the second peak. Hence, the measurements suggest plume was chopped into several clouds which moved separately with different velocities (Nodop, 1998). This strange non-physical behaviour is, however, not supported by any of the models involved in model intercomparison (Graziani et al., 1998b). Moreover, the aircraft data (Girardi et al., 1998) do not reveal any splitting of the plume. In addition, even though it has 15 been argued against any water dilution of the samples, it has also been mentioned that the instruments' catalyst (Palladium) have been repeatedly poisoned by an unknown compound, so that the catalyst needed replacement several times. Introduction    Fig. 10. Posterior measurement consistency test. Measurements (diagnosed PMCP) simulated from the reconstructed source (Fig. 6) as a function of the real measurements (observed PMCP). The dashed lines mark the discrepancy factor equal to 2.