? ? ? (S,)

(D,) ϵ D

V⃗ Vτt δ(t−

τ)/Vτ

V δ V λ K t B < ϕ, ψ >=

∫_{R^{4}} ϕψdx⃗dt

Backtracking, Eulerian inverse transport and the adjoint transport equation for linear atmospheric trace species

Jean-Pierre Issartel, Frédéric Hourdin, Bertrand Cabrit, Olivier Talagrand, Abderra hmane Idelkadi

Brouillon: January 15, 2002

Introduction

Inverse problems are taking an increasing place in the study of the atmospheric transport of trace species (water, aerosols, chemical compounds, radionucleii). Classical questions concern the assimilation of observations of the atmospheric composition, the determination of uncertain emissions (cars, industry, natural CO2 sinks,...) or the determination of the location of an unusal event such as an accident on a nuclear plant.

General inverse technics, such as the variational assimilation introduced in the context of operational meteorology by Le Dimet et Talagrand (1986, see also, ??), generally apply to such problems. Variational technics most often imply the development of the adjoint of the direct model to compute the gradient of any cost function taken on the model output parameters with respect to input parameters (sources, initial composition, model parameters...).

We were recently involved in such studies in the frame of the Comprehensive Test Ban Treaty (CTBT). Eighty stations will be settled throughout the planet to monitor the atmospheric radioactivity from illicit nuclear tests. A backtracking approach, easily implemented on our global transport model, enabled to efficiently evaluate the monitoring acuity of the detectors for a statistically significant number of meteorological sequences (?). The idea of the ”detector-model” is now classical. The applications are often based on a qualitative interpretation of Lagrangian trajectories reconstructed back in time so as to identify the source of a tracer (??).

The specificity of the questions raised by the operation of the CTBT monitoring system led us to an in-depth investigation of the domain (partly presented by ?). We first present our approach as a generalization of Lagrangian backtracking, by the introduction of the reciprocity of atmospheric transport. This intuitive approach is valid for tracers transported along the same trajectories as the air including diffusive processes or a linear decay. It is equivalent, then, to the traditional adjoint approach.

Accordingly the adjoint transport equation of these tracers may be derived from physical consid erations without knowing anything about the adjoint algebra. For other linear trace species the intuitive approach fails. Je crois qu’il faut parler du retrotransport comme transport a rebours d’une sensibilit avant de parler de fonction importance An adjoint definition of backtracking is still possible as we shall see using the “importance function” from neutron transport theory. Then we shall explain that the statistical interpretation of diffusive inverse transport is valid for determining the origin of a particle. It is false when applied to macroscopic sources the uncertain identification of which is less random than biased or underdetermined for independent reasons. Finally we present an academic illustration and a more realistic one using the results of the ETEX campaign. Our investigation will be restricted here to the interpretation of single measurements. The joint interpretation of several measurements would require additional techniques.

It is surprising that the results described by former authors are poorly known and exploited. Pudykiewicz adressed the identification of sources as a sensibility problem and proposed a gene ral adjoint solution with a statistical management of the uncertainties. Earlier Uliasz and Pielke had even recognized the adjoint model to be an inverse model with the same diffusion as the direct model. In fact the adjustment of intuitive considerations with widely accepted theoretical arguments was still to be done. Intuitions are about the equivalence of backtracking with a problem of sensibility; about the necessarily uncertain identification of a source. Diffusion is subject to wrong intuitions concerning its impossible inversion in a backtracking approach or its statistical influence on dispersion. The clarification of these questions has been hindered furthermore by purely technical difficulties.

peut-etre annoncer le plan. Le dbut du paragraphe suivant est un peu abrupt

1 The air exchange rate and perfect tracers

We consider two volumes of air, a source S at some time t_{s}, and a detector D at a later time t_{d}. The mass of air contained in S at time t_{s} is denoted M(S, t_{s}), and the mass of air contained in D at time t_{d} is denoted M(D, t_{d}). The mass of air exchanged between S and D, i.e. the mass of air which is present in both S at time t_{s} and D at time t_{d}, is denoted M_{e}(S, t_{s}, D, t_{d}). We then define the exchangerate between the source and the detector as :

ε(S, t_{s}, D, t_{d}) = \frac{M_{e}(S, t_{s}, D, t_{d})}{M(S, t_{s})M(D, t_{d})}  (1)

