Petra Seibert1
Institute of Meteorology and Physics
University of Agricultural Sciences
Türkenschanzstr. 18, A-1180 Vienna, Austria
seibert_+aT+_.ac.at,
http://www.boku.ac.at/imp/envmet/
Andreas Stohl
Lehrstuhl für Bioklimatologie und Immissionsforschung,
Technical University of Munich, Freising, Germany
http://www.fw.tum.de/LST/METEOR/stohl/astohl.html
| The derivation of source-receptor matrices for a passive, nonreacting and nondepositing tracer from backward runs of a Lagrangian particle dispersion model is shown. This technique is applied to the first release of the European Tracer Experiment ETEX. The source-receptor matrix is used to invert measurements to obtain the temporal evolution of the release with given release location as a first illustration and test of the method. |
The purpose of this study was to test the potential of inverse modelling with data from the European Tracer Experiment ETEX [2,]. In this context, inverse modelling means to derive information on a source of an atmospheric trace constituent from ambient concentration measurements, using a dispersion model for the description of the source-receptor relationship. At first, we show how a Lagrangian particle dispersion model (LPDM) can be used in a backward mode to derive source-receptor matrices. This is a novel approach. Then the technique is applied to the ETEX release and a regularised inversion is performed to deduct the temporal evolution of the release. This is, however, just a very first application, and the outlook outlines possible refinements and further evaluations.
There are a number of different approaches to inverse modelling of atmospheric trace constituents. For a linear problem (ambient concentrations proportional to the source strength), one can obtain the unknown source from the solution of a linear system of equations, with measurements as knowns and elements of the source-receptor matrix (SRM) as coefficients. Depending on the number of knowns and unknowns as well as the conditioning of the matrix, one often needs to apply a so-called regularisation (input of a-priori knowledge) [6]. The SRM can be obtained by forward or backward modelling, and with a Lagrangian or a Eulerian model. The backward approach is most efficient when the number of unknown source elements exceeds the number of measurements. Lagrangian models have the advantage compared to Eulerian models that they have a better numerical accuracy (mainly, no numerical diffusion due to the advection scheme), and that they do not suffer from a grid resolution problem near point sources. In the backward mode, this is especially important because measurements are almost always made at a point or along a line (aircraft) and with Eulerian models one would need to apply a near-source sub-model or else loose once more in accuracy.
The other approach is to search iteratively for the minimum of a cost function, made up by the misfit between measurements and model, and a regularisation term if necessary. This requires the calculation of the gradient of this cost function with respect to the source elements. Often, this approach is realised with a Eulerian model and its adjoint, the only viable approach for large nonlinear problems [1]. It is also possible to use a Eulerian model and its adjoint for the determination of a SRM by defining suitable cost functions.
In the forward mode, particles are released from prescribed source areas, and transported with the mean and turbulent velocity field. They carry a mass depending on source strength and particle release rate. It can be altered by processes such as deposition or decay. Concentrations are calculated by summing up the masses of all the particles in a grid cell, and dividing the sum by the cell's volume.
In the backward mode, used for inverse modelling, the same formalism and computer model are applied, but the particle trajectories are integrated backward in time, using a negative time step. However, there is a new interpretation: particle masses are replaced by mixing ratios c º c/r (where c is the mass concentration and r the air density) in an infinitesimal volume of air. In our Lagrangian frame of reference, c can change only by sinks such as deposition or decay and by uptake of emissions. Without deposition, the individual rate of change of the mixing ratio is
| (1) |
The sources are independent of c, and sinks not considered as we are dealing with an inert tracer2. Eq. 1 can be easily integrated, yielding the instantaneous concentration. In ETEX (and many other practical applications), measured concentrations represent a point in space and an average in time (ETEX: sampling interval 3 h). Thus, we need to average in time. If this averaging time exceeds substantially the time scale of the turbulent fluctuations, the temporal averaging alone is sufficient and additional ensemble averaging is not necessary. Next, we introduce a discretisation in time (index j) with J back trajectories in equal intervals between t1 and t2, and index i for the temporal and spatial variation of S+ (``gridded in space and time''):
| (2) |
Let us now consider that we have not just one measurement, but many (different sites and times), denoted by l. Then we see that the elements of the desired source-receptor matrix M are
| (3) |
| (4) |
If the density variations between the receptors and the sources are small, mixing ratios c may be substituted by mass concentrations, and S+ by the emitted mass per volume and time, the parameter used normally. Otherwise, c and S+ can be transformed with the air density (climatological values are sufficient).
The Lagrangian particle dispersion model FLEXPART (details see section 3.2) was used in this application. FLEXPART gives gridded concentration fields as output. How can we obtain gridded residence times from this output? The time-averaged concentration in a grid cell [`c]i as computed in FLEXPART can be written as
| (5) |
Using the residence fraction f, the average residence time Dti during DT can be calculated as
| (6) |
| (7) |
The European Tracer Experiment ETEX comprised two releases of a gaseous, nondepositing, inert tracer. We are dealing here only with the first one. This experiment can be characterised as follows [2,]:
In the context of backward modelling, confusion can easily arise because, for example, the physical receptors are used as sources in the inverse model runs, whereas the physical source location acts as a receptor location. Therefore, we introduce the following notation: Source, concentration, receptor, etc. denote the physical system or a forward simulation; the same words with a * refer to the backward simulation. Thus, for example, a receptor is used as a source*.
FLEXPART is a LPDM developed by Andreas Stohl, and tested on various tracer experiments [8]. It performed well for the ETEX-1 release (correlation coefficient between observations and simulated concentrations 0.59, fractional bias 0.01).
FLEXPART was run with the following specifications:
Thus, the size of our source-receptor matrix is (41×21×9×35) × (168×30) = 180,810 × 5,740 » 1.4 ·109 elements. In sparse matrix format, this amounts to » 300 MByte of data. This matrix is too big for direct evaluation. Therefore, we have to consider the following options for simplifications:
1) Take only the lowest layer (or a vertical integral) as potential source, as concentrations* were well mixed near the release site.
2) Assume either the location as given (sample application: nuclear power plant accident), or release time (sample application: nuclear explosion with explosion time known from nuclide ratios).
3) Start with a coarse source resolution and narrow it down iteratively while reducing the domain size, horizontally, vertically, and/or temporally.
Here, we have applied simplifications 1 and 2 (release location given).
For a first visualisation, we plot fields that would result if a backward simulation would have been made with releases* proportional to the measured concentrations at each receptor. The plot may show horizontal fields for one time interval (or a loop over the whole simulation), a time-height-section at the release site, or just the time evolution at the release site. All these evaluations can be done for each monitoring site separately, or as the sum of all of them. This corresponds to classical ``trajectory statistics'' [7], or the ``adjoint tracer'' fields presented by Pudykiewicz [4]. This is not yet an inversion, but it is nice to illustrate the backward simulation.
Figure 1: Time series of a tracer* released at the measurements sites
proportional to measured tracer concentrations, sampled at the true
release site. Each line corresponds to one measurement location
(denoted by its ETEX station code) and the sum over all measurement
intervals. Asterisks denote the maxima, and vertical
lines the begin and end of the release. Note the logarithmic
scale.
Figure 1 shows that this method gives already a good picture of the timing of the release. The maximum lies within the real release interval for most stations, but some of them (e.g., D30 = Mannheim in Germany) are clear outliers. Figure 2 gives an idea how the ``adjoint tracer'' cloud* looks like at different times. The left frame is at the time of the release and shows its maximum near the release location while the right frame is at a later time (previos stage of the backward simulation) and looks similar to the ``real'' cloud at this stage as depicted by measurements. This figure may be compared with Fig. 3 in [5].
Figure 2: Horizontal distribution of a tracer* released at the measurements sites proportional to measured tracer concentrations. Note the logarithmic scale. The real release location is marked by a hollow black dot. The left (top) panel shows the tracer during the time of the real release (23 Oct, 18-21 UTC), the right (bottom) panel at a previous stage of the inverse simulation (25 Oct, 21-24 UTC).
We perform a generalised inversion (this means a regularised solution, see [6]) of the LSE yo » Mx (Eq. 4), where yo are the observations, minimising the cost function
| (8) |
Figure 3 shows the result. We see that the start of the release was accurately captured. The end of the release was found well but is too fuzzy and slightly biased towards later times. As the peak caused by station Mannheim at 50 h shows, the inversion is disturbed to some extent by outliers.
Figure 3: ETEX Release 1, true time-dependent release rate (merged into
3 h intervals) and release rate obtained by inverse modelling (with and
without the outlier station Mannheim).
Observed concentrations and concentrations computed by following methods have been compared as a test.
The results are presented in Table 1. We see that the SRM of the inverse run yields worse results than the forward run, but results with the reconstructed source term are better than with the real source term with respect to correlation and RMS error (not, of course, for bias). The total release is underestimated by 35% inspite of an unbiased forward simulation. We can thus conclude that the backward simulation suffers from the 3-hourly temporal and the grid-cell horizontal resolution of the source. The bias is obviously caused by formulation of cost function in combination with simulation errors.
Table 1: Statistical performance indicators for the
comparison modelled - measured tracer concentrations. Meaning of
abbreviations for the method see text.
| Method | r | relative bias | relative RMSE | |x| |
| FWD | 0.59 | » 0.01 | ?(not published) | 340 kg |
| BWD-i1 | 0.44 | -0.35 | 1.84 | 208 kg |
| BWD-i2 | 0.44 | +0.05 | 2.99 | 340 kg |
| BWD-t | 0.33 | +0.13 | 3.51 | 340 kg |
Backward simulations with a LPDM have successfully been established as a tool for the generation of source-receptor matrices. This method is attractive with respect to accuracy and computational effort. It yields useful results already in this first implementation, and it can be expected that it will improve after the following steps:
In the future, this method shall be applied also to other tracer experiments such as ANATEX and CAPTEX as well as to the release from the Chernobyl NPP disaster. Finally, the method may be used for unknown sources.
1 This work is a contribution to EUROTRAC-2 subprojects GLOREAM and GENEMIS-2, funded by FWF under P1295-GEO. Meteorological data from ECMWF have kindly been provided through ZAMG (Vienna).
2 The following derivation is restricted to the case without sinks. However, it is intuitively clear that the LPDM can be applied in the same way also for the linear case with sinks, as the relative magnitude of any loss along a certain transport distance is the same in the forward and in the backward mode. The exact derivation will be presented later.
