\documentclass{article} \usepackage{a4} \usepackage{epsf} \usepackage{graphicx} \usepackage[round]{natbib} \bibliographystyle{apalike} \tolerance=8000 \pretolerance=4000 \newenvironment{listdot}{\begin{list}{$\bullet$}{\setlength{\itemsep} {0ex} \setlength{\parsep}{0ex} \leftmargin 1em}}{ \end{list}} \newenvironment{listcar}{\begin{list}{$\square$}{\setlength{\itemsep} {0ex} \setlength{\parsep}{0ex} \leftmargin 1em}}{ \end{list}} % Final simple colonne % \textwidth 17cm % \oddsidemargin 0cm % Final double colonne \setlength{\oddsidemargin}{-4mm} \setlength{\evensidemargin}{0mm} \setlength{\columnsep}{10mm} \setlength{\textwidth}{175mm} \setlength{\textheight}{229mm} \twocolumn \title{Backtracking, Eulerian inverse transport and the adjoint transport equation for linear atmospheric trace species \author{Jean-Pierre Issartel, Fr\'ed\'eric Hourdin, Bertrand Cabrit, Olivier Talagrand, Abderra hmane Idelkadi} \date{Brouillon: \today} } \begin{document} \def\A{{S}} \def\B{{D}} \def\V{{\bf v}} \def\ta{t_s} \def\tb{t_d} \def\M{M} \def\Mc{M^{\mathsf{ex.}}} \def\ME{\mathcal{M}} \def\MEc{\mathcal{M}^{\mathsf{ex.}}} \def\q{c} \def\m{q} \def\S{\xi} \def\aaa{{(\bf A)}} \def\bbb{{(\bf B)}} \def\C{\varepsilon} \def\excoeff{\varepsilon} \def\cd{cd} \def\cs{cs} \def\moy#1{\overline{#1}} \def\exbar{\bar{\varepsilon}} \def\exrbar{\bar{\varepsilon}^*} \def\grad{{\mbox{\bf grad}}} \def\der#1\def\ex{\varepsilon} \def\ex{\varepsilon} \def\exr{\varepsilon^*} \def\xe{$^{133}$Xe } \def\ba{${ ^{140}}$Ba } \def\dix#1{10$^{#1}$} \def\microbq{$\mu$Bq } \def\m#1{$^{-#1}$} % #2{\frac{\partial #1}{\partial #2}} \def\dep#1{\left(#1\right)} \def\depb#1{\left[#1\right]} \def\dem{1/2} \def\eq#1{Eq.~\ref{eq:#1}} \def\fg#1{Fig.~\ref{fg:#1}} \def\sec#1{Section~\ref{sec:#1}} \def\dq{{\dep{\delta q}}} \def\S{S} \def\D{D} \def\ts{t_s} \def\td{t_d} \def\vs{V_s} \def\vd{V_d} \def\ms{M_s} \def\md{M_d} \def\M{M} \def\Mc{M^{\mathsf{ex.}}} \def\ME{\mathcal{M}} \def\MEc{\mathcal{M}^{\mathsf{ex.}}} \def\c{c} \def\q{q} \def\dt#1{\frac{\partial #1}{\partial t}} \def\C{\varepsilon} \def\k{\kappa} \def\l{\lambda} \def\g{\gamma} \def\se{\sigma} \def\pe{\pi} \def\point{\vec{c}} \def\mesure{\mu} \def\source{\sigma} \def\inttemps{\int_{-\infinty}+{\infinty}} \def\intespacetemps{\int_{\Omega \times R}} \def\ex{c} \def\exr{c^*} \def\x{\vec{x}} \def\kz{K_z} \def\diffvert#1{\frac{1}{\rho}\frac{\partial}{\partial z}\dep{\rho \kz \frac{\partial#1}{\partial z}}} \def\intespace{\int_\Omega} \def\dv{d{\vec{x}}} \def\vent{\vec{V}} \def\div#1{{\mbox{div}}\dep{#1}} \def\ex{c} \def\exr{c^*} \def\afaire#1{{\bf #1}} \maketitle %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection*{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 \cite[1986, see also, ][]{Tala:87,Cour:87},\nocite{LeDi:86} 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 \cite[]{Hour:00}. 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 \cite[]{Ramonet, n'co}. The specificity of the questions raised by the operation of the CTBT monitoring system led us to an in-depth investigation of the domain \cite[partly presented by ][]{Hour:00}. 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. \afaire{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. \afaire{peut-etre annoncer le plan. Le début du paragraphe suivant est un peu abrupt} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{The air exchange rate and perfect tracers} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \def\c#1#2{{\overline{c}}_{#1}^{#2}} \def\cstar#1#2{{\overline{c^*}}_{#1}^{#2}} We consider two volumes of air, a $source$ $S$ at some time $\ts$, and a {\em detector} $D$ at a later time $\td$. The mass of air contained in $S$ at time $\ts$ is denoted $M(S, \ts)$, and the mass of air contained in $D$ at time $\td$ is denoted $M(D, \td)$. The mass of air exchanged between $S$ and $D$, $i. e$. the mass of air which is present in both $S$ at time $\ts$ and $D$ at time $\td$, is denoted $M_e(S, \ts, D, \td)$. We then define the $exchange rate$ between the source and the detector as : \begin{equation} \label{excoef} \excoeff(S, \ts, D, \td) = \frac{M_e(S, \ts, D, \td)}{M(S, \ts) M(D, \td)} \end{equation} The massic concentration in the detector of the fraction of air originating from the source is equal to $M(S, \ts) \excoeff(S, \ts, D,\td)$. Symmetrically, the massic concentration in the source of the fraction of air that will be present in the detector is equal to $M(D, \td) \excoeff(S, \ts, D, \td)$. 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 $\c{S,\ts}{D,\td}$ 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 $\cstar{D,\td}{S,\ts}$ in $S$ of a unit amount of tracer that would happen to be uniformly distributed in $D$ with respect to air mass. The equality \begin{equation} \label{recipro} \c{S,\ts}{D,\td}=\cstar{D,\td}{S,\ts} \end{equation} expresses what we will call the {\em 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. \afaire{ Remarque. Les mots 'plume' and 'retroplume' ne sont en fait pas utilisés dans l'article, au moins jusqu'à la Section 3.} Computation of the exchange rate requires only the knowledge of the velocity field $\vec{v(\vec{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 $\td$. The air exchanged between $S$ and $D$ is the air contained in the intersection $S_{td}\bigcap D$. Symmetrically, the detector volume $D$ is transported backward in time into the volume $D_{ts}$ at time $\ts$. 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}\bigcap D$ and $S \bigcap D_{ts}$. \afaire{Je verrais bien ici la figure que Frédéric montre régulièrement.} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Modeling backward transport} \subsection*{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 $\ts$, advecting those particles with the flow, and counting those that reach $D$ at time $\td$. One can alternatively start with a large number of particles in $D$ at time $\td$, follow them back in time and count those that originated in $S$ at time $t_s$. \afaire{ Remarque. Dans le paragraphe qui suit, j'ai ôté toute référence à la dynamique de contours qui, contrairement à ce que semblait dire le texte originel, correspond à une description lagrangienne, et non eulérienne, 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 \begin{equation} \label{eq:d} \frac{\partial \ex}{\partial t}+\V.\grad\ \ex = \sigma(\vec{x},t) \end{equation} where $\sigma(\vec{x}, t)$ is the source of tracer. In the present case where the source in uniform in $S$ and concentrated at time $\ts$, one is led to take $\sigma(\vec{x}, t)= \delta(t-\ts)/M(S, \ts)$ in $S$, $0$ outside ($\delta$ is the mass of Dirac). \eq{d} is then integrated from the initial condition $c(\vec{x}, t)$ = 0 at some time $t<\ts$). The exchange rate is then the average over $D$, with respect to air mass, of the solution $c(\vec{x}, \td)$ at time $\td$. In the backward option, the equation \begin{equation} \label{eq:recipro} -\frac{\partial \exr}{\partial t}-\V.\grad\ \exr = \pi (\vec{x},t) \end{equation} where $\pi(\vec{x}, t)= \delta(\td-t)/M(D, \td)$ in $D$, $0$ outside\\,, is integrated from the 'final' condition $c^*(\vec{x}, t)$ = 0 at some time $t>\td$). We have considered a source and a detector that are instantaneous in time. \eq{d} and \eq{recipro} still stand for non-instantaneous sources and detectors. It is sufficient to replace $\sigma(\vec{x},t)$ and $\pi (\vec{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. \afaire{Question. Le coefficient d'échange reste-t-il bien défini dans ce cas?} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection*{Turbulent mixing of perfect tracers} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% In numerical practice, \eq{d} 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. \afaire{Question; Est-ce que l'ergodicité avraiment quellquechose à voir ici? Réponse de Frédéric : 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 $\zeta$, in the forward and backward transport equations : \begin{equation}\label{eq:dir} \frac{\partial c}{\partial t} + \vec{v}.grad \ c + \zeta(c) = \sigma \end{equation} \begin{equation}\label{eq:inv} - \frac{\partial c^*}{\partial t} - \vec{v}.grad \ c^* + \zeta(c^*) = \pi \end{equation} 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{inv}. 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{d}, is lost in diffusion. What the backward integration of \eq{inv} is intended to produce is only the amount of tracer that was emitted by the source at time $\ts$ and present in the detector at time $\td$. Everything else being the same, that amount decreases with increasing difference $\td$-$\ts$, which shows that diffusion must remain diffusive for bacwrad integration of \eq{inv}. 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 sûr que le paragraphe suivant, que je laisse ici tel quel, soit vraiment nécessaire. 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 $5000 K$ down to $260 K$. The coherence of a sample of air exists only through and by the time of the sampling process. According to \eq{inv} this coherence was not yet formed in the past and according to \eq{dir} 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. \def\ave#1{\overline{#1}} \def\mave#1{\tilde{#1}} The proper introduction of turbulent closure relies on the notion of ensemble average. The ensemble average $\ave{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 $\mave{q}=\ave{\rho q}/\ave{\rho}$ and perturbations $q'=q-\mave{q}$. Then $\ave{\rho q'}=0$ Remarque. Cette dernière é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 \begin{equation} \label{eq:direct} \frac{\partial \rho}{\partial t}+\div{\rho \V} = 0 \end{equation} \eq{direct} is first transformed into its conservative form \begin{equation} \label{eq:direct} \frac{\partial \rho\ex}{\partial t}+\div{\rho \V \ex} = \rho\sigma(\vec{x},t) \end{equation} Taking the ensemble average of this equation and noting that \begin{equation} \ave{\rho \V \ex}=\ave{\rho}\mave{\V}\mave{\ex}+\ave{\rho \V'\ex'}, \end{equation} lead to \begin{equation} \label{eq:direct} \frac{\partial \ave{\rho}\mave{\ex}}{\partial t} +\div{\ave{\rho}\mave{\V}\mave{\ex}} + \div{\ave{\rho\V'\ex'}} = \ave{\rho}\sigma(\vec{x},t) \end{equation} The variables $\ave{\rho}$, $\mave{\V}$ and $\mave{\ex}$ 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 $\rho$, $\V$ and $\ex$. 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 \begin{equation} \label{eq:direct} \frac{\partial \ex}{\partial t}+\V.\grad\ \ex +\frac{1}{\rho}\div{\ave{\rho\V'\ex'}} = \sigma(\vec{x},t) \end{equation} \def\pert#1{#1'} According to the hypothesis we have made above, we assume the vertical component $\ave{\rho w'\ex'}$ of the turbulent transport $\div{\ave{\rho\V'\ex'}}$ contributes significantly to the diffusion. The mixing length view is that the value of $\ex$ in subsiding air is representative of the air at a typical distance $l$ (the mixing length) above. Thus, $\pert{\ex}=\moy{\ex}(z+l)-\moy{\ex}(z)$ for downard motions ($w'<0$) and $\pert{\ex}=\moy{\ex}(z-l)-\moy{\ex}(z)$ for upward motions, resulting for both cases in \begin{equation} \ave{\rho w'\ex'}=-\rho K_z \frac{\partial \moy{\ex}}{\partial z} \end{equation} where $K_z\simeq \ave{\| w' |} l$ is the eddy-diffusivity. The transport equation thus finally reads \begin{equation} \label{eq:directb} \frac{\partial \ex}{\partial t}+\V.\grad\ \ex -\frac{1}{\rho}\frac{\partial}{\partial z} \dep{\rho K_z \frac{\partial \ex}{\partial z}}= \sigma (\vec{x},t) \end{equation} 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. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection*{Extension to non conservative tracers} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The reciprocity of atmospheric transport extends to non conservative tracers following air motion with a linear decay or growth: \begin{equation} \frac{\partial c}{\partial t} + \vec{v}.grad \ c =-\lambda c \hspace{10mm} ( \lambda > 0 \mbox{ for a decay}) \end{equation} $\lambda$ 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, $\lambda(\vec{x}, t)$ will in general be a function of both space and time. What is important here is that $\lambda$ 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{excoef}. Still defining the exchange rate as the concentration in $D$ at $\td$ resulting from the release of a unit tracer in $S$ at $\ts$, it reads as : \begin{eqnarray} &&\varepsilon\dep{S,\ts,D,\td}=\frac{1}{\M\dep{S,\ts}\M\dep{D,\td}}\times\\ &&\int_{\gamma \mbox{\ in\ } \Gamma_{S,D}} {\exp\depb{-\int_{\ts}^{\td}{\lambda\dep{\gamma,t} dt} }\rho(\vec{x},\ts) d\Omega_S} \nonumber \end{eqnarray} where $\Gamma_{S,D}$ is the ensemble of trajectories originating in $S$ at time $\ts$ and ending in $\D$ at time $\td$, $\Omega_S$ is an elementary volume in $S$ at the origin of the trajectory $\gamma$ and $\lambda\dep{\gamma,t}$ is the value of $\lambda$ at time $t$ along $\gamma$. Since the integral bears on trajectories connecting $(S,\ts)$ and $(D,\td)$, the expression above is unmodified if the mass element $\rho(x,\ts)d\Omega_S$ in the source is changed into a mass element $\rho(x,\td)d\Omega_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,\td)$ of a unit of tracer uniformly injected in $(S,\ts)$. But it can be computed as well as the mean concentration in $(S,\ts)$ of a unit of retro-tracer uniformly injected in $(D,\td)$. The symmetry of \eq{eqnarray} 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 $\ts$ and still present in the detector at time $\td$. As in the case of diffusion, that amount decreases with increasing difference $\td$-$\ts$, which clearly shows that the concentration of retrotracer must decrease in the course of the backward integration. \subsection*{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 \begin{equation}\label{eq:directc} \frac{\partial \ex}{\partial t}+\V.\grad\ \ex -\diffvert{\ex}+\lambda\ex=\sigma (\vec{x},t) \end{equation} or by backward integration of the equation \begin{equation}\label{eq:reciproc} -\frac{\partial \exr}{\partial t}-\V.\grad\ \exr -\diffvert{\exr}+\lambda\exr= \pi (\vec{x},t) \end{equation} \eq{reciproc} is obtained from \eq{directc} by just changing the signs of velocity $\V$ and of time $t$. This suggests that, given a numerical model for \eq{directc}, 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{directc} may not possess the reciprocity property which is at the basis of the derivation of \eq{reciproc}. 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. \afaire{ Remarque. Je ne sais pas très bien quoi dire ici à propos de la diffusion. Les arguments qui précèdent suggèrent qu'il suffit, pour assurer la réciprocité, de garder la même diffusion pour l'intégration à rebours que pour l'intégration directe, alors que l'analyse adjointe conduit à une condition beaucoup plus forte, à savoir que l'opérateur de diffusion soit auto-adjoint pour le produit scalaire . Reponse de FH : dans le cas d'un modèle symmétrique, c'est auto-adjoint. Dans le cas d'un modèle dissymétrique, il faut prendre l'adjoint qui renverse la dissymétrie.} 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{reciproc}. The detection was said positive for a source in $S$ at $\ts$ if the concentration $\exr$ was larger than 1~mBq/kg at that point thus exactly reversing the role of source and detector. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Link with adjoint methods} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection*{The perfect case} The exchange rate obviously describes the sensibility of a perfect tracer measurement in $(D,t_d)$ to a source in $(S,t_s)$. The convergence of ideas of backward transport and sensibility has been mentionned long ago (CITATIONS...). It seems nevertheless that purely mathematical obstacles have hindered the clear formulation of this convergence with its limitations. The measurement operates as a scalar product on the tracer concentration and the sampling distribution. When $c$, $\sigma$, $c^*$, $\pi$ are all referred to the unit mass of air according to the logic of the exchange rate, the analytic form of the measurement product is : \begin{equation} \label{proscal} \mu (c,\pi) = \int \rho(\vec{x}, t) \ c(\vec{x},t) \ \pi(\vec{x},t) \ d \vec{x} \ dt \end{equation} The dispersion law of the perfect tracer is in fact the dispersion law of the air. Completed with appropriate boundary conditions, it leads to a one to one relation $c=L(\sigma)$ between concentration and source fields. The reciprocity relation then writes as : \begin{equation} \int \rho \ L(\sigma) \ \pi(\vec{x},t) \ d \vec{x} \ dt = \int \rho \ \sigma(\vec{x},t) \ {c^*}(\vec{x},t) \ d \vec{x} \ dt \end{equation} From this relation, valid for any source and detector volumes $S$, $D$, we conclude : $c^* = L^*(\pi)$. $L^*$ is the operator adjoint to $L$ for the measurement product. In paragraph 2 we showed how to obtain an inverse transport law out of physical considerations. We now obtain the same law as an adjoint law $L^*$. An analogous conclusion clearly stands for non-conservative tracers following air masses so that \eq{directc} and \eq{reciproc} now appear to be adjoint equations. This double derivation of the inverse law of perfect tracers is a property of the transport fluid, of the air. In particular the diffusion operator is self-adjoint : \begin{equation} \label{eq:difadj} \zeta = \zeta^* \end{equation} \eq{difadj} is no longer a specific property of inverse transport, it is a general property of diffusion.\\ The physical link between backward transport and sensibility analysis does not depend on conventions, but its transparency does and especially the self-adjoint character of the diffusion. Diffusion is not self-adjoint for badly chosen conventions such as follows. Generally the measurement is represented as the simpler product $\mu(C,\pi)=\int C \pi d\vec{x} dt$ which amounts to unsymmetrically referring forward variables $C=\rho c$, $\Sigma=\rho \sigma$ to the unit volume and backward variables $c^*$, $\pi$ to the unit mass of air. \afaire{Il faut adoucir un peu la critique.} Pudykiewicz on the one hand, Uliasz and Pielke on the other hand used the analytic form $\mu(C,\pi)=\int C \pi d\vec{x} dt$, all with negative consequences. The first used a sound parameterisation of diffusion. The symmetry of normal and adjoint transport was accordingly kept hidden to him and he did not clearly recognize his adjoint tracer to be a retrotracer. The two latter started from their understanding of the symmetry of forward and backward transports. They wrongly required their parameterisation of diffusion to be self-adjoint for $\int C \pi d\vec{x} dt$. \subsection*{The importance function} \afaire{L'interprétation du rétrotransport comme le transport à rebours dans le temps de la sensibilité de la mesure à la source s'étend au-delà du transport advectif diffusif. Il me semble que ca doit etre dit avant de parler de la fonction importance.} The straightforward definition of backward transport is a privilege of perfect tracers or of tracers that, at least, follow air trajectories. Otherwise the definition of an exchange rate according to material transfers is irreparably lost. Fortunately, for all linear tracers, the adjoint concept of ``importance'' from reactor physics replaces almost exactly the failing exchange rate. Ideas were introduced early in neutron transport theory by Weinberg, Wigner, Hurwitz, and named by Soodak. The definition below is taken from Lewins [?, 1960?] : \\ \\ Definition: {\it The importance} $c^*(\vec{x}, t)$ {\it is defined as the expected or probable contribution of one particule at} $\vec{x}$ {\it at time t to the meter reading. Thus a particle is ``important'' to the (future) observable reading.}\\ In this definition a tracer particle is not distinguished from its ``progeny'', that is to say all particles generated from it through non-conservative processes. Of course atmospheric tracers are not neutrons multiplied by collisions with big nuclei, their 'progeny' would be rather negative. The choice of notation $c^*$ is coherent with the previous developpements about the perfect tracers. It emerges readily that the exchange rate of perfect tracers $\varepsilon(S, \ts, D, \td)$ represents the importance of a particle from $S$ in $D$ as a sample at time $\td$. Though it is not obvious, the value of the importance function does depend on conventions chose n to describe tracer concentration, tracer source and sampling distribution. We shall still refer them to the unit mass of air though for non perfect tracers the theoretical interest of this convention vanishes. The practical interest often remains as there is generally a dominating perfect part in the dispersion equation. The particle at $\vec{x}$ at time $t$ is described as a source by the distribution: \begin{equation} {\sigma}_{\vec{x}, t}(\vec{x}', t') = \frac{\delta(\vec{x}'-\vec{x}, t'-t)}{\rho(\vec{x}, t)} \end{equation} If sources were referred to the unit volume, the same particle would be described by ${\Sigma}_{\vec{x}, t}(\vec{x}', t') = \delta(\vec{x}'-\vec{x}, t'-t)$. We obtain the following importance referred to the unit mass of air : \begin{eqnarray} \label{eq:et1} c^*(\vec{x}, t) &=& \int \rho \ L({\sigma}_{\vec{x}, t}) \ \pi(\vec{x}', t') \ d \vec{x}'\ dt' \\ \label{eq:et2} &=& \int \rho \ {\sigma}_{\vec{x}, t}(\vec{x}', t') \ L^*(\pi) \ d \vec{x}'\ dt' \end{eqnarray} that is to say : $ c^* = L^*(\pi)$. Hence the importance function can be calculated at once as the dispersion of an adjoint tracer. The adjunction is taken through the measurement product. The dispersion law and boundary conditions are adjoint to those of the tracer. Again $c$ and $c^*$ are tied by the reciprocity relation. Taking this relation as the general definition of inverse transport we shall say that the adjoi nt tracer is a retrotracer and we shall still use the word 'retroplume' for $c^*$. Air sampling distribution $\pi$ becomes the release distribution for the retrotracer. This amounts to saying the retrotracer is exactly the sample of air in the detector and this ai r has reached the detector from everywhere according to the tracer dispersion law instead of th e air dispersion law.\\ Let us consider as an example the following direct and inverse models for a trace species subje ct to deposition with a drift speed $\vec{d}(\vec{x},t)$ \begin{equation} \frac{\partial c}{\partial t} + (\vec{v} + \vec{d}).\vec{grad} \ c=\sigma(\vec{x},t) \end{equation} \begin{equation} -\frac{\partial c^*}{\partial t} - \vec{v}.\vec{grad} \ c^* - \frac{div \ \rho c^* \vec{d} }{\rho} = \pi(\vec{x},t) \end{equation} The adjoint concentration or importance function $c^*$ does describe the degree of involvement from each place $\vec{x}$ and date $t$ into the measurement $\mu$. This description does correspond to backtrajectories weighted in such a way that, if the source of tracer was a spot $q {\sigma}_{\vec{x}, t}$ in $\vec{x}$ at $t$, the total amount of tracer released should be $q = {\mu} / {c^*(\vec{x}, t)}$. Thus a spot source from an important position and date should be weaker. A position and date with zero importance are not connected to the measurement and no source cou ld be proposed there as an explanation for it. This is already the backtracking approach developped in paragraph 4. \subsection*{Academic example} \def\ms{m~s$^{-1}$} We first show an illustration of Eulerian backtracking. A cheminee is assumed to have release 1kg of toxic gas at the top of a 200~m tower. The pollutant is detected at 150~km downwind of the detector. The flow is supposed to be hroiztontal with a wind assymptotic to $U_{top}=$10~\ms in the free atmosphere reducing to $U_{bot}=$7~\ms at the surface. The wind is prescribed as $u(z)=U_{bot}+(U_{top}-U_{bot})(1-\exp{4z/H})$ where $H=1$~km is the height of the planerary boundary layer. The eddy-mixing coefficient is also dependant on height as $K_z=\max{K_{bot},(K_{top} z/h}$. \begin{figure*} \centerline{\includegraphics[width=14cm]{acad1.eps}} \caption{Academic illustration of direct and reverse transport\label{eq:acad1}} \end{figure*} \begin{figure*} \centerline{\includegraphics[width=14cm]{acad2.eps}} \caption{Academic illustration of direct and reverse transport\label{eq:acad2}} \end{figure*} The transport equation is integrated on a domain maid of 400 points in the $x$ direction (one point each 500~m) with 24 layers 50~m thick on the vertical. The time step is of 25~s. The time of injection of the tracer is 10 minutes (ie 24 timesteps). The measuremts are also averged on 10 minutes. Results are illustrated in Fig... The upper panel shows the direct plume at timestep ... \afaire{Les deux figures montrent l'evolution d'un panache depuis la surface en x=25 (panneau du haut) et le retro-transport depuis deux mesures effectuees en x=175. Les panneaux 2 et 3 correspondent respectivement a une retroinjection depuis la surface et depuis 700~m. Les deux figures correspondent a des valeurs differentes du coefficient de melange turbulent (plus faible pour la seconde).} The second and third panels show the retro-plumes for injection in the 4th and 14th layer respectively. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{Backtracking} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection*{Identifying spot sources} We shall call more specifically now 'backtracking' the quest, as an explanation for a positive tracer measurement, of an event source with no or weak extent in space and time. This question often arises from the atmospheric monitoring of pollutants. When a positive detection has an abnormal nature, the source is thought to be an event. It might be an industrial accident. We originally studied the detection of nuclear tests with sources that are, ideally, space time points. The backward reconstruction of individual trajectories has been popularized in the recent years. This use of Lagrangian calculations probably aims at avoiding the inversion of diffusion which is, as already mentioned, wrongly considered impossible. However, it is generally seen as a way of estimating qualitatively the origin of air masses, in order to interprete the composition measured at a station, or the origin of the pollution imported on a city \afaire{Faut-il citer quelque chose ici. Pas Vautard n'co qui utilisent justement les retro-trajectoires pour des prédictions quantitatives avec prise en compte du mélange de couche limite}. The position and date of the source are investigated in the neighbourhood of a single backtrajectory obtained out of a backward integration of the windfield departing from the detector. Often, in order to manage statistically the uncertainty of the method, a great number of backtrajectories are calculated with slightly different initial conditions. The source is then supposed to lie more probably in a region of space and time many backtrajectories go through. This approach amounts to saying that atmospheric dispersion is a non deterministic statistical process: two identical tracers released at the same place and date could be spread in different manners. This is not correct. Dispersion is deterministic, even with much turbulence. There are mainly two reasons, as explained in the next two sections, why the position and date of a source cannot be determined accurately. Firstly the identification of the source is biased by the limited validity of atmospheric models including our imperfect knowledge of the wind field and the limited validity of diffusion as a model for turbulence. Secondly, when handling a single measurement, i.e. a single datum, it is underdetermined. \subsection*{Single particles} Let's again focus on perfect tracers with the notations of paragraph 1 in order to clarify the statistical meaning of the exchange rate. Let's pick up one particle in the detector : where does it come from? The proportion in the detector $(D,t_d)$ of air particles coming from the source $(S,t_s)$ is : \begin{equation} \label{eq:proba} p_{d \leftarrow s} = \frac{M_e(S,t_s,D,t_d)}{M(D,t_d)} \end{equation} So, the probability that our particle was in volume $S$ at time $t_s$ is simply $p_{d \leftarrow s} = M(S,t_s)\varepsilon(S,t_s,D,t_d)$ or $p_{d \leftarrow s} = M(S,t_s)\cstar{D,td}{S,ts}$. The latter expression shows that the retroplume $c^*_{D,td}$ is the density of probability of the position of the particle at any moment in the past before it was sampled.\\ It is clear as well that the plume $c_{S,ts}$ is the density of probability of the future position of a particle departing from the source $(S,t_s)$. Let's now consider non perfect tracers. The interpretation of the direct plume as a density of probability is unaltered. The physical relevance of $c^*_{D,td}$ as such is more doubtfull because all tracer particles that have been sampled may have been released from the same place at the same moment. Finally, the calculation of a great number of Lagrangian backtrajectories just amounts to calculating $c^*_{D,td}$. The problem is that the statistical interpretation of the retroplume, valid for single particles, is false when considering macroscopic sources. \subsection*{Macroscopic sources} If the value observed for the measurement of a linear tracer is $\mu$, then the value $c^*(\vec{x},t)$ of the retroplume just indicates that the source could be in $(\vec{x},t)$ provided its intensity is exactly \begin{equation}\label{eq:eqk} K(\vec{x},t)=\mu / c^*(\vec{x},t) \end{equation} Our knowledge of $c^*$ is biased by the quality of our model, for instance the quality of wind data. The quality of the model includes as well the validity of diffusion as a representation of turbulence. Diffusion always undergoes the following shortage which is a source of confusion with the single particle problem. Because of the sharp frontier with integrated scales, turbulent motions may have macroscopic consequences. In other words, it is not always clear whether a source of tracer is really macroscopic. The second reason why we cannot find out exactly the position of the source is that, even if we perfectly know $c^*$, according to \eq{eqk}, it can really be anywhere. This uncertainty is not of a statistical nature, it is just an underdetermination. An instantaneous spot source corresponds to five parameters for the position, date and intensity. Five independent measurements are required, in theory, to rebuild it. We shall adress in a forthcoming paper how to handle multiple measurements. The quest for the source can be simplified in most practical situations. Firstly, it is often reasonable to look for a source at ground or sea level. Secondly, the position is generally more interseting than the date which can be ignored. It is even possible to admit an undetermined duration. Suppose the source in $\vec{x}$ had a rate of release $D(t)$ per unit time (the source function is $D(t) \sigma_{\vec{x},t}$). The total release is $Q=\int D(t) dt$ , and the measurement $\mu=\int c^*(\vec{x},t) D(t) dt$. As D is a positive function $\mu \leq \max_t c^*(\vec{x},t) \int D(t) dt$, or : \begin{equation}\label{thr1} Q \geq \frac{\mu}{\max_t c^*(\vec{x},t)} = K_{min}(\vec{x}) \end{equation} Hence, for each possible position, we obtain a minimum intensity of a source there, no matter when, compatible with the measurement. This information will be compared to the positions of industrial installations. \subsection*{ETEX} On october 1994, between the 23d, 16 UT and the 24th, 4UT, 240 kg of an inert tracer, called pm ch, were released from a ground position in Brittany, France. Time averaged concentration measurements were delivered each third hour by 168 detectors all over Europe. \section{Concluding remarks} The next step towards the coherent interpretation of several measurements is of a totally different, combinatorial nature and will not be described here. Reparler des calculs de capacite de detection. Comment aurait-on traite ca en partant directement d'une approche adjointe? \appendix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \section{General derivation of the adjoint transport equation} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \def\cout{J} \afaire{Rajouter les termes de diffusion verticale et de puits ou source lineaire pour traiter exactement l'equation du texte principal} Let us apply the adjoint method to \eq{direct} using as a cost function the measure \begin{equation} \cout=\intespacetemps \rho \mu \ex \ \dv dt \end{equation} We will consider a variation of $\cout$ linked to a variation of the source fonction $d\source$. Since \eq{directb} is linear, it also applies for the time evolution of $\delta \ex$ with source term $\delta \source$ so that \begin{equation} \delta \cout=\intespacetemps \rho \mu \delta\ex \ \dv dt \end{equation} We will say that we have determined the sensitivity of measurement to the source function if we are able to transform the above expression in one of the form \begin{equation}\label{eq:sensitivity} \delta \cout=\intespacetemps \gamma \delta \source \ \dv dt \end{equation} Note that this expression is just a generallisation to a continous space of the sensitivity \begin{equation} \delta \cout=\sum_i \gamma[(\x,t)_i] \delta q[(\x,t)_i] \end{equation} which would be obtained for instance for a discretized space time domain. To arrive to \eq{sensitivity}, a very general classical approach (Lagrangian multipliers?) consists in introducing \eq{directb} as \begin{eqnarray}\label{eq:adjun} \delta \cout&=&\intespacetemps \rho \mu \delta\ex \ \dv dt\\ &-&\intespacetemps \rho \exr \depb{\dt{(\delta \ex)}+\vent.\grad(\delta\ex) -\delta \source}\ \dv dt \nonumber \end{eqnarray} ($\exr$ to be defined). By exactly the same algebra as previously consisting in integration by part taking into account the continuity equation, one easily obtain \begin{eqnarray}\label{eq:adjdeux} \delta \cout&=&\intespacetemps \rho \mu \delta\ex\ \dv dt \nonumber \\ &+&\intespacetemps \rho \delta \ex \depb{\dt{\exr}+\vent.\grad\exr} \ \dv dt \nonumber \\ &+&\intespacetemps \rho\exr \delta \source \ \dv dt \end{eqnarray} By taking for $\exr$ the solution of \eq{reciprob}, this equation gives the sensitivity as \begin{equation} \delta \cout=\intespacetemps \rho \exr \delta \source \ \dv dt \end{equation} By denoting $<\ex,\exr>$ the scalar product $\int \rho\ex\exr$, and an operator $L$ translating \eq{directb} into $L(\ex)=0$, what is done to pass from \eq{adjun} to \eq{adjdeux} is to replace the operator $L$ by its adjoint defined by the following relationship~: $=<\ex,L^*\exr>$ for any distribution $\ex$ and $\exr$. For a differential equation, the identity means that each term must be transformed by an integration by part. In general, the algebra can be more complex involving integral on the boundaries. The complete adjoint equation for the transport problem with boundary terms is given in appendix {\bf A ECRIRE}. Non linear equations must also be linearized before applying those technics. Note to finish that we chose here a linear cost fonction $\cout$. For a nonlinear cost-function, such as the root-mean-square error computed with respect to observation, the fonction must first be differentiated with respect to concentration. The results plays the role of source term in the adjoint equation. \bibliography{/d2/hourdin/tex/biblio/fred} \end{document}