The massic concentration in the detector of the fraction of air originating from the source is equal to M(S, t_{s})ε(S, t_{s}, D, t_{d}). Symmetrically, the massic concentration in the source of the fraction of air that will be present in the detector is equal to M(D, t_{d})ε(S, t_{s}, D, t_{d}).

We now consider a perfect tracer, i.e. a tracer that exactly follows the motion of transporting air. The amount of tracer present in a given mass of air is conserved in time. The air itself, as a self tracer, is a perfect tracer. We stress that the amount of tracer can be measured in many different ways. It can be measured as a mass, but it can be measured as well as a number of particles, as an electric charge, as a form of internal energy, or in still many other ways. The only requirement for what follows is that the amount must be additive, i.e. the amount of tracer contained in two non-overlapping volumes of air A and B together must be the sum of the amounts contained in A and B respectively. It is then seen that the exchange rate is the average massic concentration \overline{c}_{S,t_{s}}^{D,t_{d}} in D of a unit amount of tracer released in S uniformly with respect to air mass. Symmetricaly, it is also the average massic concentration \overline{c^{∗}}_{D,t_{d}}^{S,t_{s}} in S of a unit amount of tracer that would happen to be uniformly distributed in D with respect to air mass. The equality

\overline{c}_{S,t_{s}}^{D,t_{d}} = \overline{c^{∗}}_{D,t_{d}}^{S,t_{s}}  (2)

expresses what we will call the reciprocity of atmospheric transport. Throughout this paper we shall use the notation c for denoting the massic concentration of a tracer transported from past into future by the flow, and the notation c^{∗} for denoting the concentration of a ’retrotracer’ transported backward in time.

Remarque. Les mots ’plume’ and ’retroplume’ ne sont en fait pas utiliss dans l’article, au moins jusqu’ la Section 3.

Computation of the exchange rate requires only the knowledge of the velocity field v⃗(x⃗, t) as a function of space and time.

The reciprocity principle can also be stated in the following terms. The source volume S is advected by the flow and is transformed into the volume S_{td} at time t_{d}. The air exchanged between S and D is the air contained in the intersection S_{td} ⋂ D. Symmetrically, the detector volume D is transported backward in time into the volume D_{ts} at time t_{s}. Stating the reciprocity principle is just saying that the same mass of air, and therefore the same amount of perfect tracer, is contained in both S_{td} ⋂ D and S ⋂ D_{ts}.

Je verrais bien ici la figure que Frdric montre rgulirement.

2 Modeling backward transport

Pure advection of perfect tracers

We now look at how backward transport can be computed in practice, and consider first the case of pure inviscid advection. In the Lagrangian description of fluid motion, which is the most commonly used for identification of sources, the exchange rate can be evaluated by injecting a large number of particles, uniformly with respect to air mass, in S at time t_{s}, advecting those particles with the flow, and counting those that reach D at time t_{d}. One can alternatively start with a large number of particles in D at time t_{d}, follow them back in time and count those that originated in S at time t_{s}.

Remarque. Dans le paragraphe qui suit, j’ai t toute rfrence la dynamique de contours qui, contrairement ce que semblait dire le texte originel, correspond une description lagrangienne, et non eulrienne, du mouvement. Reponse de FH : pour moi, c’est une approche mi chemin

In the Eulerian description, one is led to consider the massic concentration of the tracer, which is conserved for each fluid particle. The corresponding evolution equation reads

\frac{∂c}{∂t} + v.grad c = σ(x⃗, t)  (3)

where σ(x⃗, t) is the source of tracer. In the present case where the source in uniform in S and concentrated at time t_{s}, one is led to take σ(x⃗, t) = δ(t−t_{s})/M(S, t_{s}) in S, 0 outside (δ is the mass of Dirac). Eq. ?? is then integrated from the initial condition c(x⃗, t) = 0 at some time t < t_{s}). The exchange rate is then the average over D, with respect to air mass, of the solution c(x⃗, t_{d}) at time t_{d}.

In the backward option, the equation

−\frac{∂c^{∗}}{∂t} − v.grad c^{∗} = π(x⃗, t)  (4)

where π(x⃗, t) = δ(t_{d} − t)/M(D, t_{d}) in D, 0 outside ,, is integrated from the ’final’ condition c^{∗}(x⃗, t) = 0 at some time t > t_{d}).

We have considered a source and a detector that are instantaneous in time. Eq. ?? and Eq. ?? still stand for non-instantaneous sources and detectors. It is sufficient to replace σ(x⃗, t) and π(x⃗, t) with terms which are appropriate functions of time. For instance, a pump detector localised at a point in space will be represented by a function which is a dirac in space and extends over a time interval.

