\section{Rétro-transport et transport adjoint\label{sec:transport}} Cette question du rétro-transport nous a donc été posée dans un cadre militaire. Il s'agissait d'évaluer la capacité d'un réseau de stations 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. La capacité de détection peut être évaluée avec une méthode directe 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 (l'objectif retenu par le TICE était de détecter des essais de 1~kt équivalent TNT au plus 15 jours après le tir partout sur le globe). On transporte ce radio-élément sur quinze jours en regardant si à une des stations du réseau, la concentration en radio-élément dépasse le seuil de détection des stations de mesure au cours de cette période. Reste à effectuer ce calcul à partir de tous les points du globe (les points d'un maillage globale par exemple) et pour un ensemble statistiquement représentatif de situations météorologiques. Il est beaucoup plus efficace de traiter ce problème en mode rétro-transport. Dans ce cas particulier, la symétrie est complète comme on le montre ci-dessous. On peut en fait injecter la quantité de radio-élément à la station puis inverser le sens du temps dans le modèle. Les points de la planètes auxquels un essai aurait été détecté dans les 15 jours précédents sont ceux où la concentration du rétro-radio-élément dépasse le seuil de détection. Dans ce cas précis, le rapport de coût numérique entre le calcul direct et rétro 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. On commence ici par introduire la symétrie du transport atmosphérique dans le cas d'une source et d'un détecteur ponctuels (comme dans le cas des calculs effectués pour le TICE) à partir de considérations physiques. On montre ensuite comment on peut retrouver la même propriété de symétrie de façon mathématique à partir de l'approche adjointe. \def\cout{{\cal M}} \def\A{{S}} \def\B{{D}} \def\V{{\bf v}} \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\sec#1{Section~\ref{sec:#1}} \def\S{S} \def\D{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{{\bf 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{{m}^{ex}(#1,#2,#3,#4)} \def\masse#1#2{m(#1,#2)} \def\ex{c} \def\exr{{c^*}} \def\ex#1#2#3#4{\chi(#1,#2,#3,#4)} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \subsection{Symétrie temporelle du transport atmosphérique} %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% \def\cc#1#2{{\overline{c}}_{#1}^{#2}} \def\cstar#1#2{{\overline{c^*}}_{#1}^{#2}} On introduit tout d'abord la symétrie du transport atmosphérique en s'intéressant à un traceur parfait (qui suit les trajectoires fluides sans source ni puits) distribué uniformément à un instant $\ts$ dans un volume source $S$. On suppose que la détection 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 massique moyenne $\cout$ (pour mesure) de traceur dans $D$ à $\td$ peut s'écrire \begin{equation}\label{eq:recipro1} \cout=Q\ \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)$. L'expression \eqb{recipro1} est clairement symétrique et la mesure $\cout$ peut être évaluée soit comme la concentration (intensive) moyenne dans $(D,\td)$ résultant du transport direct après injection d'une quantité $Q$ de traceur en $(S,\ts)$, soit comme la concentration en $(S,\ts)$ obtenue après injection d'une quantité $Q$ de rétro-traceur en $(D,\td)$ qu'on suit en remontant à rebours dans le temps le long des trajectoires fluides. \begin{figure} \centerline{ \includegraphics[width=7cm,clip]{\local/FIGURES/schemaA1g.eps} \includegraphics[width=7cm,clip]{\local/FIGURES/schemaA2g.eps} } \centerline{ \includegraphics[width=9cm,clip]{\local/FIGURES/schemaBg.eps} } \caption{Illustration de la réciprocité du transport atmosphérique. \label{fg:schema}} {\bf A}~: Le domaine $S$, initialement occupé par le traceur au temps $\ts$, est transformé par le mouvement en un filament au temps $\td$ (gris sombre). L'air contenu dans le volume $D$ se répartit alors entre de l'air provenant de $S$ (intersection de $D$ avec le filament, en noir) et de l'air sans traceur. Si on remonte le long des trajectoires atmosphériques à partir du volume $D$ au temps $\td$, les particules dans la partie noire proviennent de $S$, et les autres doivent se répartir à l'extérieur de $S$ (filament gris clair). {\bf B}~: Idéalisation de la dilution d'un polluant injecté le matin près de la surface, dans une couche limite nocturne peu épaisse (l'air pollué est montré en gris sombre), et mélangé en milieu de journée dans la couche limite convective pour donner en soirée un air moins pollué (gris clair). La concentration observée en soirée peut être obtenue soit en injectant le polluant le matin et en le mélangeant complètement sur la verticale à midi, soit par un calcul inverse, en injectant le polluant près de la surface le soir, et en le mesurant le matin précédent, en appliquant en milieu de journée le même mélange vertical que lors du calcul direct. \end{figure} On peut donner une illustration de cette réciprocité dans deux cas extrêmes. La première illustration (partie A de la \fig{schema}) s'apparente à une vision de type advection des contours. Pour simplifier l'image, on va s'intéresser à un écoulement bidimensionnel non divergent. Dans la limite d'un écoulement non visqueux, on sait que l'advection va se contente 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 actuellement dans $D$. De façon symé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, épaisse de 500~m pour fixer les idées, dans laquelle on injecte un polluant en surface. 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 pollution 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 brassage proviennent avec une équiprobabilité des deux premiers kilomètres. Le rétro-polluant subit donc exactement la même dilution 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 \begin{equation}\label{eq:timesym} \cout=Q\ exp\depb{-\lambda\dep{\td-\ts}}\frac{\massex{S}{\ts}{D}{\td}}{\masse{S}{\ts}\masse{D}{\td}} \end{equation} 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)$. A une décroissance radioactive dans le futur est associée la même décroissance radioactive dans le passé. Ce résultat s'étend en fait facilement à n'importe quel puits linéaire, pour lequel le taux de décroissance $\lambda(\x,t)$ peut varier dans l'espace et dans le temps (réaction chimique avec un composant très peu affecté par la réaction en question, paramétrisations simples du lessivage par les pluies, ...) \footnote{On parle dans la suite uniquement de puits linéaires, plus fréquents en pratique que les sources, mais les résultats s'appliquent bien sûr à une source linéaire associée à un coefficient $\lambda<0$.}. Dans ce cas, la concentration $\cout$ peut s'exprimer de façon générale comme \begin{equation}\label{eq:noncons} \cout= \frac{Q}{\masse{S}{\ts}\masse{D}{\td}}\ \int_{\gamma \mbox{\ in\ } \Gamma_{S,D}} {\exp\depb{-\int_{\ts}^{\td}{\lambda\dep{\gamma,t} dt} }\rho(\x,\ts) d\Omega_S} \end{equation} 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$.\footnote{ Cette écriture comme une intégrale sur l'espace des trajectoires atmosphériques est inspirée d'un travail effectué par ailleurs sur le calcul des échanges radiatifs dans l'atmosphère dans lequel la puissance nette échangée entre deux volumes de gaz peut être exprimée comme une intégrale sur l'espace des chemin optiques \cite[]{Cher:96,Dufr:05PNE}. } 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}$. On peut calculer $\cout$ en intégrant vers le futur l'équation d'advection \begin{equation} \label{eq:direct1} \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 à $Q\ \delta(t-\ts)/\masse{S}{\ts}$ à l'intérieur du volume $S$, et à 0 à l'extérieur, avec la condition supplémentaire 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é. De façon symétrique, on peut calculer $\cout$ en intégrant vers le passé l'équation de rétro-transport \begin{equation} \label{eq:retro1} -\frac{\partial c^*}{\partial t}-\V.\grad\ c^*+\lambda c^* = \mu \end{equation} où $c^*\dep{\x,t}$ est à nouveau une concentration massique d'un traceur qui sera appelé {\em rétro-traceur} par la suite. La distribution de la mesure $\mu(\x, t)$ 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 symétrie du transport atmosphérique peut alors s'écrire formellement comme \begin{equation}\label{eq:princip1} \cout=\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é. A noter qu'avec les notation choisies, ce sont les couples $(\sigma/Q,c/Q)$ pour le transport direct et $(\mu,c^*)$ pour le rétro-transport [en (kg$^{-1}$~s$^{-1}$, kg$^{-1}$)] qui jouent des rôles symétriques. 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 ci-dessus sous forme de relations entre sources et concentrations, $c\equiv \L(\sigma)$ et $c^*\equiv \L(\mu)$, la réciprocité du transport atmosphérique se résume alors à la relation $\scal{\L(\sigma),\mu}=\scal{\sigma,\L^*(\mu)}$. Ceci établit, sur des bases physiques, 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{Dérivation adjointe} On utilise ici l'approche adjointe pour aboutir par un autre chemin aux résultats de la section précédente. On en profite pour traiter un cas un peu plus général. On considère une mesure de concentration de la forme \begin{equation}\label{eq:cout} \cout=\int_{\Omega\times\tau} \rho \mu c \ \dxdt \end{equation} où $\mu(\x,t)$ est la fraction de traceur observée au point $\x$ et au temps $t$ par unité de masse du fluide qui transporte le traceur et par unité de temps. On considère la dépendance de $\cout$ par rapport à l'émission $\sigma(\x,t)$, qui peut être maintenant n'importe quelle distribution spatio-temporelle, ainsi que par rapport à la distribution initiale de traceur $c(\x,t_i)$ et à un possible apport latéral sur la {\em frontière entrante} $\domegai$ de $\Omega$, \ie\ le long de la partie de la frontière pour laquelle la vitesse est dirigée vers l'intérieur du domaine ($\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 $\cout$) et n'importe quel paramètre d'entrée (ici la source, la concentration initiale et l'apport de traceurs aux frontières du domaine). Voici comment se décline la méthode. L'\eq{direct1} est introduite dans l'expression de la mesure \ref{eq:cout}~: \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 reconnaît 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:conserv1} \dt{\rho}+\div{\rho \V}=0 \end{equation} Après réarrangement 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{retro1}, 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 finalement \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$ à $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 rétro-traceur introduit précédemment \footnote{ 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 \cout$ 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 $\cout$ par rapport aux paramètres d'entrée. Ces sensibilités dépendent alors encore linéairement de la variable adjointe.}. Dans le cas où l'atmosphère est initialement dépourvue de traceur et en l'absence d'apport latéral, l'expression \eqb{adjquatre} de la mesure \eqb{cout} se réduit au dernier terme, ce qui complète la démonstration mathématique de la relation de symétrie \eqb{princip1} établie précédemment sur la base de considérations cinématiques. 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, on a \begin{equation} \frac{d}{dt} \int_{\Omega} \rho c c^* \dx=0 \end{equation} A tout instant $t$, la mesure \begin{equation}\label{eq:dtmesure} \cout=\int_{\Omega} \rho c c^* \dx \end{equation} peut alors être évaluée à 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 par le détecteur. \subsection{Lien entre symétrie et conservation} \def\scalg#1{\left\{#1\right\}} Dans le cas d'un traceur conservatif ($\lambda=0$), la partie homogène de l'équation de transport \eqb{direct1} \begin{equation}\label{eq:direct111} \frac{\partial c}{\partial t}+\V.\grad\ c=0 \end{equation} est équivalente à 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} avaient en effet remarqué que, si une équation d'évolution linéaire \begin{equation}\label{eq:directg} \frac{d c }{dt}=Lc \end{equation} conserve un produit scalaire $\scalg{,}$ dans le temps -- on entend par là que pour toute solution $c$ de l'\eq{direct111} la quantité $\scalg{c,c}$ est conservée dans le temps --, alors cette équation d'évolution est identique à son équation adjointe pour le produit scalaire conservé. Cette dernière s'écrit \begin{equation}\label{eq:adjointg} -\frac{d c^*}{dt}=L^*c^* \end{equation} où $L^*$ est l'adjoint de $L$ par rapport à $\scalg{ , }$. La réciproque est également vraie. La preuve découle d'une succession de transformations élémentaires de la dérivée temporelle de la quantité conservée \begin{eqnarray} \frac{d}{dt}\scalg{c,c}\label{eq:cons1} =\scalg{\frac{dc}{dt},c}+\scalg{c,\frac{dc}{dt}} =\scalg{Lc,c}+\scalg{c,Lc} =\scalg{\dep{L+L^*}c,c} \end{eqnarray} Dans le cas où l'équation d'évolution conserve le produit scalaire $\scalg{,}$, on a donc $\scalg{\dep{L+L^*}c,c}=0$ quelque soit $c$. Ceci implique $L+L^*=0$ et les équations \ref{eq:directg} et \ref{eq:adjointg} sont donc identiques\footnote{A noter cependant que la première transformation nécessite une certaine régularité de la distribution $c$ et n'est donc pas valable dans un cadre mathématique absolument général.}. Dans le cas de l'advection pure, 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:conserv1}) pour obtenir la symétrie exacte entre les équations \ref{eq:retro1} et \ref{eq:direct1}.