\section{Retro-transport et transport adjoint}

Cette question du rétro-transport nous a donc été posée dans un cadre
tout militaire.
Il s'agissait d'évaluer la capacité d'un réseau de station
mesurant la radioactivité atmosphérique à détecter et si possible localiser
des essais nucléaires.

Nous avions donc une source, relativement ponctuelle en espace et en temps
et des détecteurs.
Pour la capacité de détection, la question directe peut être résolue de la 
façon suivante~: en un point de la planète, on injecte la quantité de
radio-élément correspondant à un essai nucléaire typique (la spécification
de détection du réseau TICE est de détecter des essais de 1~kt équivalent TNT).
On transporte ce radio-élément sur quinze jours. Sur ces quinzes jours, on
regarde si, à une des stations du réseau, la concentration en radio-élément
dépasse le seuil des détecteurs.
Reste à effectuer ce calcul à partir de tous les points du globe (les points
d'un maillage par exemple) et pour un ensemble statistiquement représentatif
de situations météorologique.

Il est beaucoup plus efficace de traiter ce problème en mode rétro-transport.
Dans ce cas particulier, la symmétrie est complète. On peut en fait injecter
la quantité de rétro-élément à la station, inverser le sens du temps dans
le modèle. Les points de la planètes auquels un essai aurait été détecté
dans les 15 jours précédents sont ceux ou la concentration du
rétro-radio-élement dépasse le seuil de détection.
Dans ce cas précis, le rapport de coût numérique est le rapport entre
le nombre de stations et le nombre de localisations testées.
Pour un réseau d'une cinquantaine de stations et un maillage avec une
résolution de quelques centaines de kilomètres, disons 10$^4$ points,
le rapport est de l'ordre de 200.

Cette symmétrie du transport atmosphérique ce comprend aisément en
introudisant un coefficient d'échange entre source et détecteur.



\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\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\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\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\massex#1#2#3#4{{\cal M}^{ex}(#1,#2,#3,#4)}
\def\masse#1#2{{\cal M}(#1,#2)}





\def\ex{c}
\def\exr{{c^*}}
\def\ex#1#2#3#4{\chi(#1,#2,#3,#4)}



\def\afaire#1{{\bf #1}}
\maketitle


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Réciprocité du transport atmosphérique et rétro-transport}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\def\cc#1#2{{\overline{c}}_{#1}^{#2}}
\def\cstar#1#2{{\overline{c^*}}_{#1}^{#2}}