Question. Le coefficient d’change reste-t-il bien dfini dans ce cas?

Turbulent mixing of perfect tracers

In numerical practice, Eq. ?? must be discretized and can explicitly describe only the scales of the motion which are larger than the discretization scale. Diffusion operators are introduced in numerical models in order to describe the statistical effect of the unresolved scales of motion on the resolved scales. The resulting ’turbulent mixing’ is particularly important for planetary atmospheres in the so-called planetary boundary layer. As a result of diffusion, a tracer emitted by a source is dispersed over a wider area, and its concentration is decreased. Symmetrically a detector can detect a tracer originating from a wider area, but with a weakened sensitivity. These succinct considerations show that diffusion must remain diffusive for computation of the exchange coefficient through backward integration.

More precisely, the reciprocity relation applies to the complete transport, including its small scale unresolved component. In the ergodic approach turbulent motions are described as a random process.

Question; Est-ce que l’ergodicit avraiment quellquechose voir ici? Rponse de Frdric : je ne crois pas

If unresolved motions are assumed to be statistically symmetric in time, they must be parameterised by the same operator, denoted here ζ, in the forward and backward transport equations :

\frac{∂c}{∂t} + v⃗.grad c + ζ(c) = σ  (5)
−\frac{∂c^{∗}}{∂t} − v⃗.grad c^{∗} + ζ(c^{∗}) = π  (6)

This result may seem paradoxical, and the reader may legitimately wonder if the law of increase of entropy is not violated in the backward integration of Eq. ??. The answer lies in the fact that the backward integration is not intended at determining the original spatial position of each particle of tracer, and is thus not an exact reconstruction of the past from the future. The individuality of tracer particles, conserved in Eq. ??, is lost in diffusion. What the backward integration of Eq. ?? is intended to produce is only the amount of tracer that was emitted by the source at time t_{s} and present in the detector at time t_{d}. Everything else being the same, that amount decreases with increasing difference t_{d}-t_{s}, which shows that diffusion must remain diffusive for bacwrad integration of Eq. ??.

As an example, we discuss these aspects in more detail in the classical description of eddy-diffusivity, as built from mixing length theory. For simplicity, we consider only vertical diffusion, which dominates in the boundary layer. The developments below could easily be extended to multidimensional diffusion.

Remarque. Je ne suis pas sr que le paragraphe suivant, que je laisse ici tel quel, soit vraiment ncessaire. Il n’est en tout cas pas sa place.

Though simple this result is perplexing because, as diffusion is irreversible, it seems to run counter the second principle. The obstacle is only apparent. The atmosphere is an open system that permanently refreshes its organization while degrading so lar into thermal radiation cooled roughly from 5000K down to 260K. The coherence of a sample of air exists only through and by the time of the sampling process. According to Eq. ?? this coherence was not yet formed in the past and according to Eq. ?? i t will be lost in the future.

Parameterizations of eddy-diffusivity are based upon the hypothesis that unresolved turbulence consists of statistically symmetric upward and downward motions. Those motions produce a diffusion of atmospheric tracers just as random walk of molecules produces molecular diffusion. When inverting the direction of time and motions, turbulence still consists of upward and downward motions which have the same statistical effect as in the forward march in time.

The proper introduction of turbulent closure relies on the notion of ensemble average. The ensemble average \overline{q} is defined as the average of a quantity q taken over a set of random realisations of turbulent motions. From ergodic considerations, no distinction is generally made between ensemble average, spatial average at the resolution of the model or time average. For compressible flows, it is convenient to introduce a mass weighted average q̃ = \overline{ρq}/\overline{ρ} and perturbations q′ = q − q̃. Then \overline{ρq′} = 0

Remarque. Cette dernire galit n’est pas valide pour des moyennes spatiales ou temporelles, qui sont des moyennes glissantes.

Introducing the equation for conservation of the mass of transporting air

\frac{∂ρ}{∂t} + div (ρv) = 0  (7)

Eq. ?? is first transformed into its conservative form

\frac{∂ρc}{∂t} + div (ρvc) = ρσ(x⃗, t)  (8)

Taking the ensemble average of this equation and noting that

\overline{ρvc} = \overline{ρ}ṽc̃ + \overline{ρv′c′},  (9)

lead to

\frac{∂\overline{ρ}c̃}{∂t} + div (\overline{ρ}ṽc̃) + div (\overline{ρv′c′}) = \overline{ρ}σ(x⃗, t)  (10)

