\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}.