On introduit tout d'abord la symmétrie du transport atmosphérique
en s'intéressant à un traceur parfait (qui suit les trajectoires fluides
sans source ni puit) distribué uniformément à un instant $\ts$ dans un 
volume source $S$. On suppose que la détection du polluant consiste en
la mesure de la concentration moyenne de traceur dans un volume $D$
au temps $\td$. 
Pour une quantité totale injectée $Q$ (quantité extensive en kg, atomes, ...),
la concentration moyenne dans $D$ à $\td$ peut s'écrire 
$Q\times\ex{S}{\ts}{D}{\td}$ avec
\begin{equation}
\ex{S}{\ts}{D}{\td}=\frac{\massex{S}{\ts}{D}{\td}}{\masse{S}{\ts}\masse{D}{\td}}
\end{equation}
où $\masse{{\cal V}}{t}$ est la masse d'air contenue dans le volume $\cal V$ au
temps  $t$ et $\massex{S}{\ts}{D}{\td}$ est la masse d'air échangée entre
$(D,\td)$ et $(S,\ts)$.
La quantité $\ex{}{}{}{}$ qu'on appelera
{\em coefficient d'exchange} dans ce qui suit,  est clairement symmétrique
et peut être évalué soit comme la concentration (intensive)
moyenne dans $(D,\td)$
résultant dsu transport direct d'un traceur injecté en quantité unitaire dans
$(S,\ts)$, soit comme la
concentration en $(S,\ts)$ d'un rétro-traceur injecté en $(D,\td)$,
qu'on suit en remontant à rebourd dans le temps le long des trajectoires
fluides.


\begin{figure}
\centerline{
\includegraphics[width=6cm,clip]{\local/FIGURES/schema1.eps}
\includegraphics[width=10cm,clip]{\local/FIGURES/schema2.eps}
}
\caption{Illustration de la réciprocité du transport atmosphérique
\label{fg:schema}}
\end{figure}

On peut donner une illustration de cette réciprocité dans deux cas extrêmes.

La première illustration (paneau A de la \fg{schema})
s'apparente à une vision de type advection
des contours. Pour simplifier l'image,
on va s'intéresser à un écoulement bi-dimensionnel non divergent.
Dans la limite d'un écoulement non visqueux, on sait que l'advection va
se contenter de déformer (et déplacer) la surface $S$ contenant initialement
le traceur.
Supposons que la nouvelle surface obtenue au temps $\td$ après advection,
intersecte $D$.
Cette intersection contient tout l'air qui venait de $S$ et qui se trouve
actuellent dans $D$. De façon symmétrique, si on remonte depuis $(D,\td)$
le long des trajectoires
fluides, les points de l'intersection doivent revenir dans $S$ alors que les
autres points du volume $D$ vont se disperser autour (par exemple sous
forme d'un filament).


Le second cas est celui d'une couche limite nocturne de 500~m
dans laquelle on injecte un polluant. A midi, on suppose
qu'une couche limite convective se développe brassant complètement l'air
sur 2~km. Si on détecte la polution en surface après le brassage, elle
sera 4 fois plus faible qu'avant.
De la même $façon$, si on marque l'air contenu dans les 500 premiers mètres
après le brassage et qu'on remonte les trajectoires individuelles des
particules fluides, ces particules avant le brassages proviennent
avec une équiprobabilité des deux premiers kilomètres. Le rétro-polluant
subit donc exactement la même dillution avant le brassage que le polluant
direct après. On retrouve bien que, pour une même injection, polluant direct
dans $D$ et rétro-polluant dans $S$ ont la même concentration.

\subsection{Extension à des sources et puits linéaires}

Définie de cette façon, la réciprocité s'étend facilement à des puits et
sources linéaires.
Dans ce cas, l'échange n'est plus régit uniquement par le taux d'échange d'air
comme dans le cas conservatif.
Il faut tenir compte de la création
ou de la destruction de traceur le long des trajectoires.

Dans le cas d'un radioélément avec un taux de décroissance 
$\lambda$,
si le même taux de décroissance est appliqué pour les transports direct
et inverse, la même concentration 
$Q\ex{S}{\ts}{D}{\td}exp\depb{-\lambda\dep{\td-\ts}}$
sera obtenue lors de la mesure en $(D,\td)$ d'un radio-élément injecté en
$(S,\ts)$ ou en mesurant en $(S,\ts)$ la concentration d'un rétro-radio-élément
injecté en $(D,\td)$.
Une décroissance radioactive dans le futur devient une décroissance
radioactive dans le futur.

Pour le problème concrêt qui nous avait été posé dans un premier temps,
à savoir l'évaluation de la capacité de détection du réseau TICE, 
la solution numérique du problème direct aurait consisté à injecter
10$^{15}$~Bq de Xenon (rejet supsposé d'un tir de okt équivalent TNT)
indépendamment dans chaque maille du modèle et à regarder au cours du
temps toutes les
stations $D_i$ pour lequelles on aurait dépassé un certain seuil
(le seuil de détection du réseau Xenon est typiquement fixé à \afaire{Mettre
au propre les aspects CTBT}).
Dans l'approche retro, il suffisait d'effectuer le tir aux stations,
puis de remonter dans le temps en appliquant la même décroissance
radioactive. Les lieux détecté étaient tous les points du globe pour
lesquels la concentration des rétro-panaches ainsi obtenus dépassaient 
le même seuil de détection.


Ce résultat s'étend en fait facilement à n'importe quel puit ou source
linéaire, pour lequel le taux de croissance ou décroissance peut
varier dans l'espace et dans le temps $\lambda(\x,t)$ (réaction
chimique avec un composant dominant pas ou très peu affecté par la réaction
en question, paramétrisations simples du lessivage par les pluies, ...).
Dans ce cas, on peut utiliser un coefficient d'échange généralisé,
défini colmme la concentration au niveau du détecteur
d'un traceur (non conservatif) unitaire réparti uniformément dans $S$ à $\ts$
ou, de façon équivalente, comme la concentration en $(S,\ts)$ d'un rétro-traceur
unitaire injecté en $(D,\td)$,
suivi en remontant le long des trajectoires d'air et en lui appliquant la
même source où le même puit que dans le calcul direct.
Ce coefficient généralisé s'écrit
\begin{eqnarray}\label{eq:noncons}
&&\ex{S}{\ts}{D}{\td}=\frac{1}{\masse{S}{\ts}\masse{D}{\td}}\times\\
&&\int_{\gamma \mbox{\ in\ } \Gamma_{S,D}}
{\exp\depb{-\int_{\ts}^{\td}{\lambda\dep{\gamma,t} dt}
}\rho(\x,\ts) d\Omega_S}
\nonumber
\end{eqnarray}
où $\Gamma_{S,D}$ est l'ensemble des trajectoires qui ont leur origine
dans $S$ au temps $\ts$ et leur extrémité en $D$ à $\td$, 
$d\Omega_S$ est un volume élémentaire dans $S$, à l'origine de la trajectoire
$\gamma$ et 
$\lambda\dep{\gamma,t}$ est la valeur de $\lambda$ au temps $t$ le long de
$\gamma$. Comme l'intégrale porte sur des trajectoires reliant $(S,\ts)$ et
$(D,\td)$, l'expression ci-dessus est inchangée si l'élément de masse
$\rho(\x,\ts)d\Omega_S$ dans le volume source est remplacé par
$\rho(\x,\td)d\Omega_D$ au niveau du détecteur. Pour $\lambda=0$,
l'intégrale se réduit à $\massex{S}{\ts}{D}{\td}$.

Le coefficient d'échange ainsi défini peut se calculer soit par une intégration
vers l'avant de l'équation d'advection 
\begin{equation}
\label{eq:direct}
\frac{\partial c}{\partial t}+\V.\grad\ c+\lambda c = \sigma
\end{equation}
où $c\dep{\x,t}$ est la concentration massique du traceur et
$\sigma\dep{\x,t}$ est la distribution de la source, égale à
$\delta(t-\ts)/\masse{S}{\ts}$ à l'intérieur du volume $S$,
et à 0 à l'extérieur, avec la condition suppléementaire que $c\dep{\x,t}=0$
à un instant $t_i<\ts$, et qu'il n'y a pas d'apport de traceur au travers des
frontières du domaine $\Omega$ considéré (condition satisfaite automatiquement
pour un domaine global comme celui de LMDZ).

De façon symmétrique, le coefficient d'échange peut être calculé en intégrant
à rebourd dans le temps l'équation
\begin{equation}
\label{eq:retro}
-\frac{\partial c^*}{\partial t}-\V.\grad\ c^*+\lambda c^* = \mu
\end{equation}
$c^*\dep{\x,t}$ est à nouveau une concentration massique d'un traceur
qui sera appelé {\em retro-traceur} par la suite.
$\mu(\x, t)$, la distribution de la mesure, vaut
$\delta(t-\td)/\masse{D}{\td}$ à l'intérieur de $D$ et 0 en dehors, avec
la condition supplémentaire que $c^*\dep{\x,t}=0$ pour un temps
$t_f>\td$, et qu'il n'y a pas non plus d'apport de rétro-traceur le long
des frontières du domaine $\Omega$.

La reciprocité du transport atmosphérique peut alors s'écrire formellement comme
\begin{equation}\label{eq:princip1}
\ex{S}{\ts}{D}{\td}=\int_{\Omega\times\tau} \rho \mu c \dxdt
                   =\int_{\Omega\times\tau} \rho \sigma c^* \dxdt
\end{equation}
où $\tau=[t_i,t_f]$ est le domaine temporel considéré.


Du fait de la linéarité du transport, les équations ci-dessus s'étendent
également à des émissions $\sigma$ et mesures $\mu$ non locales,
que ce soit dans le temps ou dans l'espace.

		   \def\L{{\cal L}}

Si on réécrit à présent les équations cid-essus sous forme de
relations entre sources et concentrations, $c=\L(\sigma)$ et $c^*=\L(\mu)$,
la réciprocité du transport atmosphérique se résume alors à
$\scal{\L(\sigma),\mu}=\scal{\sigma,\L^*(\mu)}$
ce qui établit, sur des bases purement physique, que les équations pour le
transport direct et rétro sont adjointes l'une de l'autre
(ou que les opérateurs $L$ et $L^*$ sont adjoints l'un de l'autre)
pour le produit scalaire pondéré par la masse de l'air
\begin{equation}\label{eq:scalarproduct}
 \scal{\phi,\psi}=\int_{\Omega\times\tau} \rho \phi\psi \dxdt
\end{equation}


\subsection{Transport adjoint}




\def\cout{J}

Nous utilisons ici l'approche adjointe pour aboutir par un autre chemin
aux résultats de  la section précédente.

La mesure obtenue dans le détecteur $D$, peut, quelque soit
la distribution $c$ du traceur, s'écrire sous la forme
\begin{equation}\label{eq:cout}
\cout=\int_{\Omega\times\tau} \rho \mu c  \ \dxdt
\end{equation}
Cette mesure dépend donc linéairement de la concentration $c$,
qui elle-même dépend linéairement de l'émission $\sigma$ ainsi
que du flux de traceur entrant, le cas échéant.
Ce flux entrant a lieu le long de la {\em frontière entrante}
$\domegai$ du domaine $\Omega$ considéré,
\ie\   la partie de la frontière $\domega$
 le long de laquelle la vitesse $\V$ est dirigée vers l'intérieur du
 domaine
$\Omega$ ($\V.\n<0$ où $\n$ est le vecteur normal sortant).

La méthode adjointe fournit une approche générale pour expliciter
le lien entre une observable quelconque (ici $J$) et n'importe quel
paramètre d'entrée (ici la source et l'entrée de traceurs aux
frontières du domaine).
Voici comment se décline la méthode.

L'\eq{direct} est introduite dans  
L'\eq{cout} qui est réécrit
\begin{eqnarray}
\cout&=&\int_{\Omega\times\tau} \rho \mu c \  \dxdt\nonumber
\\\label{eq:adjun}
&-&\int_{\Omega\times\tau} \rho c^* \depb{\dt{c}+\V.\grad c
+\lambda c-\sigma}\ \dxdt
\end{eqnarray}
où $c^*\dep{\x,t}$ (qui est fondamentalement un multiplicateur de
Lagrange) est à déterminer.

Si on transforme par intégration par partie la partie advective de
l'\eq{adjun}:
\begin{eqnarray}\label{eq:adjdeux}
I&=&\int_{\Omega\times\tau} \rho c^* \depb{\dt{c}+\V.\grad c} \dxdt\\
 &=&\int_{\Omega}\depb{\rho c^* c}_{t_i}^{t_f} \dx\nonumber
+\int_\tau \depb{\rho\V c^* c.\n}_{\domega} dt\\
&-& \int_{\Omega\times\tau} c \depb{\dt{\rho c^*}+\div{\rho\V c^*}}\dxdt
\end{eqnarray}
on reconnait la forme conservative de l'équation de continuité pour 
$c^*$, qui peut être transformée en forme advective en utilisant 
l'équation de continuité pour l'air~:
\begin{equation}\label{eq:conserv}
\dt{\rho}+\div{\rho \V}=0
\end{equation}
Après réarangement des termes, on obtient
\begin{eqnarray}\label{eq:adjtrois}
\cout&=&\int_{\Omega\times\tau} \rho c^* \sigma \ \dxdt\nonumber\\
 &-&\int_{\Omega}\depb{\rho c^* c}_{t_i}^{t_f} \dx
-\int_\tau \depb{\rho\V c^*  c.\n}_{\domega} dt\nonumber\\
&+&\int_{\Omega\times\tau} \rho c \depb{\dt{c^*}+\V.\grad c^*-\lambda c^*+\mu}
\ \dxdt
\end{eqnarray}
Si on prend pour $c^*$ la solution de l'\eq{retro}, avec la condition
que $c^*=0$ à l'instant $t_f$ ainsi que le long de la frontière sortante
$\domegao$ ($\V.\n>0$), on obtient finaelement
\begin{equation}\label{eq:adjquatre}
\cout=
 \int_{\Omega} {\rho c^* c}_{|_{t_i}} \dx
-\int_\tau \depb{\rho\V c^*  c.\n}_{{\domega}_i} dt
+\int_{\Omega\times\tau} \rho  c^*  \sigma \ \dxdt
\end{equation}
qui explicite $\cout$ comme fonction des conditions initiales
$c$ at $t_i$, du flux rentrant de traceur à la frontière
$\rho\V c.\n_{|\domegai}$ et de l'émission de traceur $\sigma$.
Les facteurs multiplicatifs associés dépendent linéairement de la
variable adjointe $c^*$, qui s'avère être identique au 
{\em retro-traceur} introduit précédement.

Dans le cas plus général où soit l'équation d'évolution soit l'observable
sont non linéaires, il faut considérer les perturbations de premier ordre
$\delta c$, $\delta\sigma$ et $\delta J$ au lieu de
$c$, $\sigma$ et $\cout$ \cite[se reporter par exemple à ][]{Tala:87}.
Dans ce cas, l'analogue de l'\eq{adjquatre} fournit la sensibilité
de $J$ par rapport aux paramètres d'entrée.
Ces sensibilités dépendent alors encore linéairement de la variable
adjointe.

Si on revient au problème qui nous intéresse, avec un traceur injecté
au temps $\ts$ (avec $c=0$ avant $\ts$) et pas d'apport de traceur
par les frontières,
l'\eq{adjquatre} se réduit à
\begin{equation}
\cout=\int_{\Omega\times\tau}\rho\mu c\dxdt=
            \int_{\Omega\times\tau}\rho \sigma c^*\dxdt
\end{equation}
ce qui complète la preuve mathématique du principe de réciprocité 
(\ref{eq:princip1})
exposé plus haut à partir de considérations phyiques.

A noter qu'avec la même algèbre, on montre que, pour un traceur conservatif
($\lambda=0$ et $\depb{\rho\V c.\n}_{{\domega}_i}=0$),
et pour des instants compris strictement entre émission et mesure,
\begin{equation}
\frac{d}{dt} \int_{\Omega} \rho c c^* \dx=0
\end{equation}
A tout instant $t$, $\int_{\Omega} \rho c c^* \dxdt$ est en fait une
évaluation de la mesure à partir de la distribution $c$ de traceur
à l'instant $t$ et de la distribution ($c^*$) au même instant de l'air qui sera
échantillonné plus tard au niveau du détecteur.

\def\scalg#1{\left[#1\right]}
On voit que la partie advective de l'équation
\begin{equation}
\frac{\partial c}{\partial t}+\V.\grad\ c=0
\end{equation}
est équvalente à son propre adjoint pour le produit scalaire pondéré
par la masse d'air ~(\ref{eq:scalarproduct}).
Cette propriété découle directement de la conservation de la
masse d'air ou de traceur.
\cite{Tala:87} avait remarqué que, si uen équation dévolution linéaire
conserve un scalar produit dans le temps, alors l'équation directe
et adjointe sont identiques. Ce résultat ce montre en écrivant
l'équation directe sous la forme
\begin{equation}\label{eq:directg}
\frac{d c }{dt}=Lc
\end{equation}
On écrit de même l'équation adjointe pour un produit scalaire
$\scalg{ , }$ sous la forme
\begin{equation}\label{eq:adjointg}
-\frac{d c^*}{dt}=L^*c^*
\end{equation}
Supposons que le produit scalaire est conservé au cours du temps
par l'\eq{directg}, ce qui veut dire que pour toute solution $c$
de l'équation \eq{directg}, la quantité
$\scalg{c,c}$ est constante dans le temps, alors on peut écrire
\begin{eqnarray}
0
&=&\frac{d}{dt}\scalg{c,c}\\
&=&\scalg{\frac{dc}{dt},c}+\scalg{c,\frac{dc}{dt}}\\
&=&\scalg{Lc,c}+\scalg{c,Lc}\\
&=&\scalg{\dep{L+L^*}c,c}
\end{eqnarray}
Ceci étant vrai quelque soit $c$, impliquer qu'on a $L+L^*=0$,
et que les équations \ref{eq:directg} and \ref{eq:adjointg} sont
donc identiques.

Dans le cas de l'advection pure, le la concentration massique du traceur $c$
est conservée pour n'importe quel élément de masse
$dm=\rho\dx$. De ce fait, la quantité
$\int_\Omega\rho c^2\dx$  est conservée dans le temps.
On obtient donc l'identité des équations directe et adjointe du transport
 pour le produit scalaire pondéré par la masse d'air
comme un cas particulier du résultat de \cite{Tala:87}.

Dans les dérivations algébriques présentées précédemment,
c'est la présence de la densité de l'air $\rho$ dans la
seconde intégrale de l'\eq{adjun} qui permet de tirer avantage de la conservation
de la masse
(\ref{eq:conserv}) pour obtenir la symmétrie exacte entre
les équations \ref{eq:retro} et \ref{eq:direct}.