The variables \overline{ρ}, ṽ and c̃ are ensemble averages. In view of the ergodic interpretation mentioned above, they can be considered as local space-time averages, and are thus somewhat similar to the explicit variables of numerical models. From now, they will be denoted ρ, v and c.

Introducing once again the mass conservation, (valid in the same form for an individual realisation and for the ensemble average), one obtains the averaged advection equation

\frac{∂c}{∂t} + v.grad c + \frac{1}{ρ}div (\overline{ρv′c′}) = σ(x⃗, t)  (11)

According to the hypothesis we have made above, we assume the vertical component \overline{ρw′c′} of the turbulent transport div (\overline{ρv′c′}) contributes significantly to the diffusion. The mixing length view is that the value of c in subsiding air is representative of the air at a typical distance l (the mixing length) above. Thus, c′ = \overline{c}(z + l) − \overline{c}(z) for downard motions (w′ < 0) and c′ = \overline{c}(z − l) − \overline{c}(z) for upward motions, resulting for both cases in

\overline{ρw′c′} = −ρK_{z} \frac{∂\overline{c}}{∂z}  (12)

where K_{z} ≃ \overline{∥w′|}l is the eddy-diffusivity. The transport equation thus finally reads

\frac{∂c}{∂t} + v.grad c − \frac{1}{ρ} \frac{∂}{∂z} (ρK_{z} \frac{∂c}{∂z} ) = σ(x⃗, t) (13)

Note that the ”closure” consisting in expressing K_{z} as a function of the large scale meteorological variables can be quite complex.

Inverting the direction of w′ in the above analysis changes nothing as soon as upward and downward turbulent mass fluxes are statistically symmetric. So exactly the same diffusion operator must be kept, without changing sign, for the backward integration. The same analysis would of course apply to horizontal components.

“Local closure” or “eddy-diffusion”, as just described, is not the only way of parameterizing subgrid scale transport. In particular, state-of-the-art parameterizations of convection are based upon an explicit representation of convective mass fluxes. To be precise, those “mass flux” parameterizations are generally sub-grid scale on the horizontal but not on the vertical. In this kind of parameterizations, the transport is disymmetric. For instance, there is generally a distinct representation of strong up-drafts, with vertical winds having values from a few to tens of meters per second, which occur in the center of convective towers. These updrafts are balanced by slower sinking, either nearby the column, or at larger scales. In those parameterizations, the sign of mass fluxes must be changed explicitly in the backward integration : for instance, the air must be carried down to the surface very fast in the center of the tower. However, even convective mixing tends to diffuse tracer both in the future and in the past.

Extension to non conservative tracers

The reciprocity of atmospheric transport extends to non conservative tracers following air motion with a linear decay or growth:

\frac{∂c}{∂t} + v⃗.grad c = −λc  (λ > 0 for a decay) (14)

λ can be for instance the inverse of the radioactive time-life of a radionuclide, a rate of scavenging by rain, or still the rate of a chemical reaction. In the last two cases, λ(x⃗, t) will in general be a function of both space and time. What is important here is that λ must be independent of the tracer concentration.

In order to define an exchange rate that takes into account the variations of tracer concentration along a trajectory, it is necessary to modify appropriately the ’exchanged mass’, so as to increase or decrease it in proportion with the increase or decrease of tracer concentration. This leads to the following definition, in which the numerator is modified from Eq. ??. Still defining the exchange rate as the concentration in D at t_{d} resulting from the release of a unit tracer in S at t_{s}, it reads as:

ε (S, t_{s}, D, t_{d}) = \frac{1}{M (S, t_{s}) M (D, t_{d})} × (15)
∫_{γ in Γ_{S,D}} exp [− ∫_{t_{s}}^{t_{d}} λ (γ, t) dt] ρ(x⃗, t_{s})dΩ_{S}

where Γ_{S,D} is the ensemble of trajectories originating in S at time t_{s} and ending in D at time t_{d}, Ω_{S} is an elementary volume in S at the origin of the trajectory γ and λ (γ, t) is the value of λ at time t along γ. Since the integral bears on trajectories connecting (S, t_{s}) and (D, t_{d}), the expression above is unmodified if the mass element ρ(x, t_{s})dΩ_{S} in the source is changed into a mass element ρ(x, t_{d})dΩ_{D} in the detector.

Again there are two reciprocal interpretations of the exchange rate. It can first be computed as the mean concentration in (D, t_{d}) of a unit of tracer uniformly injected in (S, t_{s}). But it can be computed as well as the mean concentration in (S, t_{s}) of a unit of retro-tracer uniformly injected in (D, t_{d}). The symmetry of Eq. ?? shows that the same attenuation or increase must be applied in the backward integration as in the forward integration. A decay of the direct tracer in the forward integration imposes a decay of the retro-tracer in the backward integration, and not the growth that would restore the original source.

Once again, we obtain a result that may seem paradoxical. And here again, the answer lies in the fact that the output of the backward integration is not the original amount of tracer, but the amount emitted by the source at time t_{s} and still present in the detector at time t_{d}. As in the case of diffusion, that amount decreases with increasing difference t_{d}-t_{s}, which clearly shows that the concentration of retrotracer must decrease in the course of the backward integration.

Numerical implementation

The exchange coeficient between a source and a detector, as it results from large scale advection, vertical turbulent mixing and linear decay or growth, can finally be computed either by forward in time integration of the equation

\frac{∂c}{∂t} + v.grad c − \frac{1}{ρ} \frac{∂}{∂z} (ρK_{z} \frac{∂c}{∂z} ) + λc = σ(x⃗, t) (16)

or by backward integration of the equation

−\frac{∂c^{∗}}{∂t} − v.grad c^{∗} − \frac{1}{ρ} \frac{∂}{∂z} (ρK_{z} \frac{∂c^{∗}}{∂z} ) + λc^{∗} = π(x⃗, t)

(17)

Eq. ?? is obtained from Eq. ?? by just changing the signs of velocity v and of time t. This suggests that, given a numerical model for Eq. ??, the corresponding backward integration could be performed by just changing the signs of v and t in the model. Because however of truncation errors, a discretized model of Eq. ?? may not possess the reciprocity property which is at the basis of the derivation of Eq. ??.

Anticipating on the discussion of the following section, discretizing the adjoint equation is not necessarily equivalent to taking the adjoint of the discretized direct equation. Considering for instance advection, backward trajectories obtained by just changing the signs of v and t in the discretized direct model may not be identical with forward trajectories.

Remarque. Je ne sais pas trs bien quoi dire ici propos de la diffusion. Les arguments qui prcdent suggrent qu’il suffit, pour assurer la rciprocit, de garder la mme diffusion pour l’intgration rebours que pour l’intgration directe, alors que l’analyse adjointe conduit une condition beaucoup plus forte, savoir que l’oprateur de diffusion soit auto-adjoint pour le produit scalaire . Reponse de FH : dans le cas d’un modle symmtrique, c’est auto-adjoint. Dans le cas d’un modle dissymtrique, il faut prendre l’adjoint qui renverse la dissymtrie.

A very straightforward application of this backward transport equation was used to evaluate the detection capability of the CTBT atmospheric network. For a given 1 kt test (specification of CTBT), or 10^{15} Bq of Xe133 (one Becquerel corresponds to one disintegration per second), the test was detected at the station if the concentration was above 1 mBq/kg. In the model, we rather released the 10^{15} Bq at a station of the network, integrated Eq. ??. The detection was said positive for a source in S at t_{s} if the concentration c^{∗} was larger than 1 mBq/kg at that point thus exactly reversing the role of source and detector.

2.1 Illustrations

La Fig. ?? montre un test numrique de cette rciprocit du transport atmosphrique sur la campagne ETEX. Plus prcisment, on a trac pour neufs des stations, la concentration observe au cours des 72 heures de la campagne et la concentration simule avec les approches directe et rtro. Pour les calculs inverses, les coefficients d’change sont d’abord calculs en injectant continuement sur 3 heures (Δt) partir de t_{k} = t_{0} + kΔt un traceur unit pour chaque station D_{i}. Le transport est calcul rebours dans le temps jusqu’ t_{0}. Pour chaque couple (i, k), on conserve le coefficient d’change χ_{i,k}, valeur moyenne entre t_{0}+12 h et t_{0} de la concentration du traceur dans la maille de Monterfile (d’o a t mis le PMCH). Finallement, la concentration “retro” montre sur la Fig. ?? (pointills) est simplement c^{_{i}(t_{k})=q×χ_{i,k}}.

Comme on l’a dit plus haut, l’accord entre observations et simulations directes est comparable aux meilleures simulations prsentes dans les tudes d’intercomparaison (??). Cependant, ce qui nous in-